Monte-Carlo replicas in full
The parametric bootstrap, the linear-model equivalence theorem proved, where it breaks.
NNPDF made this method the modern default: instead of computing curvature, refit the noise. Generate many pseudo-datasets that fluctuate around the real one exactly as the experiment's covariance says data should fluctuate, fit each one, and read uncertainties off the spread of the fitted ensemble. No Hessian, no $\Delta\chi^2$ surfaces, embarrassingly parallel. This chapter builds the method completely: the bootstrap logic, a full proof of the theorem that makes it exact for linear models — including the demonstration that it then coincides with the Bayesian posterior — the precise way nonlinearity breaks it, the practical machinery of early stopping and post-fit selection with the noise-layer arithmetic our own runs (PPDF-5, PPDF-13) verify to three digits, and the failure mode that flat directions inflict on it — different in kind from the Hessian's, and quieter.
The parametric bootstrap
The experiment reports data $D \in \mathbb{R}^N$ with covariance $C$. The method manufactures $K$ replica datasets,
fits each one, $\theta^{(k)} = \operatorname{argmin}_\theta\, \chi^2\big(\theta;\, D^{(k)}\big)$, and takes ensemble statistics: the central PDF is the replica mean, the band the replica standard deviation, and any observable is evaluated per replica and summarized the same way. The logic is frequentist: if data fluctuate around truth with covariance $C$, then refitting datasets that fluctuate around our data with the same $C$ should map out how the fit responds to statistically legitimate wiggles. Whether that response distribution is the uncertainty you want is exactly the question the rest of this chapter answers.
Let the model be linear, $T(\theta) = X\theta$ with fixed $N\times p$ matrix $X$. The replica-$k$ objective is $\chi^2_k(\theta) = (X\theta - D - \eta^{(k)})^{\mathsf T} C^{-1} (X\theta - D - \eta^{(k)})$. Setting the gradient to zero, $2X^{\mathsf T}C^{-1}(X\theta - D - \eta^{(k)}) = 0$, and writing $M \equiv X^{\mathsf T} C^{-1} X$ (invertible if the design has full rank — flag that caveat now):
$$ \theta^{(k)} = M^{-1} X^{\mathsf T} C^{-1}\big(D + \eta^{(k)}\big) = \hat\theta + M^{-1} X^{\mathsf T} C^{-1}\, \eta^{(k)} $$
where $\hat\theta = M^{-1}X^{\mathsf T}C^{-1}D$ is the fit to the unfluctuated data. The map $\eta \mapsto \theta$ is linear, so $\theta^{(k)}$ is Gaussian with mean $\hat\theta$ and covariance
$$ \operatorname{Cov}\big[\theta^{(k)}\big] = M^{-1} X^{\mathsf T} C^{-1}\, \underbrace{\mathbb E\big[\eta\,\eta^{\mathsf T}\big]}_{=\,C}\, C^{-1} X\, M^{-1} = M^{-1} M M^{-1} = M^{-1} $$
$$ \boxed{\;\theta^{(k)} \sim \mathcal N\big(\hat\theta,\; (X^{\mathsf T} C^{-1} X)^{-1}\big)\;} $$
Now do the Bayesian calculation on the original data. With flat prior, $p(\theta \mid D) \propto \exp\!\big[-\tfrac12(X\theta - D)^{\mathsf T}C^{-1}(X\theta - D)\big]$. Expand the exponent about $\hat\theta$; the linear term vanishes by the normal equations, leaving
$$ p(\theta \mid D) \;\propto\; \exp\!\Big[-\tfrac12\,(\theta-\hat\theta)^{\mathsf T} M\, (\theta-\hat\theta)\Big] \;=\; \mathcal N\big(\hat\theta,\; M^{-1}\big) $$
Same mean, same covariance: the replica ensemble samples the Bayesian posterior exactly. And since the $\chi^2$ is exactly quadratic here with Hessian $H = 2M$, the Hessian covariance $2H^{-1} = M^{-1}$ agrees too. In the linear-Gaussian world there is one uncertainty and three ways to compute it — the triple equivalence §5.3 elevates to a theorem. Every disagreement between methods is therefore a measurement of departure from this world.
Let $T(\theta)$ be nonlinear and expand about the data fit: $T(\hat\theta + \delta) = T(\hat\theta) + J\delta + \tfrac12\,\delta^{\mathsf T} K\, \delta + \dots$ with $J$ the Jacobian and $K$ the model curvature. Solving the replica normal equations order by order in $\eta$ gives $\delta^{(k)} = A\eta^{(k)} + Q\big(\eta^{(k)}, \eta^{(k)}\big) + \dots$ with $A = M^{-1}J^{\mathsf T}C^{-1}$ linear and $Q$ a quadratic form built from $K$. Two consequences: the ensemble mean shifts, $\mathbb E[\delta] = \mathbb E[Q(\eta,\eta)] \neq 0$ — a bias of order (curvature × data variance) that no amount of replicas removes — and the ensemble distribution skews, in a way that need not match the skew of the Bayesian posterior (which weights parameter space by the likelihood, not by the pushforward of data noise through the fitting map). Costantini et al. (arXiv:2404.10056) study exactly this gap quantitatively and show the replica distribution and the posterior part company at second order — the two agree only in the linear/Gaussian limit proved above. §5.3's worked example makes the difference concrete in one dimension.
The machinery in practice: stopping, selection, and the loss arithmetic
Real replica fits (NNPDF's, and ours) add two pieces of procedure. First, cross-validation early stopping: each replica's data are split (ours: 80/20), the optimizer minimizes on the training fraction and halts when validation loss stops improving (patience 20, min-delta $10^{-5}$ in PPDF-5). This regularizes: it stops the fit from chasing the noise component that the finite parametrization could absorb, effectively reducing the fitted degrees of freedom below $p$. What it obscures is equally real: the stopped endpoint is no longer the minimum of any stated objective, so the ensemble's statistical meaning now depends on optimizer, initialization, split, and patience — pure procedure, nowhere in the likelihood. Second, post-fit selection: replicas with anomalous final loss are cut. The cut needs arithmetic, not folklore, and the arithmetic is worth deriving.
Fit a replica $D + \eta$ and evaluate its training loss per point at the endpoint. For a linear model the residual is $r = r_D + (\mathbb 1 - P)\eta$, where $r_D = D - X\hat\theta$ is the data's own best-fit residual and $P = XM^{-1}X^{\mathsf T}C^{-1}$ the projector onto the model span ($P$ is idempotent and self-adjoint in the $C^{-1}$ metric). Then, since the cross term averages to zero and $\mathbb E[\eta^{\mathsf T}(\mathbb 1-P)^{\mathsf T}C^{-1}(\mathbb 1-P)\eta] = \operatorname{tr}(\mathbb 1 - P) = N - p$:
Now count layers in each of our settings. L0 closure: $D = T(\theta_{\rm truth})$ exactly, so $\chi^2_{\min}(D) = 0$ and loss/pt $\approx 1 - p_{\rm eff}/N \approx 1$ — one noise layer, the replica's own. Measured (PPDF-5): 1.009. L1 closure: $D = T(\theta_{\rm truth}) + \varepsilon$ carries real-style noise, so $\chi^2_{\min}(D)/N \approx 1 - p_{\rm eff}/N$, and the replica adds an independent layer $\eta$: loss/pt $\approx 2 - 2p_{\rm eff}/N \approx 2$. Measured: 2.000. (Early stopping keeps $p_{\rm eff}$ well below $p$ and leaves a little unminimized residual, pushing the number back up toward the layer count — the two effects nearly cancel here.) Real data (PPDF-13, MSHT20 model): $\chi^2_{\min}/N = 1.002$ (PPDF-12), so expect $\approx 1 + 1.002 \approx 2.0$; measured: 2.05. This arithmetic is why PPDF-5's L1 run needed the post-fit cut moved to 3.0 (the framework default 1.5 assumes one layer and would have executed the entire healthy ensemble), and why PPDF-13's cut was pre-declared at 3.5: comfortably above the expected 2.0 plus its replica-to-replica scatter, low enough to catch genuine optimizer failures — 83/100 passed, and one further replica exported non-finite parameters and was dropped, leaving 82. Declaring the threshold before looking is what separates a selection rule from tuning.
Three honest virtues: no second derivatives (no $10^7$-entry Hessian of a neural network, no indefiniteness crisis); each replica is an independent job, so the method scales horizontally with perfect efficiency; and non-Gaussian data features propagate automatically, where the Hessian is congenitally symmetric. The price: the ensemble is a distribution over outcomes of a procedure — initialization, optimizer, split, stopping, selection — and only in the linear-Gaussian limit does the procedure's fingerprint provably vanish from the answer. Where data constrain tightly, the likelihood grips every replica and the fingerprint is invisible. Where data go silent, the fingerprint is the band. Our three-way comparison (PPDF-15) displays both regimes on the same plot.
The flat-direction failure — quieter than the Hessian's
Recall §5.1: along a flat direction ($Jv \approx 0$) the Hessian method diverges or goes undefined — a loud failure. Replicas fail silently. The likelihood does not depend on $v^{\mathsf T}\theta$, so no replica's noise draw exerts any pull along $v$: the optimizer's endpoint component along $v$ is simply wherever initialization put it, minus whatever drift the optimizer's dynamics and early stopping happened to produce. The ensemble spread along $v$ therefore measures the initialization distribution filtered through the procedure — a perfectly reproducible number with no inferential content. It can come out small (bands that look confident about the unconstrained: PPDF-6 watched MC valence bands hug the truth purely because initializations centred near it) or large (bands inflated by wide inits), with no rule connecting either to a statement about the proton. On real data this is exactly what PPDF-15 measures: replica-to-Bayesian width ratios of 8.1 ($\Sigma$), 10.3 ($T_8$), 0.6 ($V$) — over- and under-coverage, in the same fit, with central offsets reaching $7\sigma$. When are replicas excellent? When their theorem's premises approximately hold: data-rich problems where $d_{\rm eff}$ approaches $p$ (no silent directions), mild nonlinearity over the scale of the noise, and a compute budget that wants horizontal scaling. Global fits with $\sim$4600 points are far closer to that regime than our deliberately information-starved 648-point laboratory — which is precisely why the laboratory is diagnostic.
Worked example — the equivalence theorem by hand
Given: straight-line model $T(x) = \theta_1 + \theta_2 x$ at $x = 0, 1, 2$; data $D = (1.00,\, 1.50,\, 2.00)$; $C = \sigma^2\mathbb 1$ with $\sigma = 0.1$. So $X^{\mathsf T} = \big(\begin{smallmatrix}1&1&1\\0&1&2\end{smallmatrix}\big)$ and the data fit is exact: $\hat\theta = (1.0,\, 0.5)$.
Step 1 — analytic posterior. $M = X^{\mathsf T}X/\sigma^2 = 100\big(\begin{smallmatrix}3&3\\3&5\end{smallmatrix}\big)$; $M^{-1} = \frac{1}{600}\big(\begin{smallmatrix}5&-3\\-3&3\end{smallmatrix}\big)$. So $\sigma(\theta_1) = \sqrt{5/600} = 0.091$, $\sigma(\theta_2) = \sqrt{3/600} = 0.071$, correlation $-3/\sqrt{15} = -0.77$.
Step 2 — the replica map, explicitly. $\delta\theta^{(k)} = M^{-1}X^{\mathsf T}C^{-1}\eta^{(k)}$ works out to $\delta\theta_1 = (5\eta_1 + 2\eta_2 - \eta_3)/6$, $\;\delta\theta_2 = (\eta_3 - \eta_1)/2$. Sanity check: $\operatorname{Var}(\delta\theta_1) = \sigma^2(25+4+1)/36 = 5/600$ ✓.
Step 3 — five draws $\eta^{(k)} \sim \mathcal N(0, 0.1^2\mathbb 1)$: $(+0.12, -0.05, +0.08)$, $(-0.21, +0.14, +0.03)$, $(+0.07, +0.19, -0.11)$, $(-0.04, -0.13, -0.02)$, $(+0.16, -0.08, +0.22)$. The map gives $\delta\theta_1^{(k)} = \{+0.070, -0.133, +0.140, -0.073, +0.070\}$ and $\delta\theta_2^{(k)} = \{-0.020, +0.120, -0.090, +0.010, +0.030\}$.
Step 4 — ensemble vs posterior. Sample standard deviations: $s(\theta_1) = 0.113$ vs analytic 0.091; $s(\theta_2) = 0.077$ vs 0.071. A 5-replica standard deviation carries relative scatter $\approx 1/\sqrt{2(K-1)} = 35\%$, so both agree well within sampling noise — and the agreement is not approximate in structure: Step 2 shows each $\theta^{(k)}$ is an exact draw from $\mathcal N(\hat\theta, M^{-1})$, the posterior. At $K = 100$ (our production size) the band converges to the analytic answer to $\sim 7\%$.
Apply $\;\mathbb E[\text{loss/pt}] \approx \chi^2_{\min}/N + 1\;$ (corrections: $-p_{\rm eff}/N$ per absorbed layer, $\lesssim$ a few %, partly cancelled by incomplete minimization under early stopping):
L0 (PPDF-5, LH model): $\chi^2_{\min}/N = 0$; layers = replica noise only → expect $\approx 1$. Measured 1.009.
L1 (PPDF-5): $\chi^2_{\min}/N \approx 1$ (data layer $\varepsilon$, unabsorbable); + replica layer $\eta$ → expect $\approx 2$. Measured 2.000.
Real data (PPDF-13, MSHT20): $\chi^2_{\min}/N = 1.002$ from PPDF-12 — nature's "noise layer" relative to the model; + replica layer → expect $\approx 2.00$. Measured 2.05.
The ledger closing to within a few % across three regimes is a nontrivial audit: it says the pseudo-data generation, the optimizer, the stopping rule and the covariance are all statistically consistent — the same role the closure-test identities play for the Bayesian pipeline (§4.4). Note what it does not audit: the flat directions, which contribute nothing to any loss and are exactly where the method's answer is procedure, not inference.