pith. sign in

arxiv: 2605.20048 · v1 · pith:53656YL5new · submitted 2026-05-19 · ✦ hep-ex

Gaussian Process Eigenmodes for Statistical and Systematic Uncertainties in Template Fits

Pith reviewed 2026-05-20 03:45 UTC · model grok-4.3

classification ✦ hep-ex
keywords Gaussian processestemplate fitsBarlow-Beestoneigenmodessystematic uncertaintiesMonte Carlo statisticsHistFactorylog-Gaussian Cox process
0
0 comments X

The pith

Gaussian process eigenmodes unify statistical and systematic uncertainties in template fits while containing Barlow-Beeston as a limit case.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper proposes replacing the per-bin Barlow-Beeston gamma factors and systematic interpolation parameters in HistFactory templates with a small number of Gaussian-constrained amplitudes derived from eigenmodes. The approach fits log-Gaussian Cox processes to Monte Carlo data to create smooth functional representations and uses their posterior covariance augmented by rank-1 updates for each systematic variation to form a unified eigenmode basis. Truncating to the leading eigenmodes reduces the number of parameters independently of the number of bins. A sympathetic reader would care because this addresses the scaling problems in multi-dimensional analyses and with limited Monte Carlo samples. The paper proves that the construction includes the standard Barlow-Beeston method as a limiting case and that the Gaussian process posterior variance is bounded above by the Barlow-Beeston variance in every bin.

Core claim

The central discovery is that the posterior covariance of a log-Gaussian Cox process fitted to Monte Carlo histogram data, when augmented with rank-1 updates for each systematic shape variation, yields an eigenmode basis that encodes both statistical and systematic template uncertainties. Truncating this basis to the dominant modes replaces the full set of per-bin parameters with a reduced set of Gaussian amplitudes. This construction contains the Barlow-Beeston approach as a limiting case, and the resulting Gaussian process posterior variance is bounded above by the Barlow-Beeston variance at every bin.

What carries the argument

The unified eigenmode basis obtained from the decomposition of the augmented posterior covariance matrix of the log-Gaussian Cox process, which encodes both Monte Carlo statistical uncertainties via the posterior and systematic uncertainties via rank-1 updates.

If this is right

  • The number of uncertainty parameters no longer scales with the number of bins in the template.
  • Statistical and systematic uncertainties are handled within a single coherent framework rather than separately.
  • The method enables the use of smooth functional forms instead of discrete histograms for templates.
  • Truncation to leading modes preserves the essential uncertainty information with reduced computational cost.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • This could allow more efficient fitting procedures in high-dimensional analyses at particle colliders.
  • The variance bound implies that uncertainty estimates may be more precise than in the standard approach in some cases.
  • Similar eigenmode techniques might be applied to other uncertainty sources in experimental physics.

Load-bearing premise

The log-Gaussian Cox process posteriors fitted to the Monte Carlo data accurately represent the template uncertainties so that the eigenmode truncation preserves the essential statistical and systematic information without loss of fidelity.

What would settle it

A numerical test comparing the full per-bin Barlow-Beeston likelihood to the truncated eigenmode version on the same Monte Carlo samples, checking if the fit results and uncertainty bounds match within the expected approximation error from truncation.

Figures

Figures reproduced from arXiv: 2605.20048 by Vincent Alexander Croft.

