Nested sampling, derived
Volume-shrinkage statistics, the evidence estimator and its error budget, slice sampling at high dimension.
Section 4.1 ended with an integral: $Z = \int L(\theta)\,\pi(\theta)\,d\theta$ over 52 dimensions. A grid of just 10 points per axis would need $10^{52}$ likelihood calls; the observable universe does not contain enough silicon. Nested sampling — Skilling's 2004 algorithm, run in this lab through UltraNest — solves the problem by a change of perspective so clean it deserves to be derived rather than described: it converts the $d$-dimensional integral into a one-dimensional one over prior volume, then measures that volume statistically. This chapter derives the identity, the shrinkage statistics, the error budget, and the constrained-sampling machinery (slice sampling) that makes it work at $d = 52$ — ending with a nested-sampling run you can do by hand and a decoded line from our actual logs.
The change of variables: likelihood level sets
Define, for every likelihood level $\lambda$, the prior volume above that level:
$V$ is monotone decreasing: $V(0) = 1$ (all of the prior), $V(L_{\max}) = 0$. It is the cumulative distribution of likelihood under the prior, and it will absorb all 52 dimensions at once.
Write the likelihood at each point as an integral over levels ("layer cake"): $L(\theta) = \int_0^\infty \mathbf 1\big[\lambda < L(\theta)\big]\, d\lambda$. Substitute and swap the order of integration (Fubini — everything is positive):
$$ Z = \int \pi(\theta) \int_0^\infty \mathbf 1\big[L(\theta) > \lambda\big]\, d\lambda\, d\theta = \int_0^\infty \underbrace{\int_{L(\theta)>\lambda} \pi(\theta)\, d\theta}_{V(\lambda)}\; d\lambda = \int_0^\infty V(\lambda)\, d\lambda $$
Now use the geometric fact that the area under a monotone curve equals the area under its inverse: the region $\{(\lambda, V) : 0 < V < V(\lambda)\}$ can be sliced vertically ($\int V\,d\lambda$) or horizontally ($\int L\,dV$, where $L(V)$ is the inverse function — the likelihood level enclosing prior volume $V$). Hence
$$ Z \;=\; \int_0^1 L(V)\, dV $$
A one-dimensional integral of a monotone decreasing function. Every dimension of $\theta$-space has been hidden inside the map $V \mapsto L(V)$. $\blacksquare$
We cannot compute $L(V)$ analytically, but we can sample it. Maintain $N$ "live points" drawn from the prior subject to $L > \lambda_{\text{current}}$. Key fact: if $\theta$ is uniform (w.r.t. $\pi$) inside the current constrained region of volume $V_{\text{cur}}$, then its own enclosed volume $u = V(L(\theta))/V_{\text{cur}}$ is uniform on $[0,1]$ — the probability-integral transform. At each iteration we remove the worst live point (lowest $L$, hence largest enclosed volume): the new volume is $V_{\text{new}} = t\, V_{\text{cur}}$ where $t = \max(u_1, \dots, u_N)$. The largest of $N$ uniforms has density $p(t) = N t^{N-1}$, so:
$$ \mathbb E[\ln t] = \int_0^1 \ln t \; N t^{N-1} dt = \Big[t^N \ln t\Big]_0^1 - \int_0^1 t^{N-1} dt = -\frac1N, \qquad \operatorname{var}[\ln t] = \frac{1}{N^2} $$
(the variance from $\mathbb E[(\ln t)^2] = 2/N^2$ by the same integration by parts, twice). Successive shrinkages are independent, so after $k$ removals
$$ \ln V_k = \sum_{i=1}^k \ln t_i \;\approx\; -\frac{k}{N} \;\pm\; \frac{\sqrt k}{N} \qquad\Longrightarrow\qquad V_k \approx e^{-k/N} $$
The volume shrinks by a known statistical factor per step — geometrically, steadily, regardless of dimension. This is the engine: likelihoods are measured exactly, volumes are known probabilistically. $\blacksquare$
With the sequence $(V_k, L_k)$ in hand ($L_k$ = the discarded point's likelihood), the integral is a weighted sum (trapezoid in practice):
$$ Z \;\approx\; \sum_k L_k\, \Delta V_k, \qquad \Delta V_k = V_{k-1} - V_k \approx \frac{V_k}{N} $$
The dominant uncertainty is not in the $L_k$ (exact) but in the $\ln V_k$, which random-walk with standard deviation $\sqrt k / N$. The posterior bulk — where $L\,\Delta V$ peaks — sits at prior volume $\ln V_* \approx -H$, where $H = \int p \ln (p/\pi)\, d\theta = \mathbb E_{\text{post}}[\ln L] - \ln Z$ is the information (KL divergence from prior to posterior): the number of e-folds of compression the data demand. Reaching it takes $k_* \approx NH$ steps, by which time the accumulated volume jitter is $\sqrt{k_*}/N = \sqrt{H/N}$. Since the whole sum rides on those volumes,
$$ \operatorname{var}[\ln Z] \;\approx\; \frac{H}{N} $$
— more live points or a gentler compression, smaller error; and since cost per shrinkage step is $\mathcal O(1)$ likelihood calls at fixed efficiency, halving $\sigma(\ln Z)$ costs $4\times$ the compute. $\blacksquare$
The algorithm, stated plainly
(1) Draw $N$ live points from the prior; set $Z = 0$, $V_0 = 1$. (2) Find the worst live point: record
$(L_k, V_k \approx e^{-k/N})$, add $L_k \Delta V_k$ to $Z$. (3) Replace it with a fresh point drawn from the
prior subject to $L > L_k$ (the hard part — below). (4) Repeat until the termination criterion:
the largest possible remaining contribution, $L_{\max}^{\text{live}} \cdot V_k$, falls below a set fraction
(frac_remain) of the accumulated $Z$. (5) Posterior samples come free: every discarded point,
weighted by its share of the integral,
— exactly the weighted samples that §4.1's pushforward identity consumes. One run, both deliverables: the evidence for model comparison (§4.3) and the posterior for bands. That two-for-one is why this lab's entire pipeline is built on nested sampling rather than on MCMC plus a separate evidence estimator.
A single MCMC chain equilibrates within a mode; crossing a deep likelihood valley requires an exponentially unlikely excursion, so a second mode 30 units of $\ln L$ away may simply never be visited — and no diagnostic run within one mode will tell you. Nested sampling never crosses valleys: it approaches from below. The constraint $L > \lambda$ starts at $\lambda = 0$, where the live ensemble blankets the whole prior, and rises adiabatically; as the constrained region splits into islands, live points are already resident on each, and each island's volume is tracked by the same shrinkage statistics. Modes are not found by travel but by inheritance. (The failure mode is quantitative instead: if $N$ is too small an island can lose all its residents by unlucky discards — one reason our production runs use 300–1200 live points, and why UltraNest's reactive mode grew our population mid-run, below.)
Sampling under the constraint — the real engineering
Step (3) hides the entire difficulty: draw from the prior restricted to $L > \lambda$, where that region is a shrinking, curved, possibly disconnected sliver of a 52-dimensional box.
The naive method — draw from the full prior, keep if $L > \lambda$ — has acceptance rate exactly $V(\lambda)$, which the algorithm itself drives to $e^{-k/N}$: exponentially small by design. Even sampling from a bounding box around the live points fails geometrically: if the constrained region is a blob of radius $r_{\text{shell}}$ inside a bounding box of side $r_{\text{box}}$, the acceptance is the volume ratio
$$ \text{acc} \;\sim\; \Big(\frac{r_{\text{shell}}}{r_{\text{box}}}\Big)^{d} $$
At $d = 52$, a bounding box just $30\%$ too wide per axis costs $1.3^{52} \approx 8\times10^5$ — a million wasted likelihood calls per accepted point. Any method whose per-axis sloppiness multiplies is dead at high $d$; we need moves whose cost is additive in dimension. $\blacksquare$
Slice sampling is that method. To draw from a 1-D density under a constraint: from the current point,
place an interval of width $w$ around it; step out — extend the interval in units of $w$ until both ends
lie outside the constrained region; then sample uniformly inside, shrinking the interval toward the
current point on each rejection until a point satisfying $L > \lambda$ is found. Every "rejection" does useful
work (it shrinks the bracket), so a 1-D draw costs a handful of likelihood calls regardless of how awkwardly the
slice is placed. In $d$ dimensions, UltraNest applies this along a sequence of directions (random, or informed by
the live-point covariance); decorrelating from the start point needs $\mathcal O(d)$ such 1-D moves. Our runs use
nsteps = 156 $= 3 \times 52$ — three full direction sweeps per replacement, a standard safety factor:
too few steps leaves the new point correlated with a discarded one, which silently biases volumes (the shrinkage
derivation assumed independent uniform draws).
Population self-adjustment. UltraNest's reactive sampler monitors the evidence error budget as it runs; in our MSHT20 real-data fit it judged $\sqrt{H/N}$ too large mid-run and grew the population from 300 toward 1200 live points, exactly the $N$-scaling Derivation 3 prescribes. The observed cost showed up where the theory says it must: overall efficiency (accepted points per likelihood call) landed between 0.16% and 0.9% across our runs — each posterior sample at $d=52$ costs of order $156$ slice evaluations plus stepping-out overhead, i.e. hundreds of calls.
Reading the logs
A real status line from an UltraNest run in this lab (see the run notes attached to PPDF-12):
Decoded: iteration $= k$, the number of shrinkage steps, so $\ln V \approx -k/N$;
ncalls the total likelihood evaluations, giving efficiency $41500/25.9\times10^6 \approx 0.16\%$;
logz the running $\ln \sum L_k \Delta V_k$, already close to the final $-381.9 \pm 0.29$;
remainder_fraction the termination watchdog $L_{\max} V_k / Z$ — at 0.05% the not-yet-integrated
volume can shift $\ln Z$ by at most $\sim 5\times10^{-4}$, so the run is essentially done;
Lmin..Lmax the $\ln L$ span of the current live ensemble — 26 units here, the width of the bump in
the figure above. When Lmax plateaus while Lmin climbs toward it, the ensemble is
compressing onto the peak: the endgame.
Worked example 1 — nested sampling by hand
Geometry first: $L$ decreases in $\theta$, so the prior volume above level $L(\theta)$ is just $V = \theta$: the exact profile is $L(V) = e^{-V^2/0.08}$, and the exact evidence is $Z = \int_0^1 e^{-\theta^2/0.08} d\theta = \sigma\sqrt{\pi/2}\,\mathrm{erf}(1/\sigma\sqrt2) = 0.25066$, $\ln Z = -1.384$.
First shrinkage steps ($V_k = e^{-k/100}$, the expected trajectory): $k{=}1$: $V = 0.99005$, $L = e^{-0.99^2/0.08} = 4.8\times10^{-6}$, contribution $\approx 4\times10^{-8}$. $k{=}2$: $V = 0.98020$, $L = 6.1\times10^{-6}$. $k{=}3$: $V = 0.97045$, $L = 7.7\times10^{-6}$. Partial $Z$ after three steps: $\sim 1.6\times10^{-7}$ — we are integrating the likelihood desert, exactly as the bump figure predicts.
Checkpoints every 50 steps (trapezoid between rows; $\theta_k = V_k$):
$k{=}50$: $V=0.6065$, $L=0.0101$, running $Z = 0.0020$
$k{=}100$: $V=0.3679$, $L=0.1842$, $Z = 0.0252$
$k{=}150$: $V=0.2231$, $L=0.5367$, $Z = 0.0773$
$k{=}200$: $V=0.1353$, $L=0.7954$, $Z = 0.1358$
$k{=}250$: $V=0.0821$, $L=0.9192$, $Z = 0.1814$
$k{=}300$: $V=0.0498$, $L=0.9695$, $Z = 0.2119$
$k{=}350$: $V=0.0302$, $L=0.9887$, $Z = 0.2311$
$k{=}400$: $V=0.0183$, $L=0.9958$, $Z = 0.2429$; terminate and add the remaining slab
$\approx V_{400} \times 1 = 0.0183$.
Result: $Z \approx 0.2612$, $\ln Z \approx -1.343$, vs exact $-1.384$. The $+0.04$ discrepancy is the coarse 50-step trapezoid (a per-step sum lands within $10^{-3}$); the statistical error from Derivation 3 is $\sqrt{H/N}$ with $H = \mathbb E_{\text{post}}[\ln L] - \ln Z \approx -0.50 + 1.38 = 0.88$ nats, i.e. $\sigma(\ln Z) \approx 0.094$ — our deterministic-trajectory answer sits well inside one sigma of exact. Note the shape of the run: 400 steps for $-\ln V \approx 4$ e-folds at $N = 100$, with the bulk of $Z$ arriving between $k \approx 100$ and $250$, i.e. around $-\ln V \approx H \approx 0.9$–$2.5$. Everything scales from here: our real fits compress tens of e-folds at $N = 300$–$1200$, hence the tens of thousands of iterations in the log line above.
Worked example 2 — the error budget of a real run
Inputs: information $H \approx 40$ nats (the posterior occupies $\sim e^{-40}$ of the 52-dimensional prior box — a typical compression once $\sim$14 directions are constrained, §4.3), baseline population $N = 300$.
Prediction: $\sigma(\ln Z) \approx \sqrt{H/N} = \sqrt{40/300} = \sqrt{0.133} = 0.37$.
Reported: $\ln Z = -381.9 \pm 0.29$ (PPDF-12). The measured error is slightly better than the $N=300$ estimate — consistent with the reactive population growth toward 1200 live points through the compression phase ($\sqrt{40/1200} = 0.18$ would be the pure-1200 limit; a run that spends its early e-folds at 300 and its bulk at higher $N$ lands in between). Sanity rule used across our result cards: an evidence error budget of $\sim 0.3$, comfortably below $0.5$, is what makes verdicts like $\Delta\ln Z = -109$ (PPDF-9) hundreds of sigma — the sampling noise is never the story.