Determination of proton PDF uncertainties with Markov Chain Monte Carlo

TL;DR

This work uses Markov Chain Monte Carlo (MCMC) to sample the full posterior distribution of a 15-parameter proton PDF model, enabling principled Bayesian uncertainty quantification and direct propagation to observables. By fitting to a broad dataset of deep inelastic scattering and Drell–Yan data with approximate NNLO theory, the authors compare MCMC-derived uncertainties to the Hessian method, revealing substantial non-Gaussian features in several PDF components and demonstrating how MCMC can define a statistically meaningful tolerance via the $\\chi^2$ distribution. The study shows that, when Gaussianity holds, MCMC and Hessian results largely agree, but in non-Gaussian regimes the Hessian approach can misestimate uncertainties, highlighting the value of a posterior-based analysis for percent-level precision SM tests. The paper also discusses practical limitations of the current MCMC implementation and outlines future directions, including potential adoption of Hamiltonian Monte Carlo to improve sampling efficiency for high-dimensional PDF fits.

Abstract

We present an analysis of parton distribution functions (PDFs) of the proton using Markov Chain Monte Carlo (MCMC) methods. The MCMC approach naturally implements Bayes' theorem and thus provides a means to directly sample the underlying probability distribution - in this case the probability distribution of the PDF parameters. This allows for a straightforward propagation of the resulting uncertainties into any PDF-dependent observable, preserving their simple probabilistic interpretation. In our analysis we include a broad set of deep inelastic scattering data from HERA, BCDMS and NMC experiments along with the Drell-Yan, $W$ and $Z$ boson data from LHC and Tevatron experiments, which combined with theoretical calculations at next-to-next-to-leading order in QCD allow for realistic determination of PDFs. The main focus of this analysis is to explore alternative methods for PDF uncertainty estimation that are more firmly grounded in statistical principles. We show that the flexibility of the Bayes framework, allowing e.g. to account for non-Gaussianity or inconsistencies of data sets, is crucial to extract realistic uncertainties when such assumptions are not fulfilled. We also demonstrate that MCMC allows one to determine the $Δχ^2$ value corresponding to a given confidence level in the sample, which can in turn be used as a statistically well-founded tolerance criterion used in the Hessian method, thus addressing one of its main long-standing drawbacks.

Figures (9)

• Figure 1: Kinematic coverage of the experimental data sets in the $(x,Q^2)$-plane. The data sets are displayed after cuts and leading order kinematics have been assumed for the DY data sets. The gray area indicates the global $Q^2$- and $W^2$-cut.
• Figure 2: The convergence of six of the 36 generated chains towards the same minimum. The left plot shows the $(\chi^2-\chi^2_{min})$-value for each sample and the right plot the parameter $p_4^{u_v}$ versus the Monte Carlo time. The dashed lines indicate the times when a reset was performed and the grey regions the steps when the static/fixed proposal distribution was used. Note that only the first 160 000 points are shown (for the remaining 319 000 points the behaviour is similar to the last shown part).
• Figure 3: The samples/marginal distribution for the parameter $p_4^{d_v}$ (left axis) along with a one-dimensional $\chi^2$-scan (right axis). The gray-shaded areas indicate the regions excluded by the flat prior.
• Figure 4: The integrated autocorrelation time $\tau_{\rm int}$ as a function of the thermalization time $t_{\rm therm}$ (after thinning with $\eta=3\,000$, motivated in the next section). At $t_{\rm therm}>120\,000$ the autocorrelation time exhibits a plateau. Repeating this exercise for each parameter individually, we find the same plateau starting in the range of $t_{\rm therm}\in[50\,000,130\,000]$.
• Figure 5: The integrated autocorrelation time $\tau_{\rm int}(W)$ (upper plot, see \ref{['eq:tau_int']}) and the autocorrelation function $\rho(W)$ (lower plot, see \ref{['eq:rho_W']}) as a function of the summation window $W$ in units of Monte Carlo time. In principle $\tau_{\rm int}$ is defined by summing over $\rho(W)$ until infinity. Since this is impossible we resort to the $\Gamma$-method that defines an algorithm to find the optimal value of $W$ to truncate the summation. This is referred to as the optimal summation window, $W_{\rm opt}$, indicated by the vertical red dashed lines. The horizontal red dashed lines show the estimated $\rho(W_{\rm opt})$ and $\tau_{\rm int}(W_{\rm opt})$. Going from left to right we have thinned the original chain with factors $\eta=\{500,1500,3000\}$ and thus reducing $\tau_{\rm int}$.