Figure 1
Figure 1. Figure 1: FIG. 1. Numerical confirmation of the variance bound Eq. ( [PITH_FULL_IMAGE:figures/full_fig_p003_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: FIG. 2. Experiment B truth distribution: three back [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: FIG. 3. The binning dilemma for the Experiment A rare-resonance search. The same [PITH_FULL_IMAGE:figures/full_fig_p006_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: FIG. 4 [PITH_FULL_IMAGE:figures/full_fig_p007_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: FIG. 5. Pull distribution at [PITH_FULL_IMAGE:figures/full_fig_p008_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: FIG. 6. Fit diagnostics for a single Experiment B pseudo-experiment at [PITH_FULL_IMAGE:figures/full_fig_p009_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: FIG. 7. Eigenvalue spectrum for Experiment B. Circles: stat [PITH_FULL_IMAGE:figures/full_fig_p011_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: FIG. 8. Leading four eigenmodes of the combined covariance for Experiment B, visualised as [PITH_FULL_IMAGE:figures/full_fig_p012_8.png] view at source ↗
Figure 11
Figure 11. Figure 11: FIG. 11. Asymptotic CL [PITH_FULL_IMAGE:figures/full_fig_p013_11.png] view at source ↗
Figure 10
Figure 10. Figure 10: FIG. 10. Modifier chain recovery: (a) ˆµ [PITH_FULL_IMAGE:figures/full_fig_p013_10.png] view at source ↗
Figure 12
Figure 12. Figure 12: FIG. 12. The efficient information fraction [PITH_FULL_IMAGE:figures/full_fig_p015_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: FIG. 13. Background estimates for a single Dataset A pseudo [PITH_FULL_IMAGE:figures/full_fig_p015_13.png] view at source ↗
read the original abstract

Template histograms are the foundation of statistical inference at the Large Hadron Collider. The HistFactory likelihood encodes template uncertainty through per-bin Barlow-Beeston gamma factors for Monte Carlo statistical error and through interpolation-based modifiers for systematic shape variations. These two mechanisms scale with the number of bins, which becomes problematic for multi-dimensional analyses and for templates constructed from limited Monte Carlo samples. We propose the use of eigenmode decomposition for efficiently estimating statistical and systematic uncertainties when replacing histogram templates with smooth functional representations derived from log-Gaussian Cox process posteriors fitted to the Monte Carlo data. The posterior covariance, augmented by rank-1 updates for each systematic shape variation, provides a unified eigenmode basis that encodes both statistical and systematic template uncertainty. Truncating to the leading eigenmodes replaces the full set of per-bin gamma factors and interpolation parameters with a small number of Gaussian-constrained amplitudes. We prove that this construction contains Barlow-Beeston as a limiting case and that the Gaussian Process posterior variance is bounded above by the Barlow-Beeston variance at every bin.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript proposes replacing per-bin Barlow-Beeston gamma factors and systematic interpolation parameters in HistFactory-style template fits with a reduced set of Gaussian-constrained amplitudes obtained from eigenmode decomposition of the posterior covariance of a log-Gaussian Cox process fitted to Monte Carlo data, augmented by rank-1 updates for each systematic shape variation. The central claims are a mathematical proof that the construction recovers the Barlow-Beeston method as a limiting case and that the Gaussian Process posterior variance is bounded above by the Barlow-Beeston variance at every bin.

Significance. If the derivations hold, the approach could meaningfully improve the scalability of template fits for multi-dimensional analyses or those with limited Monte Carlo statistics by reducing the number of nuisance parameters while preserving a unified treatment of statistical and systematic template uncertainties. The explicit proof of the limiting case and the variance bound constitutes a clear technical strength.

major comments (2)
  1. [§4.2] §4.2, the limiting-case argument: the demonstration that Barlow-Beeston is recovered requires an explicit statement of the order of limits (GP correlation length to zero versus number of retained eigenmodes to infinity); the current derivation leaves ambiguous whether truncation precedes or follows the limit.
  2. [§5.3, Eq. (32)] §5.3, Eq. (32): the variance-bound proof is stated for the full posterior covariance; it is not shown that the inequality survives projection onto the leading eigenmodes, which is load-bearing for the claim that the truncated model remains conservative.
minor comments (2)
  1. [§2.3] The notation for the rank-1 systematic updates in §2.3 would be clearer if an explicit matrix expression were provided rather than a verbal description.
  2. [Figure 4] Figure 4 caption should report the cumulative variance fraction captured by the retained eigenmodes to allow readers to assess truncation fidelity.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive comments on our manuscript. We address each major comment below.

read point-by-point responses
  1. Referee: [§4.2] §4.2, the limiting-case argument: the demonstration that Barlow-Beeston is recovered requires an explicit statement of the order of limits (GP correlation length to zero versus number of retained eigenmodes to infinity); the current derivation leaves ambiguous whether truncation precedes or follows the limit.

    Authors: We agree that an explicit statement of the order of limits is needed for clarity. In the revised §4.2 we will specify that the Barlow-Beeston limit is recovered by first sending the GP correlation length to zero (recovering independent per-bin Poisson statistics) and only then taking the number of retained eigenmodes to infinity. This ordering prevents any premature truncation from affecting the limiting case. revision: yes

  2. Referee: [§5.3, Eq. (32)] §5.3, Eq. (32): the variance-bound proof is stated for the full posterior covariance; it is not shown that the inequality survives projection onto the leading eigenmodes, which is load-bearing for the claim that the truncated model remains conservative.

    Authors: The bound var_GP(j) ≤ var_BB(j) holds for the full posterior covariance at each bin j. The truncated covariance is the partial sum of the leading eigenmode contributions and is therefore a positive semi-definite reduction of the full covariance matrix. It follows immediately that var_trunc(j) ≤ var_full(j) ≤ var_BB(j) for every bin, so the inequality survives truncation. We will add a short remark in §5.3 making this observation explicit. revision: yes

Circularity Check

0 steps flagged

No circularity: mathematical proof of limiting case and variance bound is independent of fitted inputs

full rationale

The paper's central result is a mathematical proof that the proposed GP eigenmode construction recovers Barlow-Beeston as a limiting case and that GP posterior variance is bounded above by Barlow-Beeston variance at every bin. This follows directly from the properties of the posterior covariance of the log-Gaussian Cox process (with rank-1 systematic updates) and its eigen-decomposition. No step reduces a prediction to a fitted parameter by construction, no self-citation chain bears the load of the uniqueness or derivation, and the proof does not rely on renaming or smuggling an ansatz. The derivation is self-contained against the stated assumptions about the LGCP model; questions of model adequacy for MC data are separate from circularity in the claimed mathematical relationship.

Axiom & Free-Parameter Ledger

2 free parameters · 1 axioms · 0 invented entities

The central claim rests on modeling Monte Carlo templates via log-Gaussian Cox process posteriors whose covariance yields a useful eigenbasis; this introduces GP fitting as the core modeling step beyond standard histogram methods.

free parameters (2)
  • Gaussian Process hyperparameters
    Kernel and other hyperparameters of the log-Gaussian Cox process fitted to MC data; these are chosen or optimized to define the posterior covariance.
  • Number of retained eigenmodes
    Truncation threshold for keeping leading modes; this controls the approximation level.
axioms (1)
  • standard math Standard properties of Gaussian Processes, covariance operators, and eigen decompositions hold for the posterior.
    Invoked when forming the unified eigenmode basis from the posterior covariance plus rank-1 updates.

pith-pipeline@v0.9.0 · 5702 in / 1209 out tokens · 48235 ms · 2026-05-20T03:45:34.316788+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

46 extracted references · 46 canonical work pages · 8 internal anchors

  1. [1]

    Cranmer, G

    K. Cranmer, G. Lewis, L. Moneta, A. Shibata, and W. Verkerke,HistFactory: A tool for creating statisti- cal models for use with RooFit and RooStats, Tech. Rep. CERN-OPEN-2012-016 (2012)

  2. [2]

    Asymptotic formulae for likelihood-based tests of new physics

    G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Asymp- totic formulae for likelihood-based tests of new physics, Eur. Phys. J. C71, 1554 (2011), arXiv:1007.1727

  3. [3]

    ATLAS Collaboration, Observation of a new particle in the search for the Standard Model Higgs boson with the 10 ATLAS detector at the LHC, Phys. Lett. B716, 1 (2012), arXiv:1207.7214

  4. [4]

    CMS Collaboration, Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys. Lett. B716, 30 (2012), arXiv:1207.7235

  5. [5]

    ATLAS Collaboration, A detailed map of Higgs boson interactions by the ATLAS experiment ten years after the discovery, Nature607, 52 (2022)

  6. [6]

    V. Croft and Histimator Collaboration,Histimator: Gaussian process templates for HistFactory likeli- hoods (2026), adds parallel pseudo-experiment driver run pseudo experiments, projected-pull diagnostic, ShapeFactor and InterpCode-4 limit tests

  7. [7]

    Barlow and C

    R. Barlow and C. Beeston, Fitting using finite Monte Carlo samples, Comput. Phys. Commun.77, 219 (1993)

  8. [8]

    J. S. Conway, Incorporating Nuisance Parameters in Likelihoods for Multisource Spectra, inPHYSTAT 2011 (2011) arXiv:1103.0354 [physics.data-an]

  9. [9]

    Alexe, J

    B. Alexe, J. Bendavid, L. Bianchini, and A. Brus- chini, Under-coverage in high-statistics counting experi- ments with finite MC samples, (2024), arXiv:2401.10542 [physics.data-an]

  10. [10]

    Møller, A

    J. Møller, A. R. Syversveen, and R. P. Waagepetersen, Log Gaussian Cox Processes, Scand. J. Statist.25, 451 (1998)

  11. [11]

    P. J. Diggle, P. Moraga, B. Rowlingson, and B. M. Tay- lor,Spatial and Spatio-Temporal Log-Gaussian Cox Pro- cesses: Extending the Geostatistical Paradigm, Vol. 28 (2013) pp. 542–563

  12. [12]

    Modeling Smooth Backgrounds and Generic Localized Signals with Gaussian Processes

    M. Frate, K. Cranmer, S. Kalia, A. Vandenberg-Rodes, and D. Whiteson, Modeling Smooth Backgrounds and Generic Localized Signals with Gaussian Processes, (2017), arXiv:1709.05681 [physics.data-an]

  13. [13]

    Gandrakota, A

    A. Gandrakota, A. Lath, A. V. Morozov, and S. Murthy, Model selection and signal extraction using Gaussian Process regression, JHEP02, 230

  14. [14]

    I. Frid, L. Barak, S. Jairam, M. Kagan, and R. Hyneman, Log Gaussian Cox Process Background Modeling in High Energy Physics, (2025), arXiv:2508.11740 [hep-ex]

  15. [15]

    Whitehorn, J

    N. Whitehorn, J. van Santen, and S. Boeser, Penalized splines for smooth representation of high-dimensional Monte Carlo datasets, Comput. Phys. Commun.184, 2214 (2013)

  16. [16]

    Lindgren, H

    F. Lindgren, H. Rue, and J. Lindstr¨ om, An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation ap- proach, J. R. Stat. Soc. B73, 423 (2011)

  17. [17]

    K. S. Cranmer, Kernel estimation in high-energy physics, Comput. Phys. Commun.136, 198 (2001)

  18. [18]

    H. P. Dembinski and A. S. W. Abdelmotteleb, A new maximum-likelihood method for template fits, Eur. Phys. J. C82, 1043 (2022)

  19. [19]

    C. E. Rasmussen and C. K. I. Williams,Gaussian Pro- cesses for Machine Learning(MIT Press, 2006)

  20. [20]

    de Boor,A Practical Guide to Splines, revised ed

    C. de Boor,A Practical Guide to Splines, revised ed. (Springer, 2001)

  21. [21]

    P. H. C. Eilers and B. D. Marx, Flexible Smoothing with B-splines and Penalties, Stat. Sci.11, 89 (1996)

  22. [22]

    Kimeldorf and G

    G. Kimeldorf and G. Wahba, A Correspondence Be- tween Bayesian Estimation on Stochastic Processes and Smoothing by Splines, Ann. Math. Stat.41, 495 (1970)

  23. [23]

    H. P. Dembinskiet al.,iminuit: A Python interface to MINUIT (2020)

  24. [24]

    A. L. Read, Presentation of search results: The CLs tech- nique, J. Phys. G28, 2693 (2002)

  25. [25]

    Trial factors for the look elsewhere effect in high energy physics

    E. Gross and O. Vitells, Trial factors for the look else- where effect in high energy physics, Eur. Phys. J. C70, 525 (2010), arXiv:1005.1891 [hep-ex]

  26. [26]

    Gandrakota, A

    A. Gandrakota, A. Cueto, and E. Gross, GP-based estimation of the look-elsewhere effect, (2022), arXiv:2206.12328 [physics.data-an]

  27. [27]

    M. K. Titsias, Variational Learning of Inducing Variables in Sparse Gaussian Processes,AISTATS 2009, Proc. Mach. Learn. Res.5, 567 (2009)

  28. [28]

    A. G. Wilson and H. Nickisch, Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP), ICML 2015, Proc. Mach. Learn. Res.37, 1775 (2015)

  29. [29]

    R. Brun, F. Rademakers,et al., ROOT: An Object- Oriented Data Analysis Framework (2024)

  30. [30]

    P. J. Bickel, C. A. J. Klaassen, Y. Ritov, and J. A. Well- ner,Efficient and Adaptive Estimation for Semiparamet- ric Models(Johns Hopkins Univ. Press, 1993)

  31. [31]

    R. D. Cousins, J. T. Linnemann, and J. Tucker, Eval- uation of three methods for calculating statistical sig- nificance when incorporating a systematic uncertainty into a test of the background-only hypothesis for a Pois- son process, Nucl. Instrum. Meth. A595, 480 (2008), arXiv:physics/0702156

  32. [32]

    Neyman and E

    J. Neyman and E. L. Scott, Consistent Estimates Based on Partially Consistent Observations, Econometrica16, 1 (1948)

  33. [33]

    S. A. Murphy and A. W. van der Vaart, On profile like- lihood, J. Amer. Statist. Assoc.95, 449 (2000)

  34. [34]

    W. K. Newey, The Asymptotic Variance of Semiparamet- ric Estimators, Econometrica62, 1349 (1994)

  35. [35]

    Chernozhukov, D

    V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins, Double/debiased machine learning for treatment and structural parame- ters, Econom. J.21, C1 (2018)

  36. [36]

    P. D. Dauncey, M. Sherwood, S. Sherwood, and A. sherwood, Handling uncertainties in background shapes: the discrete profiling method, JINST9, P04005, arXiv:1408.6865. Appendix A: Posterior calibration The GP posterior provides a point estimate and an uncertainty band for the true template rate. For these to be useful in a likelihood, the uncertainty must ...

  37. [37]

    In the stat-only case the eigenvalues decay smoothly fromλ 1 = 0.30 across the 40 bins

    Eigenvalue spectrum and mode shapes Figure 7 shows the eigenvalue spectrum of the GP pos- terior covariance for Experiment B, with and without systematics. In the stat-only case the eigenvalues decay smoothly fromλ 1 = 0.30 across the 40 bins. Adding the four systematic directions inserts a dominant calibration mode at ˜λ1 = 2.50 followed by three modes o...

  38. [38]

    Figure 9 maps the truncation-versus-coverage tradeoff on Dataset B, fitting the eigenmode model atk∈ {1,2,3,4,6,8,10,15,20,40}on 300 pseudo-experiments perkatµ true = 1

    Mode-truncation behaviour and full-rank equivalence The analytical argument of Section III C predicts that full-rank GP eigenmodes reproduce Barlow–Beeston in- ference. Figure 9 maps the truncation-versus-coverage tradeoff on Dataset B, fitting the eigenmode model atk∈ {1,2,3,4,6,8,10,15,20,40}on 300 pseudo-experiments perkatµ true = 1. Inside the plateau...

  39. [39]

    The parameter of interest isµ

    The Poisson semiparametric model A binned signal-plus-background analysis is a semi- parametric inference problem in the sense of Bickel, Klaassen, Ritov, and Wellner [30]. The parameter of interest isµ. The nuisance is the background shape g= (g 1, . . . , gN), anN-vector whose true value is un- known. The observed data and the MC auxiliary mea- surement...

  40. [40]

    The semiparametric information bound The Fisher information matrix of the model (D1)–(D2) at the true parameter values (µ0,g 0) has block structure I= Iµµ I T µg Iµg Igg ! ,(D4) with entries Iµµ = X j s2 j νj ,[I µg]j = sj νj , [Igg]jk =δ jk 1 νj + τ gj ,(D5) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Signal strength ¹ 0.0 0.2 0.4 0.6 0.8 1.0CLs CLs = 0:05 Expected §2¾ ...

  41. [41]

    (D3) maximisesℓ(µ,g) jointly overµandg

    Barlow–Beeston achieves the bound The profile likelihood estimator in the full model of Eq. (D3) maximisesℓ(µ,g) jointly overµandg. At each 14 fixedµ, the first-order condition∂ℓ/∂g j = 0 yields the profiled background Bj ≡n j +a j −(1+τ)µs j , ˆgj(µ) = Bj + q B2 j + 4(1+τ)τ µs jaj 2(1+τ) , (D9) which is the closed-form expression of Conway [8]. Sub- stit...

  42. [42]

    signal direction

    The GP plug-in as a two-step estimator The GP template is constructed in two steps. First, a GP is fitted to the MC countsa, yielding a smooth background estimate ˆgGP j =e ˆf(x j)wj/τ. Second, this estimate is plugged into the data likelihood andµis es- timated from ˆµplug = arg max µ X j nj ln(µsj + ˆgGP j )−(µs j + ˆgGP j ) . (D10) The background ˆgGP ...

  43. [43]

    The resulting estimator ˆµDML is unbiased and √n-consistent

    DML: unbiasedness without efficiency The DML framework of Chernozhukovet al.[35] con- structs an explicit Neyman-orthogonal moment condition by adding a correction that exploits the unbiased his- togram estimator ˜aj =a j/τ: ψorth = X j sj nj νj −1 + X j sj νj ˆgGP j −˜aj .(D13) Under a perturbation ˆgGP =g+δg, the bias from the first sum and the correcti...

  44. [44]

    Their discrete profiling method treats the function-form index as a discrete nuisance parameter, profiled alongside µby taking the likelihood envelope across candidate func- tions

    Connection to discrete profiling Daunceyet al.[36] addressed a related problem: un- certainty in the functional form of a smooth background. Their discrete profiling method treats the function-form index as a discrete nuisance parameter, profiled alongside µby taking the likelihood envelope across candidate func- tions. The coverage and bias of this metho...

  45. [45]

    The first quantifies how much of the parameter-of-interest information remains after the nui- sance has been profiled (Fig

    Quantifying the significance deficit The semiparametric analysis is more than a structural caveat: it decides which method is appropriate at which τ, and the quantitative answers come from two comple- mentary pictures. The first quantifies how much of the parameter-of-interest information remains after the nui- sance has been profiled (Fig. 12). The secon...

  46. [46]

    Regime summary The analysis above delineates where each method has the advantage. In the statistically limited, single-channel regime, BB’s joint profiling achieves the semiparametric efficiency bound and the GP plug-in does not; the per- bin gamma factors are not wasted parameters but the mechanism of efficient inference. In the systematically limited re...