One card per registered result; IDs are immutable. Each card has four tabs:
Results (what we produced and what you're looking at), Explanation (the physics and the method),
Code (exactly how to reproduce it), Notes (the story of how it happened). Status:
verified ·
provisional (caveats attached) ·
lesson (a failure that taught us something).
1 · Calibration — fitting planted protons whose true shape we know.
Synthetic data, real machinery: these results test the instrument, not the proton.
PPDF-1verified · replication
Clean drill passed: fed noise-free synthetic data, the fit recovers the planted proton near-exactly.
What you're looking at: three numbers, no plot — this result is the numbers.
The first says the best fit misses the 648 synthetic data points by 0.0002 of an error bar each: since the data
carried zero noise, anything much larger would have signalled a bug. The paper we're replicating quotes
χ² ~ 10⁻⁴–10⁻⁵ for this same test (Table 3.1) — we land in the same regime. The third number is a health
identity (explained in the next tab) that only holds when the sampling is behaving; here it holds to 3%.
The physics. A proton's structure is encoded in PDFs — for each parton type, x·f(x) says
how the proton's momentum is shared (wiki). Our synthetic "measurements" are
deep-inelastic structure functions F₂(x, Q²): the rate of electrons scattering off the proton at momentum
fraction x and resolution Q. Theory connects the two linearly: prediction = FK ⊗ f, where the FK table
precomputes the QCD matrix element and DGLAP evolution (wiki).
The model being fitted: the Les Houches benchmark shape — per flavour combination,
x·f(x) = A x^α (1−x)^β (1 + ε√x + γx) — 13 free parameters after sum rules. The truth was planted using this
same shape, so exact recovery is possible.
The statistics. Misfit is χ²(θ) = [D − T(θ)]ᵀ C⁻¹ [D − T(θ)] with C the experimental
covariance (t0 prescription, wiki). Bayes: posterior p(θ|D) ∝ e^(−χ²/2) π(θ), uniform
prior boxes π. Nested sampling (UltraNest, 200 live points) explores the 13-dim posterior and computes the
evidence Z = ∫ e^(−χ²/2) π dθ. On noiseless data the optimum is χ² = 0; the health identity
⟨χ²⟩ = χ²_min + complexity diagnoses honest sampling (wiki).
# environment (one-time)
micromamba create -n colibri-dev -f colibri/environment.yml
micromamba run -n colibri-dev pip install -e colibri/colibri/examples/les_houches_example
# run the fit (auto-downloads FK tables + PDF sets on first run)
micromamba run -n colibri-dev les_houches_exe lh_fit_closure_test.yaml
2 Jul, afternoon. First full flight of the machinery, and it flew on the first attempt — the
env build was the only slow part. The 8m22s wall time included downloading four decades of experimental data.
The moment we trusted the statistics was seeing posterior-avg χ² match the Bayesian complexity — that identity
only holds when the posterior is behaving. Every later result stands on this card.
Provenance: theory 40000000 · truth PDF LH_PARAM_20250519 ·
artifacts output/lh_fit_closure_test/ · 2 Jul 2026, Apple-Silicon CPU.
PPDF-2verified · our own check
The full uncertainty band is honest: the planted truth stays inside it at 100% of checkpoints, in every fitted flavour.
What you're looking at: four panels, one per fitted flavour combination (Σ = all quarks
summed; g = gluon; V, V3 = valence combinations), each showing x·f versus momentum fraction x (log scale,
10⁻³ → 1). The teal band is our fit's 68% uncertainty; the dashed black curve is the planted
truth. The test: the dashed line must stay inside the band — it does, everywhere. The grey region
(x < 0.01) marks where our datasets carry no measurements. Note V/V3: enormous bands — the fit honestly
admitting the data barely constrains valence separation, not a failure.
100%
of 120 x-points with truth inside the 1σ band, all four flavours
0.42σ
worst excursion anywhere (gluon)
0.6–0.7%
posterior-mean deviation from truth (Σ, g, median)
Why this check exists. PPDF-1 proved the best fit lands on the truth; this checks
the stronger claim that the whole probability distribution is faithful. A fit can have a perfect central
value and lying error bars — and in precision physics the error bar is the product.
The method. Each posterior sample θ becomes curves x·f(x) for all flavours. At each of
120 x-points we compute the pull = [mean(x) − truth(x)] / σ(x), with σ the band half-width
(wiki). Faithful uncertainty means |pull| ≲ 1 nearly everywhere and truth-coverage
near the Gaussian 68%. The truth curve is evaluated outside the fitting code — straight from its LHAPDF
grid via validphys — so the check shares no code with the thing being checked.
The flavour basis. Curves live in the evolution basis
(wiki): Σ and g dominate F₂ and come out tight; V and V3 barely enter
electromagnetic DIS, so the posterior correctly widens toward the prior — the
flat-direction phenomenon that becomes central in PPDF-4/6.
# reload fitted model + posterior samples, evaluate truth via validphys
micromamba run -n colibri-dev python scripts/verify_closure_truth.py
# core of the comparison:
model = dill.load(open(f"{RUN}/pdf_model.pkl", "rb"))
curves = jax.vmap(model.grid_values_func(xgrid))(posterior_samples) # (1715, 14, 120)
truth = convolution.evolution.grid_values(pdf, FK_FLAVOURS, xgrid, [Q0])
pull = (curves.mean(0) - truth) / curves.std(0)
This card exists because "min χ² is tiny" leaves a loophole: the central curve could be right
while the band lies. The V/V3 panels were the day's teaching moment — enormous bands that are correct,
because the data genuinely doesn't know. Honest ignorance vs false confidence became the theme of the whole
three-method study (PPDF-6).
Provenance: script on the PPDF-1 posterior ·
outputs output/lh_fit_closure_test/verification/ · 2 Jul 2026.
PPDF-3verified · replication
Noisy drill passed: under realistic experimental noise the fit reaches χ²/N ≈ 1 — it fits the signal, not the noise.
0.977
posterior-avg χ²/N (min 0.971)
3.88
Bayesian complexity (PPDF-1 gave 4.12)
−321.8 ± 0.15
evidence logZ, converged
What you're looking at: the same fit as PPDF-1 but with realistic noise injected into
the synthetic data — and the target inverts: 0.977 ≈ 1 is the pass. Perfection (χ² ≈ 0) would mean the
fit memorized the noise; well above 1 would mean it missed real structure. The paper's benchmark is 1.00–1.01;
a single noise throw scatters this number by ±0.056, so 0.977 is statistically on target.
The setup. Level-1 closure: pseudodata D = T(θ_truth) + η, with η drawn from the
real experimental covariance matrix C — the same correlated systematics, normalizations and statistical
errors the actual SLAC/BCDMS measurements carry. Statistically indistinguishable from real data, except the
underlying truth is known.
Why χ²/N ≈ 1 is the target. If the model family contains the truth, the best fit removes
the signal and what remains is noise: E[χ²_min] ≈ N − n_eff, so χ²/N ≈ 1 for N ≫ n_eff. Much below 1: the model
absorbs noise (overfitting); above: underfitting or mis-modelled errors. Expected fluctuation of χ²/N is
√(2/N) ≈ 0.056 at N = 648 — the yardstick for "consistent with 1".
Same machinery, opposite criterion — that's what makes L0+L1 a complete calibration:
L0 catches bugs (demands perfection), L1 catches dishonest statistics (forbids perfection).
# identical runcard to PPDF-1, ONE line changed:closure_test_level: 1# noise drawn from experimental covariance
micromamba run -n colibri-dev les_houches_exe lh_fit_closure_test_L1.yaml
2 Jul, evening. The satisfying part was the direction of the number: 0.977, slightly under
1.00 and inside the noise scatter. At 0.5 we'd have caught noise-memorization; at 2.0, under-fitting. It landed
exactly where an honest machine should. With both drills passed, we declared the machinery ready for real
data.
The Hessian method finds the same optimum as the Bayesian sampler — then fails exactly where theory predicts: on flat directions.
0.04%
Bayesian-vs-Hessian agreement on the L1 minimum (628.9 vs 629.2)
1.9×10⁻⁸
L0 total χ² — gradient descent polishes to machine-level recovery
~6σ
Hessian gluon band's worst miss of the truth at L1
What you're looking at: two validations and one failure, all quantified. The agreement
numbers say two completely different algorithms (gradient descent vs nested sampling) found the same best fit —
strong evidence neither is buggy. The 6σ says the Hessian's error band is nonetheless unfaithful in
places: it under-covers the known truth on the gluon. Internally, the down-valence parameters ran to ~10⁹ with
χ² unchanged — the smoking gun of a flat direction.
The Hessian method (used by the MSHT and CT groups): find the minimum θ₀, expand
quadratically — χ²(θ) ≈ χ²₀ + ½ (θ−θ₀)ᵀ H (θ−θ₀) — and quote uncertainties along the eigenvectors of the
curvature matrix H at Δχ² = T² (tolerance T = 1 here). The construction assumes the likelihood is a bowl.
Where it must fail. Electromagnetic F₂ measures flavour sums; the valence split
barely affects predictions. Along such flat directions the curvature is
~0: H has near-zero eigenvalues, 1/√λ blows up, and the quadratic picture is meaningless — the minimizer's
parameters drift to absurd values (10⁹) at constant χ², and the quoted band becomes an artifact. The Bayesian
posterior has no such assumption: it simply widens to the prior. This is Sec. 2.4 of the Colibri paper,
reproduced live.
Method note: gradient descent (Adam, 3 restarts × 20k epochs) with the same likelihood,
data and t0 as the Bayesian fit — only the inference philosophy differs, which makes the comparison clean.
# gradient-descent methodology block (settings OURS — the paper publishes none)
optimizer_settings:
optimizer: adam
optimizer_hyperparams: {learning_rate: 2e-3}
max_epochs: 20000
param_initialiser_settings: {type: normal, random_seed: 0} # means/stds: dicts by param
hessian_settings: {tolerance: 1.0, iter_init: 3, rng_seed: 1}
actions_: [run_hessian_fit]
The eigenvector members with entries of ±10⁹ looked like a bug for about ten minutes; then χ²
came back identical along those directions and it clicked: not a bug — the phenomenon. Also root-caused
en route: the initializer wants means/stds as dicts keyed by parameter name, not lists. Watching the failure
mode live beat reading about it.
The Monte-Carlo replica method lands on its theoretical expectations to three decimals: loss/point 1.009 (clean) and 2.000 (noisy).
1.009
mean final training loss/pt over 100 replicas, L0 (theory ≈1)
2.000
same at L1 (theory ≈2 — two noise layers)
100/100
replicas passing post-fit selection at both levels
What you're looking at: a statistical bullseye. Each replica's final loss has a sharp
theoretical expectation from pure noise arithmetic (next tab); measuring 1.009 and 2.000 against predictions of
≈1 and ≈2, over 200 independent fits, says the replica machinery is statistically faithful.
The replica method (NNPDF's approach): create K pseudo-datasets D_k = D + η_k with η_k
drawn from the experimental covariance, fit each independently, and quote the spread of the K fits as the
uncertainty — bootstrap in spirit. Each fit uses gradient descent with cross-validation early stopping (hold
out 20% of points; stop when held-out loss stops improving).
The loss arithmetic. Fitting replica k, the best fit removes the signal; what remains is
the replica noise → loss/pt ≈ 1 at L0 (one layer: the MC fluctuation η_k). At L1 the "data" already carries
closure noise and the replica adds its own: two independent layers → loss/pt ≈ 2. Landing on both means the
noise model, covariance handling and early stopping all behave.
The known limitation (why this lab exists): the replica ensemble provably matches the
Bayesian posterior only for models linear in parameters (arXiv:2404.10056). These parametrizations are
not linear — hence the three-method comparisons of PPDF-6 and the MSHT20 phase.
# 100 replicas, 6-way parallel, then post-fit selection
seq 1 100 | xargs -P 6 -I{} micromamba run -n colibri-dev \
les_houches_exe lh_mc_L0.yaml -rep {} -o lh_mc_L0
micromamba run -n colibri-dev mc_postfit lh_mc_L0 # L0: default cut 1.5
micromamba run -n colibri-dev mc_postfit lh_mc_L1 -c 3.0# L1: two noise layers
The L1 post-fit failure ("0 of 100 replicas pass") briefly looked like a broken campaign. The
resolution was the two-layer arithmetic: expected loss ≈ 2 sits above the default cut of 1.5 by construction.
Raising it to 3.0 is a documented consequence of the noise model, not tuning; the 5σ outlier rejection stayed
on.
Three inference engines on identical planted data: all agree where the data constrains — and split exactly where theory predicts.
What you're looking at: the noisy-drill fits by all three methods, overlaid. Each panel:
x·f versus x for one flavour; teal = Bayesian, orange = MC replicas, purple = Hessian
(bands = 1σ), dashed = planted truth. Left and centre (Σ, gluon — well measured by F₂): the three bands
nearly coincide and all cover the dashed line; small-x offsets are the shared noise throw, not disagreement.
Right panel (V, the valence flat direction): the tell — the teal Bayesian band balloons (honest: the data
doesn't know), while orange and purple stay thin and hug the truth largely by luck of where their fits
started. Same data, same model — only the statistical philosophy differs.
Why they agree in the middle. Where data constrains, the likelihood is sharply peaked
and nearly Gaussian; all three methods reduce to the same answer — the bowl assumption holds, the linear-model
condition approximately holds, the posterior is Gaussian. Agreement there is the expected null result and
validates all three implementations at once.
Why they split at the edges. On a flat direction there is no bowl (Hessian undefined),
no linearization (replica equivalence broken, arXiv:2404.10056), and no data pull (posterior → prior). Each
method's character shows: Bayesian = honest ignorance; Hessian/MC = accidental confidence. On planted data we
can adjudicate — the truth is known — and the wide band is the correct answer.
Why it matters: real analyses quote these bands as "the proton's uncertainty". This
picture is the cleanest demonstration of when the standard shortcuts are safe and when they quietly aren't.
# overlay all three posteriors/ensembles against the truth
micromamba run -n colibri-dev python scripts/compare_methods.py
# Hessian band from eigenvector +/- pairs:
plus, minus = members[0::2], members[1::2]
sigma = sqrt(sum(((plus - minus)/2)**2, axis=0))
The synthesis card — the picture that explains the project to a non-expert in ten seconds.
Marked verified because every band traces to an independently verified card (PPDF-1, 3, 4, 5). Its real-data
sibling — same three engines, no truth line to lean on — is the MSHT20 phase's final deliverable.
2 · Toy models vs real measurements — the same machinery pointed at actual
SLAC + BCDMS data (648 points, 1970s–80s experiments). No truth line exists here; models are judged by fit quality
and by the Bayesian evidence.
PPDF-7verified · first real data
The 13-parameter toy meets real measurements and lands at χ²/N = 5.60 — a measurement of the model's rigidity, not a machinery failure.
5.60
min χ²/N on real data (3626.5 / 648)
−1835.1 ± 0.16
evidence logZ — the baseline all later models are judged against
~35 min
wall time, laptop CPU
What you're looking at: the first fit where nobody knows the right answer. 5.60 means
the model misses the average data point by 2.4 error bars — far from the ≈1 of a good description. The
calibration drills (PPDF-3 scored 0.977 on this exact machinery) exonerate the pipeline, so this number
measures the model: a 2005-era benchmark shape is too stiff for 648 real measurements. The logZ became
the incumbent for every challenger.
The data. Real F₂ structure-function measurements: SLAC (late 1970s) and BCDMS (1980s,
CERN muon beam) on hydrogen and deuterium targets — 648 points after kinematic cuts, spanning roughly
0.01 < x < 0.75. The proton/deuteron pairing gives (weak) sensitivity to the u–d flavour split. Full
correlated systematics via the NNPDF covariance implementation; t0 prescription for normalizations.
What a big χ² means here. Three suspects for χ²/N ≫ 1: (a) the model can't bend into
the data's shape (rigidity), (b) internal data tensions, (c) theory settings (perturbative order, deuteron
corrections, cuts). Calibration eliminated the machinery itself. Distinguishing (a) from (b)/(c) required a far
better model — the arc that ends at PPDF-12, where (b) and (c) are acquitted too.
The evidence baseline. logZ = ∫ e^(−χ²/2) π(θ) dθ integrates fit quality over everything
the prior allows (wiki). −1835.1 is meaningless alone; it exists to be
compared against, at identical data/theory/t0, by every later model.
# identical to PPDF-1's runcard, ONE line changed:closure_test_level: false# REAL measurements, no pseudodata
micromamba run -n colibri-dev les_houches_exe lh_fit_realdata_dis.yaml
2 Jul, night. Sampling efficiency halved relative to closure runs — the first sign real data
fights back. The correct reading of 5.6 took a beat: the drills prove the pipeline, so the misfit is a
measurement of model rigidity. It took until PPDF-12 to prove that reading right.
Grid model v1: a stalled sampler, then a model structurally blind to its own data — and an evidence referee that caught it immediately.
8.02
min χ²/N — worse than the 13-parameter toy despite 24 parameters
ΔlogZ ≈ −808
evidence penalty vs the toy baseline — missing physics, priced
290×
speed-up after configuring the slice sampler
What you're looking at: a disqualified result kept on the books. A more flexible model
scoring worse than a rigid one is the anomaly; the diagnosis (next tab) found the model couldn't
represent the physics the data measures. The −808 shows the evidence catching that within a single fit.
Lesson 1 — sampling. Nested sampling above ~10 dimensions needs a step (slice) sampler;
without the config block, UltraNest falls back to region rejection sampling, whose acceptance collapses
exponentially with dimension: 137M likelihood calls overnight at 0.003% efficiency. One config block restored
healthy sampling — 290× faster.
Lesson 2 — flavour masking. The grid model's flavour_mapping zeroes
unlisted flavours in the FK convolution. Among the masked: T3 = u⁺ − d⁺ — the combination that makes a proton
differ from a neutron, hence proton F₂ differ from deuteron F₂. Our dataset is precisely p and d measurements:
the model was structurally unable to describe the difference between its own data subsets, at any parameter
values. Published grid examples are closure tests where truth and model share the mask (it cancels); real data
was off-label.
Why the evidence caught it: logZ integrates fit quality — a model that cannot
reach the data's structure has uniformly poor likelihood, and no parameter count compensates.
The evidence prices missing physics, not complexity.
# the flaw (v1): masked flavours -> FK columns zeroed -> p/d degenerateflavour_mapping: ["\Sigma", "g", "V"]# no V3, T3, T8...# the stall + the fix (REQUIRED at ndim >~ 10):
ultranest_settings:
SliceSampler_settings: {nsteps: 100}# 290x speedup
3 Jul, morning. The stall was root-caused in Colibri's ultranest_fit.py (no step
sampler attached without the block), not guessed at. The masking flaw surfaced only after convergence, when
χ²/N = 8.02 made no sense for a more flexible model; tracing the FK masking revealed the p–d blindness.
Pedagogically the best failure the lab has produced.
First model contest on real data: the evidence prefers 13 honest parameters over 36 constrained ones by ΔlogZ ≈ −109 — verdict held loosely, caveats attached.
5.84
min χ²/N, 36-param grid with full flavours (toy: 5.60)
−1943.9 ± 0.28
logZ → ΔlogZ ≈ −109 vs the toy baseline
16.95
Bayesian complexity — the data uses ~17 of 36 parameters
What you're looking at: the flavour fix worked exactly as diagnosed (8.02 → 5.84), yet
the flexible grid still can't match the rigid toy, and the evidence prefers the toy decisively. Both models
flooring near χ²/N ≈ 5.6–5.8 — where global fits describe these same datasets at ≈1 — was the open question
this card handed to the MSHT20 phase (answered in PPDF-12).
Bayes-factor logic. Two models on identical data: ΔlogZ is the odds update the data
provides — ≈5 is already strong (~150:1), so −109 is overwhelming on its face. But the evidence integrates over
the prior, and verdicts inherit prior choices (wiki).
The caveats that keep this provisional: (1) the grid's canonical prior confines it to a
±5σ band around the published NNPDF4.0 proton while the toy roamed free uniform boxes — an uneven field;
(2) 6 x-nodes per flavour may be too coarse to bend where the data needs; (3) the grid carried no sum rules.
What survives the caveats: the workflow — two full posteriors and evidences on identical
data/theory/t0, compared like-for-like. That machinery is the point, and it's what the MSHT20 phase then used
to settle the floor question.
3 Jul, afternoon. The verdict we hold loosely; what we don't: the evidence-comparison workflow
ran end-to-end on real data (project goal 3), and the shared χ²/N floor became the question that shaped the
MSHT20 phase — where it was answered: the floor was the models.
3 · MSHT20 — porting a big-three, real-world parametrization into the Bayesian
machinery: prove the port, drill it on planted data, then fit real measurements. Ends at Maria's goal: Hessian vs
MC vs Bayesian uncertainties on the same real data.
PPDF-10verified · parity gate
The MSHT20 parametrization, ported to Colibri and proven exact: 96 parity checks against its reference implementation.
numerical momentum sum on 3,000 grid points (analytic: exact)
What you're looking at: a software result with physics stakes. Four levels of checks —
special-function integrals, every shape function, all four sum-rule computations, and the assembled flavours —
each compared point-for-point against FPPDF's own code, at the published parameter values and five
random ±10% perturbations. Everything the fits will ever call agrees to better than a part per billion. The
random draws caught a hidden ×A_g factor invisible at publication values — the reason the gate exists.
The parametrization (MSHT20, Eqs. 2–8 of arXiv:2012.04684): each basis distribution is
x·f(x) = A (1−x)^η x^δ (1 + Σᵢ aᵢ Tᵢ(y)) with Chebyshev polynomials Tᵢ in y = 1 − 2√x — a rigid but
well-crafted shape whose Chebyshev terms add controlled wiggle room. Basis: u_V, d_V, sea S, s₊, s₋, d̄/ū,
two-term gluon, charm (zero here: perturbative-charm convention). Four sum rules are imposed
analytically — valence counting (∫u_V = 2, ∫d_V = 1), momentum (Σ ∫x·f = 1, fixing A_g), zero net
strangeness (fixing the s₋ crossing x₀) — via closed-form Euler-Beta integrals, which keeps the model exact and
fast.
The port method. Transcribe FPPDF's implementation function-by-function into JAX
(differentiable, GPU-ready, what Colibri consumes), then refuse to run any fit until the transcription proves
itself: parity at machine precision across parameter space, plus physics sanities (momentum sum, valence
counting). Convention note: FPPDF's fitting basis uses S = 2(ū+d̄)+s₊ (NNPDF-style) vs the MSHT20 paper's
S = ū+d̄+s+s̄ — we match the code, bit-for-bit.
# run the parity gate (FPPDF installed as reference)
micromamba run -n colibri-dev python analysis/msht20/parity_test.py
# PARITY REPORT — 96 checks, worst rel dev = 1.257e-05 (non-gating Ic8)# RESULT: PASS — port is exact within tolerance# the bug the jittered draws exposed (FPPDF pdfs.py, int_g1_msht):
out = out * agp # multiplies by the A_g slot — 1.0 at publication values
7–8 Jul. Maria's email landed mid-build and confirmed the plan line-for-line — data, code
source, and the sharpened final goal: Hessian vs MC vs Bayesian on MSHT20. Priors, charm convention and closure
criteria were declared in the log before any fit launched (blind discipline). A port tested only at the
publication point would have carried a silent bug into every fit that followed.
MSHT20 passes its own planted-proton drills at 52 dimensions — and reveals that our DIS data sees only ~1/5 of its shape.
0.0074
clean-drill min χ²/N — planted truth recovered
0.979
noisy-drill avg χ²/N (target ≈1 within ±0.056)
~11 / 52
Bayesian complexity — directions the data actually constrains
What you're looking at: the same two-drill calibration as the toy model, at 52
dimensions. Clean drill: near-zero recovery (residual = sampler resolution — nested sampling explores, it
doesn't polish). Noisy drill: 0.979 ≈ 1, honest under realistic noise. Truth-coverage across five flavours:
65–100% against the 68% Gaussian ideal — calibrated bands, neither padded nor thin. The health identity
⟨χ²⟩ = χ²_min + complexity held exactly in both drills.
In-family truth — the key design choice. The planted proton is the MSHT20 model itself
at its published parameter values, not an external PDF set. With out-of-family truth, a nonzero L0 χ² could mean
machinery bugs or model mismatch — conflated. In-family, any residual is machinery by construction: the
drill isolates what it's meant to test.
The complexity finding. Bayesian complexity ≈ ⟨χ²⟩ − χ²_min counts the parameter
directions the data pins down (wiki). 10–12 of 52 quantifies the information
mismatch between a fixed-target DIS subset and a global-fit parametrization: F₂ on p/d constrains the singlet
and gluon well, valence weakly, strangeness barely. Prediction made and kept: honest, wide posteriors in those
directions on real data.
# closure truth generated FROM the model itself (in-family):
closure_test_pdf: colibri_model
closure_test_model_settings:
model: msht20_model
parameters: { uv_del: 0.3436, uv_eta: 3.7487, ... } # 52 publication values
micromamba run -n colibri-dev msht20_exe msht_closure_L0.yaml # then L1
8 Jul. An order of magnitude heavier than the toy's drills (~5M likelihood evaluations for L0
alone), but the PPDF-8 sampler lesson held: no stalls, clean convergence budgets. After the toy's
over-conservative 100% coverage, seeing 65–100% around the 68% ideal was the posterior behaving exactly as
statistics demands. Real-data fit launched the same hour both drills passed.
What you're looking at: the number this arc was built toward. Same 648 real
measurements, same theory, cuts and t0 as every real-data fit above — only the model changed, and the misfit
collapsed from 5.6 error bars per point to statistically perfect. Two toy-era suspects (data tensions, theory
settings) are acquitted in one stroke: the data was describable all along; the earlier models couldn't bend.
A Bayesian posterior for MSHT20 now exists — posterior samples and corner plot in the artifacts.
Why MSHT20 succeeds where the others couldn't. Against the toy (13 params): Chebyshev
terms give controlled mid-x flexibility exactly where BCDMS is most precise. Against the grid (36 params):
flexibility placed by expertise — smooth shapes with physical small-/large-x behaviour x^δ(1−x)^η,
analytic sum rules, correct flavour content — beats flexibility scattered on x-nodes. Parameter count was never
the story; where the freedom lives is.
What χ²/N = 1.002 asserts. Goodness-of-fit is a pure likelihood statement — no prior
involved: the best-fit MSHT20 curve misses the average point by exactly one error bar, matching what global
fits achieve on these datasets. The prior caveat touches only the evidence comparison: our boxes are
centered on MSHT's published values (from their global fit), which flatters logZ. But ΔlogZ ≈ +1453 is
dominated by raw likelihood (Δχ²/2 ≈ 1488 alone) — no defensible prior re-design flips a gap that size.
What it sets up. This posterior is the Bayesian leg of Maria's three-method comparison;
the MC-replica ensemble and Hessian eigenvector bands on the identical setup complete it (next card).
# identical to the closure runcard with one line changed:closure_test_level: false# REAL measurements
micromamba run -n colibri-dev msht20_exe msht_realdata_dis.yaml
8 Jul. The number came out almost embarrassingly clean: 1.002. Four days of χ²/N ≈ 5.6 across
two models had made data tensions feel likely; the answer was simpler — neither model deserved the data. The
phase's namesake deliverable exists: a Bayesian posterior for MSHT20. MC and Hessian legs running.
mean final loss/pt (theory ≈2: one noise layer + χ²/N≈1 residual)
2m36s
per replica on laptop CPU, 6-way parallel
What you're looking at: the replica leg of the three-method comparison. The loss
arithmetic landing on its prediction (≈2) says the ensemble machinery behaved; what the ensemble means
as an uncertainty is exactly what PPDF-15 tests.
Same method as PPDF-5, now on real measurements: 100 jittered copies of
the data (noise from the experimental covariance), one gradient-descent fit each, spread = uncertainty. On real
data the expected final loss/pt is ≈1 (replica noise) + χ²/N of the model against the real data (≈1 from
PPDF-12) ≈ 2 — the pre-declared selection cut (3.5) follows from that arithmetic plus
the usual outlier rejection, fixed before results existed.
8 Jul. Ran while the Bayesian posterior was still being verified — the three legs of the
comparison were produced independently, no leg saw another's results before PPDF-15 brought them together.
The 17 rejected replicas are gradient-descent stragglers; one exported row was non-finite and is dropped with
a note rather than silently.
The Hessian method structurally fails here: at 52 parameters on 648 DIS points, the curvature at the minimum is indefinite — there are no meaningful eigenvector bands to quote.
0.9845
χ²/N at the minimum — found identically four times (the fit itself is fine)
11 / 52
eigen-directions with negative (≈zero) curvature — no bowl to measure
±70,000
largest eigenvector member excursions in parameter space
What you're looking at: a clean minimum with unusable curvature. The best fit is
excellent; the quadratic approximation around it is undefined along most directions — roughly 38 of 52
parameters are simply not constrained by this dataset (matching the Bayesian complexity ≈14 of
PPDF-12). The purple explosion in the PPDF-15 plot is this table made visible.
Why it breaks. Hessian errors assume χ² is a bowl: H must be positive-definite, bands
scale as 1/√λ. Where data leaves directions flat, λ → 0 (our finite-difference resolution measures them as
consistent with zero, some marginally negative), and 1/√λ diverges — bands lose meaning. This is
PPDF-4's finding taken to its extreme: at 13 parameters two directions were flat; at 52
parameters on the same data, most are. MSHT's published global fits avoid this by fitting ~4,600 points and
freezing insufficiently-constrained parameters — our setup exposes the method's data-sufficiency requirement.
Computational note (its own lesson): three autodiff routes to this matrix hit
pathological XLA compile times (>40 min each — fused jax.hessian, HVP columns, jitted-gradient FD).
The working method: central finite differences of likelihood values in eager mode — zero compilation,
~5,600 evaluations, accuracy far beyond what T=1 bands need.
8 Jul. Half a day of XLA fighting (three methods, all compile-bound — each attempt logged)
ended with the least glamorous method winning: plain finite differences, no compiler involved. The
indefiniteness then turned a computational chore into the sharpest finding of the comparison. Worth reporting
upstream to Colibri: their Hessian path will bite anyone past ~20 parameters on CPU.
Hessian vs MC vs Bayesian on the same real data: they agree only where data constrains — and only the Bayesian bands carry a calibration certificate.
What you're looking at: five flavour panels, x·f vs x at Q₀; teal = Bayesian,
orange = MC replicas, purple = Hessian, bands = 1σ. Grey: x < 0.01, outside our data. Panels
are y-clipped to the Bayes/MC scale — the purple Hessian band exceeds every frame (that's the
PPDF-14 breakdown, visible). In the gluon panel at mid-x — the direction this data
genuinely measures — teal and orange agree closely. In Σ and T8 the orange replica band runs ~8–10× wider than
teal with centrals up to ~7σ apart: in unconstrained directions the replica spread reflects initialization
luck, not information.
1.8× · <1σ
MC vs Bayes on the gluon: width ratio and central offset — agreement where constrained
0.6–10×
MC/Bayes width ratios across flavours — no consistent rule off-constraint
What this tests. Maria's stated final goal: same model (MSHT20), same 648 real
measurements, same theory and covariance — only the uncertainty philosophy differs. Real data allows no truth
line, so methods are judged against each other plus one external fact: the Bayesian machinery holds a closure
certificate at exactly this dimensionality (PPDF-11: calibrated coverage), which
neither alternative has here.
The pattern, and why. Where the likelihood is sharply peaked (gluon at mid-x), all
honest methods converge — Gaussian likelihood makes replica-equivalence and posterior-Gaussianity approximately
true. Where directions go flat, each philosophy shows its failure mode: Hessian bands diverge (no curvature),
replica bands land wherever gradient descent starts and stops (initialization-shaped), Bayesian bands widen to
the declared prior (honest, but prior-dependent — stated on the tin). The practical conclusion: at this
dataset-to-parametrization information ratio, only the Bayesian answer is defensible — and the fix for the
others is more data, which is exactly the full-DIS run in flight.
Scope caveat: this is a statement about methods at THIS information ratio — not about
MSHT20's published errors, which come from a global dataset with parameter freezing.
# one function turns any parameter set into PDF curves at Q0 (same model for all three)
f = model.grid_values_func(xgrid) # theta -> (14 flavours, 120 x-points)
curves_of = lambda P: np.array(jax.vmap(f)(P))
# the three uncertainty objects, from each method's own output
bay = curves_of(posterior_samples) # 17,496 Bayesian samples (PPDF-12)
mc = curves_of(replica_params[finite]) # 82 MC replicas (PPDF-13; 1 non-finite dropped)
plus, minus = curves_of(eig_members[0::2]), curves_of(eig_members[1::2])
h_sig = np.sqrt(np.nansum(((plus - minus) / 2)**2, axis=0)) # Hessian: quadrature over eigen-pairs (PPDF-14)
bands = {"Bayesian": (bay.mean(0), bay.std(0)),
"MC replicas": (mc.mean(0), mc.std(0)),
"Hessian": (h_central, h_sig)}
# the two comparison metrics, per flavour, in the data region x > 0.01
width_ratio = np.nanmedian(s_method / s_bayes) # how wide vs the calibrated band
offset = np.nanmax(|m_method - m_bayes| / s_bayes) # central drift in Bayesian sigmas
micromamba run -n colibri-dev python scripts/compare_msht_methods.py
# summary CSV: output/msht_threeway/threeway_summary.csv
8 Jul. The deliverable Maria named in her email, produced ~30 hours after the email arrived.
The expected result was "three bands, mild differences"; the actual result is starker and more useful — a
live demonstration of when each standard method can and cannot be trusted, with the Bayesian calibration
certificate as the differentiator. The full-DIS rerun of this comparison (5× the data) will show whether the
alternatives recover.
