Bayesian bands and the comparison theory
The Gaussian limit where all three coincide, flat-direction formalism, honest-ignorance asymmetry.
Two chapters, two methods, and a suspicious pattern: whenever the model was linear and the data Gaussian, both gave the same answer. That is not a coincidence — it is a theorem, and this chapter proves it: near a well-constrained minimum, Hessian bands, replica ensembles, and the Bayesian posterior are three computations of one object. The interesting science is therefore entirely in the failure of the premises. We formalize the three regimes where the methods split — nonlinearity, multimodality, and flat directions — derive what each method does in each regime, and build the comparison theory that PPDF-15 puts to experimental test. Along the way we meet the asymmetry that decides our lab's verdict: the Bayesian answer is not always different, but it is the only one that announces when it is dominated by assumptions.
The Gaussian-limit theorem: one object, three routes
Assume a single minimum $\hat\theta$ with positive-definite curvature and data informative enough that the likelihood decays before nonlinearity matters. Expand the log-likelihood ($\chi^2 = -2\ln\mathcal L +$ const):
$$ \ln\mathcal L(\theta) = \ln\mathcal L(\hat\theta) - \tfrac{1}{4}\,\delta\theta^{\mathsf T} H\, \delta\theta + \mathcal O(\delta\theta^3) $$
Route 1 (Bayes). With a prior that is smooth and wide on the likelihood's scale (so it is locally constant), the posterior is $p(\theta\mid D) \propto \mathcal L(\theta)$, i.e. exactly $\mathcal N\!\big(\hat\theta,\, 2H^{-1}\big)$.
Route 2 (Hessian). §5.1's Derivation 2 read the same quadratic form and returned $C_\theta = 2H^{-1}$; the eigenvector sets at $\Delta\chi^2 = 1$, pushed through the master formula, reproduce $\sigma_{\mathcal O} = \sqrt{\nabla\mathcal O^{\mathsf T}\, 2H^{-1}\, \nabla\mathcal O}$ — the standard deviation of the linearized observable under Route 1's Gaussian. Equal, term by term.
Route 3 (replicas). In this regime the model is effectively linear across the likelihood's support, $T(\theta) \approx T(\hat\theta) + J\delta\theta$, so §5.2's equivalence theorem applies with $X \to J$: the ensemble is $\mathcal N\!\big(\hat\theta,\, (J^{\mathsf T}C^{-1}J)^{-1}\big)$. But for a (locally) linear model the exact Hessian is $H = 2J^{\mathsf T}C^{-1}J$ (§5.1's Gauss–Newton identity with vanishing curvature term), hence $(J^{\mathsf T}C^{-1}J)^{-1} = 2H^{-1}$. Equal again. $\blacksquare$
One Gaussian, $\mathcal N(\hat\theta,\, 2H^{-1})$; three algorithms. Corollary: method agreement is a measurement — it certifies you are in the regime; disagreement localizes which premise failed, and where. That is the entire design of our three-way experiment.
Regime (i): nonlinearity — skew
Keep one minimum but let $T(\theta)$ curve within the likelihood's support. The posterior $\propto e^{-\chi^2(\theta)/2}$ is now skewed — it leans toward the side where the model responds weakly (flat response = slow $\chi^2$ growth = fat tail). The three methods react differently. Bayes integrates the actual $\chi^2$, so the skew is captured by construction — given the model; no symmetry was ever assumed. Hessian replaces the tilted distribution by its osculating Gaussian: symmetric-wrong, with the asymmetric master formula a first-order patch that still cannot bend the quadratic. Replicas produce a skewed ensemble too, but §5.2's Derivation 3 showed it is the wrong skewed object: the pushforward of data noise through the fitting map, mean-biased by $\mathbb E[Q(\eta,\eta)] \neq 0$ and distributed differently from the posterior at second order (Costantini et al., arXiv:2404.10056). The first worked example makes all three concrete on one line of algebra.
How big is the effect? Make it dimensionless: with model curvature $K = \partial^2 T/\partial\theta^2$ and a $1\sigma$ parameter excursion $\delta\theta \sim \sigma_\theta$, the quadratic term competes with the linear one at the ratio $\kappa \equiv |K|\,\sigma_\theta / (2|J|)$. Skew corrections to all three methods enter at $\mathcal O(\kappa)$ and the replica mean-bias at $\mathcal O(\kappa^2)$ — so $\kappa$ is the single dial that tunes regime (i), and it grows with the parameter uncertainty, not the data noise alone. This is why the same parametrization can be "effectively linear" in a data-rich global fit ($\sigma_\theta$ small, $\kappa \ll 1$) and badly nonlinear in our 648-point laboratory, where weakly constrained directions have $\sigma_\theta$ of prior size and $\kappa \gtrsim 1$ with no expansion parameter left.
Regime (ii): multimodality — partition versus weight
Two minima with comparable $\chi^2$ — degenerate parametrization branches, sign-flip solutions — and the picture fractures. The Hessian is a strictly local instrument: it expands around whichever mode the optimizer delivered and reports that mode's width, silent about the other's existence. Replicas split: each replica's optimizer falls into the basin its initialization (and noise draw) assigns it, so the ensemble weights modes by basin geometry and init distribution — not by posterior mass $\int_{\rm mode}\mathcal L\,\pi$, which is what inference requires; a narrow-but-deep mode with a small basin gets systematically under-visited. Nested sampling (§4.2) was engineered for precisely this: its live-point population climbs all modes concurrently (cluster-tracking implementations follow them explicitly), and the evidence integral weights each by its actual mass. This regime matters less in our current fits than the other two — but degenerate branches are generic in flexible parametrizations, and only one of the three methods even has a concept of "relative weight of solutions".
Regime (iii): flat directions — the information matrix, formally
Define the information matrix $\;I \equiv J^{\mathsf T} C^{-1} J\;$ (the Fisher information of the Gaussian likelihood in the linearized model; Gauss–Newton curvature $= 2I$). A direction $v$ is flat iff $v \in \operatorname{null}(I)$, i.e. $Jv = 0$: moving along $v$ changes no prediction, hence no likelihood. Choose coordinates $\theta = (u, w)$ with $w$ spanning $\operatorname{null}(I)$, so $\mathcal L(D\mid\theta) = \mathcal L(D\mid u)$ exactly. With a factorized prior $\pi(\theta) = \pi_u(u)\,\pi_w(w)$, Bayes' rule gives
$$ p(u, w \mid D) \;=\; \frac{\mathcal L(D\mid u)\, \pi_u(u)\, \pi_w(w)} {\int \mathcal L(D\mid u')\,\pi_u(u')\,du' \int \pi_w(w')\,dw'} \;=\; p(u \mid D)\; \pi_w(w) $$
The posterior factorizes, and its $w$-marginal is exactly the prior: the analysis returns your assumptions, unmodified, on the combinations the data cannot see — which is the correct answer, honestly labelled. Meanwhile: the Hessian on $\operatorname{null}(I)$ has $\lambda = 0$ exactly (undefined covariance), or $\lambda \lessgtr 0$ at numerical resolution once residual noise meets model curvature — PPDF-14's eleven negative directions; and replicas, receiving zero likelihood pull along $w$, return the initialization-plus-procedure distribution of §5.2 — an answer with the form of an uncertainty and the content of a config file.
In regime (iii) no method conjures information that does not exist; the difference is disclosure. A Bayesian result ships with its prior, so "posterior width $\approx$ prior width along direction $w$" is a checkable, per-direction flag reading assumption-dominated here — and the scalar summary of those flags is exactly the Bayesian complexity: PPDF-12's $d_{\rm eff} = 13.7$ of 52 says the data genuinely measured $\sim$14 parameter combinations and the priors are load-bearing for the rest. A Hessian band after freezing, or a replica band in a flat direction, carries no flag: it renders as the same shaded region as a data-driven band, and nothing in the deliverable distinguishes them. The practical diagnostic this suggests — overlay prior width against posterior width, eigendirection by eigendirection — costs one plot and converts a philosophical dispute into a table. That, not any claim that Bayes "knows more", is what our lab's verdict rests on.
Coverage: the meta-criterion above all three
Is there a method-neutral referee? One: frequentist coverage — over repeated experiments, does the 68% band contain the truth 68% of the time? The catch is that coverage requires knowing the truth, so it is measurable only on planted-truth closure tests (§4.4). That is what PPDF-11 bought for the Bayesian leg at exactly our dataset-to-parametrization ratio: L0/L1 closure of the 52-parameter MSHT20 model on 648-point pseudo-data, max pull $2.63\sigma$, per-flavour $1\sigma$ coverage 65–100% against the 68% nominal — a calibration certificate, valid at that information ratio and claiming nothing beyond it. The corresponding certificates for Hessian and replicas at this ratio were attempted in PPDF-4/5/6 and came back negative (Hessian gluon pull $6\sigma$; replica bands init-shaped in flat directions). Honesty about the Bayesian route's own limits, stated plainly: (1) everything is conditional on the model — $\chi^2/N \approx 1$ is evidence of adequacy, not proof, and no posterior certifies its own model; (2) prior design is real scientific work — flat in which variables? how wide? — and regime (iii) makes the answer consequential. The Bayesian move does not eliminate that problem; it relocates it into the open, where prior–posterior comparison can measure it. The other two methods leave the same problem implicit in a freezing list or an initialization scale.
Worked example — one nonlinear measurement, three answers
Given: one measurement $y = 0.090 \pm 0.050$ of $T = \theta^2$, prior flat on $\theta \ge 0$. Best fit $\hat\theta = \sqrt{0.09} = 0.300$.
Hessian. $\chi^2 = (\theta^2 - y)^2/\sigma^2$, so $\chi^{2\prime\prime} = (12\theta^2 - 4y)/\sigma^2 = 8y/\sigma^2 = 288$ at the minimum. $\sigma_H = \sqrt{2/288} = 0.083$: symmetric interval $[0.217,\, 0.383]$.
Bayes. $p(\theta \mid y) \propto \exp[-(\theta^2 - 0.09)^2/0.005]$. Evaluate at the Hessian endpoints: $p(0.217)/p(0.383) = e^{-0.368}/e^{-0.643} = 1.32$ — the posterior is 32% taller at the lower endpoint: a heavier left shoulder toward $\theta = 0$, because the quadratic model responds feebly at small $\theta$ and the likelihood punishes slowly there. The correct 68% region is visibly asymmetric, roughly $[0.20,\, 0.37]$; only this route produces it.
Replicas. $\theta^{(k)} = \sqrt{y + \eta_k}$, $\eta \sim \mathcal N(0, \sigma^2)$. Jensen's inequality bites: $\mathbb E[\theta^{(k)}] \approx \sqrt{y}\,\big(1 - \sigma^2/8y^2\big) = 0.300 \times (1 - 0.039) = 0.288$ — the ensemble is skewed (good) but its centre is biased low by $0.14\sigma_H$ (bad), the $\mathbb E[Q(\eta,\eta)]$ term of §5.2 in the flesh, irreducible by more replicas.
Now slide $y \to 0$ (the measurement lands on the model's fold): Hessian curvature $8y/\sigma^2 \to 0$, band $\to \infty$ — undefined. Replicas: every $\eta_k \le 0$ replica ($\approx$ half) minimizes exactly at $\theta = 0$, the rest at $\sqrt{\eta_k}$: a point-mass-plus-tail ensemble, rms 0.141 but shaped by nothing physical. Bayes: $p \propto e^{-\theta^4/0.005}$, a smooth quartic-flattened peak with $\sigma_{\rm post} = \big[\Gamma(\tfrac34)/\Gamma(\tfrac14)\big]^{1/2}(2\sigma^2)^{1/4}\cdot(\dots) = 0.155$ — finite, smooth, and correct given the model. One toy, all three regimes-(i)-and-(iii) behaviours.
Given: one measurement $D = 2.00 \pm 0.10$ of $T = \theta_1 + \theta_2$; priors $\theta_i \sim \mathcal N(0, 2^2)$ independent. Rotate: $u = (\theta_1+\theta_2)/\sqrt2$, $w = (\theta_1-\theta_2)/\sqrt2$; then $T = \sqrt2\,u$, the prior stays $\mathcal N(0,4)\times\mathcal N(0,4)$, and $J = (\sqrt2,\, 0)$ in $(u,w)$ — so $I = J^{\mathsf T}C^{-1}J$ has $\operatorname{null}(I) = \operatorname{span}(w)$ exactly.
Posterior, explicitly (Derivation 2): $p(u, w \mid D) = p(u\mid D)\; \pi_w(w)$ with $p(u \mid D) \approx \mathcal N(1.414,\, 0.0707^2)$ (likelihood 28× narrower than prior — data-dominated) and $p(w\mid D) = \mathcal N(0, 2^2)$, the prior verbatim, width 2.0. The Bayesian deliverable states both numbers and flags the second as prior.
Hessian: $H = \operatorname{diag}(400,\, 0)$ in $(u,w)$; $C_\theta = 2H^{-1}$ does not exist. "Freeze $w$" quotes $\sigma_w = 0$: a delta prior in place of the honest 2.0.
Replicas: the gradient has no $w$-component, so $w^{(k)} = w_{\rm init}^{(k)}$ for every replica. Initialize $\mathcal N(0, 0.5^2)$: band 0.5 — under-covers the honest width by 4×. Initialize $\mathcal N(0, 8^2)$: over-covers by 4×. Same likelihood, same data, opposite verdicts — the "no rule" of PPDF-15, where real-fit flat combinations show MC/Bayes width ratios from 0.6 ($V$, under) to 8.1 ($\Sigma$) and 10.3 ($T_8$, over), with central offsets to $7\sigma$.