The χ², correlated systematics, and t0
Covariance anatomy, nuisance-parameter equivalence derived, the d'Agostini bias and its fix.
Every number this lab has ever put on a result card — every $\log Z$, every posterior band, every "$\chi^2/N = 1.002$" — flows through one scalar function: the $\chi^2$. This chapter builds it properly: the covariance matrix with its correlated-systematic anatomy, the theorem that makes nuisance parameters and matrix inversion the same thing, the d'Agostini bias that quietly poisons naive fits to multiplicative errors, and the $t_0$ prescription that every runcard in this lab carries. Then: how to actually read a $\chi^2/N$.
The Gaussian likelihood
Experiments publish central values $D_I$ and an uncertainty breakdown; in the Gaussian approximation (excellent for our high-statistics DIS sets) the sampling distribution of the data given a theory $T(\theta)$ is multivariate normal, and the log-likelihood is, up to a constant,
All the subtlety hides in $C$. If it were diagonal, $\chi^2$ would be the familiar sum of squared pulls. It is not diagonal, and the off-diagonal structure is not a detail — it is most of the information about how the data can be wrong.
Anatomy of the covariance
Each uncertainty source in the commondata files (§3.3) contributes according to its correlation type. Statistical and uncorrelated systematics fill the diagonal; each correlated source $a$ — a beam-energy calibration, a detector efficiency, a normalization — shifts many points coherently and contributes a rank-one outer product:
The systematic part has rank $\le N_{\rm sys}$: the experiment is telling you its points can move together in only $N_{\rm sys}$ specific patterns. BCDMS proton alone ships several such sources across 333 points — so $C$ is a dense $648 \times 648$ matrix, built once, Cholesky-factorized once, and reused for all $10^8$ likelihood calls (§3.3's budget).
There is a second, physically transparent way to handle a correlated source: give it a strength $\lambda$, shift the theory by $\lambda \sigma_I$, and pay a unit-Gaussian penalty. For one source on top of diagonal errors $s_I$:
$$ \chi^2(\theta, \lambda) \;=\; \sum_I \frac{(r_I - \lambda\, \sigma_I)^2}{s_I^2} \;+\; \lambda^2 $$
Minimize over $\lambda$: $\;\partial_\lambda \chi^2 = 0$ gives $-2\sum_I \sigma_I (r_I - \lambda \sigma_I)/s_I^2 + 2\lambda = 0$, hence
$$ \lambda^\ast \;=\; \frac{\sum_I \sigma_I r_I / s_I^2}{1 + \sum_I \sigma_I^2 / s_I^2} \;=\; \frac{\sigma^{\mathsf T} S^{-1} r}{1 + \sigma^{\mathsf T} S^{-1} \sigma}, \qquad S \equiv \mathrm{diag}(s_I^2) $$
Substitute back (expand, collect the two $\lambda^{\ast 2}$ terms):
$$ \chi^2(\theta, \lambda^\ast) \;=\; r^{\mathsf T} S^{-1} r \;-\; \frac{\big(\sigma^{\mathsf T} S^{-1} r\big)^2}{1 + \sigma^{\mathsf T} S^{-1} \sigma} $$
Now compare with the Sherman–Morrison identity for the inverse of a rank-one update:
$$ \big(S + \sigma \sigma^{\mathsf T}\big)^{-1} \;=\; S^{-1} \;-\; \frac{S^{-1} \sigma\, \sigma^{\mathsf T} S^{-1}}{1 + \sigma^{\mathsf T} S^{-1} \sigma} \quad\Longrightarrow\quad r^{\mathsf T} (S + \sigma\sigma^{\mathsf T})^{-1} r \;=\; \chi^2(\theta, \lambda^\ast) $$
Identical, term for term. Profiling nuisance strengths with Gaussian penalties is exactly the same computation as inverting the outer-product covariance — one source shown here; $N_{\rm sys}$ sources follow by iterating the update (Sherman–Morrison–Woodbury). So the covariance formulation isn't a lossy shortcut: the optimal systematic shifts are still in there, computable afterwards via $\lambda^\ast$ (that is exactly how experiments' "shifted data" plots are made).
What correlated shifts do to a fit
The equivalence explains the single most counterintuitive feature of correlated $\chi^2$'s: coherent offsets are cheap. Take two points with 1% statistical errors and a shared 5% normalization source, and let the theory sit uniformly 4% below both — four statistical standard deviations! Uncorrelated treatment: $\chi^2 = 4^2 + 4^2 = 32$, catastrophic. With the correlation, use Derivation 1's numbers ($r = (4,4)$, $s = (1,1)$, $\sigma = (5,5)$ in % units): $\sigma^{\mathsf T} S^{-1} r = 40$, $1 + \sigma^{\mathsf T} S^{-1}\sigma = 51$, so
The fit "buys" a 0.78$\sigma$ pull of the normalization source and explains both points for less than one unit of $\chi^2$. Offsets are what the systematic is for; what stays expensive is shape: pull the two residuals apart ($r = (4, -4)$) and the source can't help at all. This is why a fit can sit visibly below a dataset's points and still claim an excellent $\chi^2$ — and why per-dataset $\chi^2$ tables must be read together with the fitted shifts.
Multiplicative errors and the d'Agostini bias
Normalization uncertainties are multiplicative: "$\pm 1\%$ luminosity" means $\pm 1\%$ of the value. Of which value, though — the measurement, or the truth? The naive choice, $\sigma_I = f \cdot D_I$, plants a trap discovered by d'Agostini in fits to correlated data: points that fluctuate down get smaller absolute uncertainties, hence larger weights, and the fit is pulled systematically low.
One true value $T$, two independent measurements with purely multiplicative fractional error $f$, treated naively as uncorrelated with $\sigma_i = f D_i$. The $\chi^2$-optimal combination is the inverse-variance weighted mean:
$$ \hat T \;=\; \frac{\sum_i D_i / (f D_i)^2}{\sum_i 1/(f D_i)^2} \;=\; \frac{1/D_1 + 1/D_2}{1/D_1^2 + 1/D_2^2} \;=\; \frac{D_1 D_2\, (D_1 + D_2)}{D_1^2 + D_2^2} $$
Write a symmetric fluctuation, $D_{1,2} = T(1 \pm \epsilon)$: then $D_1 D_2 = T^2(1-\epsilon^2)$, $D_1 + D_2 = 2T$, $D_1^2 + D_2^2 = 2T^2(1 + \epsilon^2)$, and
$$ \hat T \;=\; T\, \frac{1 - \epsilon^2}{1 + \epsilon^2} \;\approx\; T\,(1 - 2\epsilon^2) \;<\; T $$
always low, at second order in the fluctuation — the noise, which should average away, instead rectifies into a one-sided bias because it leaked into the weights. Note $f$ itself dropped out of the ratio: the bias is driven by the data scatter (in a real fit, statistical scatter feeds it too), and it does not shrink with more data — 300 BCDMS points with a shared normalization just bias more coherently. In global PDF fits the effect reaches several percent on normalizations: fatal at modern precision.
The $t_0$ prescription
The cure (NNPDF, 2009): never let the data set its own multiplicative error bars. Compute the
magnitudes of multiplicative sources from the predictions of a fixed reference, $T^{(t_0)}_I$ — in our runcards,
t0pdfset: NNPDF40_nnlo_as_01180:
The reference is frozen during the fit — the weights can no longer respond to fluctuations, so the rectification mechanism is dead. If the reference is poor, iterate: fit, use the result as the new $t_0$, refit; convergence is fast because the bias is second-order. Two consequences matter for us. First, mildly counterintuitive: the "best-fit" $\chi^2$ is now computed with a covariance built from a different PDF than the one being fitted — that is the point, not a bug. Second, and crucial for this lab's whole program: the covariance defines the likelihood, and evidence only compares likelihoods computed on the same footing. Our entire comparison chain — LH toy (PPDF-7), grid (PPDF-9), MSHT20 (PPDF-12) — used the identical $t_0$ set. Change $C$ between fits and $\Delta \log Z$ silently mixes "better model" with "different question": the $+1453$ would be meaningless.
Reading a $\chi^2/N$, honestly
If the model is right and $C$ is faithful, $\chi^2$ at the true parameters follows a $\chi^2_N$ distribution: mean $N$, variance $2N$, so
Now read our ledger against that yardstick. MSHT20 real data: $1.002$ — dead centre (PPDF-12). L1 closure: $0.977$ and $0.961$ — within one scatter unit, as pseudodata must be (PPDF-3, PPDF-11). The LH toy on real data: $5.60$, i.e. $(5.60 - 1)/0.056 \approx 80$ scatter units — not "a mediocre fit" but a categorical rejection of the model (PPDF-7). Fine print for correlated data: the $\pm\sqrt{2/N}$ applies at the true parameters; after fitting $d$ effective directions the expected minimum drops to $\sim (N - d)/N$ (our Bayesian complexity $\approx 14$ in PPDF-12 predicts $0.978$ — the observed $1.002$ is comfortably within scatter). And always open the per-dataset table: a global 1.0 can hide a 0.7 and a 1.4 in tension, with shared systematics pulled to broker the compromise — the fitted $\lambda^\ast_a$ are part of the fit's story, not bookkeeping.
A good global $\chi^2/N$ says the residuals are consistent with the error budget — nothing more. It does not certify the parametrization outside the data (that's the prior's territory, §4.1), does not exclude compensating wrong physics (theory 40000000's assumptions are frozen into every prediction, §3.3), and by itself cannot arbitrate between models — PPDF-9 and PPDF-7 differ by only $0.24$ in $\chi^2/N$ but by $109$ units of $\log Z$. In this lab the division of labour is explicit: $\chi^2/N$ answers "is this model consistent with this data?", the evidence answers "which model earns the data?", and closure tests (§4.4) answer "are the uncertainties calibrated?". Three different questions; keep three different instruments.
Setup: 3 points, residuals $r = (0.60, -0.20, 0.45)$, statistical errors $s = (0.50, 0.40, 0.60)$, one correlated source $\sigma = (0.30, 0.24, 0.36)$ (a 6% normalization on values $\sim (5,4,6)$).
Route 1 — profile the nuisance. $\sigma^{\mathsf T} S^{-1} r = \frac{0.30 \times 0.60}{0.25} + \frac{0.24 \times (-0.20)}{0.16} + \frac{0.36 \times 0.45}{0.36} = 0.72 - 0.30 + 0.45 = 0.87$; likewise $\sigma^{\mathsf T} S^{-1} \sigma = 0.36 + 0.36 + 0.36 = 1.08$, so $\lambda^\ast = 0.87 / 2.08 = 0.4183$. Shifted residuals $r - \lambda^\ast \sigma = (0.4745,\ -0.3004,\ 0.2994)$ give $\chi^2 = \frac{0.4745^2}{0.25} + \frac{0.3004^2}{0.16} + \frac{0.2994^2}{0.36} + \lambda^{\ast 2} = 0.9007 + 0.5640 + 0.2490 + 0.1750 = \mathbf{1.8886}$.
Route 2 — build $C$ and invert. $C = \mathrm{diag}(0.25, 0.16, 0.36) + \sigma\sigma^{\mathsf T} = \left(\begin{smallmatrix} 0.3400 & 0.0720 & 0.1080 \\ 0.0720 & 0.2176 & 0.0864 \\ 0.1080 & 0.0864 & 0.4896 \end{smallmatrix}\right)$; solving $C z = r$ and taking $r^{\mathsf T} z$ gives $\chi^2 = \mathbf{1.8886}$ — or via Sherman–Morrison directly: $r^{\mathsf T} S^{-1} r - (\sigma^{\mathsf T} S^{-1} r)^2 / (1 + \sigma^{\mathsf T} S^{-1}\sigma) = 2.2525 - 0.87^2/2.08 = 2.2525 - 0.3639 = 1.8886$. Agreement to all digits — and note how the correlation forgave the coherent part of the residuals: naive uncorrelated $\chi^2$ was 2.25.
The lab's standard toy, expanded: one quantity, true value 1.00; two measurements with 10% multiplicative errors, fluctuated to $D_1 = 1.10$, $D_2 = 0.90$.
Naive (data-based) weights: $\sigma_1 = 0.110$, $\sigma_2 = 0.090$, so $w_1 = 82.6$, $w_2 = 123.5$ — the low point gets 49% more weight for no reason but its own downward fluctuation:
$$ \hat T = \frac{1.10 \times 82.6 + 0.90 \times 123.5}{82.6 + 123.5} = \frac{202.0}{206.1} = 0.9802 \;=\; \frac{1 - \epsilon^2}{1 + \epsilon^2}\Big|_{\epsilon = 0.1} $$
a $-2\%$ bias, exactly Derivation 2's $-2\epsilon^2$.
$t_0$ fix: take any reasonable reference, say $T^{(t_0)} = 1.00$: then
$\sigma_1 = \sigma_2 = 0.10$, the weights are equal by construction, and $\hat T = (1.10 + 0.90)/2 =
\mathbf{1.000}$. Robustness check — start from a bad reference, even the biased $T^{(t_0)} = 0.9802$: the
magnitudes are again equal ($0.098$ each), so the very first iteration returns $1.000$ and stays there.
The fixed point is the unbiased answer; the naive scheme's fixed point is not. Scale this toy by 300 — BCDMS's
correlated normalization across hundreds of points — and the $2\%$ becomes a coherent distortion of every PDF
in the fit. That one line in our runcards, use_t0_covmat: True, is what stands between
PPDF-12's $\chi^2/N = 1.002$ and a subtly, systematically low proton.