Propagating data noise through the fit: the Monte Carlo replica distribution

Technical Summary: Propagating Data Noise through the Fit: The Monte Carlo Replica Distribution

Problem Statement
The Monte Carlo (MC) replica method is a standard procedure in high-energy physics for quantifying parameter uncertainties in global fits of Parton Distribution Functions (PDFs) and Standard Model Effective Field Theory (SMEFT) Wilson coefficients. The method involves generating many noise-perturbed copies of experimental data (pseudodata) and performing independent fits to each. The resulting empirical distribution of the best-fit parameters is used as the uncertainty estimate.

While it is well-established that this method reproduces the Bayesian posterior distribution exactly when the theoretical model is linear in its parameters, it is known to deviate in the nonlinear case. However, the nature and magnitude of this deviation have not been characterized in a closed form. Without such a characterization, practitioners cannot easily determine whether the bare MC replica distribution is a faithful approximation of the Bayesian posterior or if computationally expensive corrective sampling (e.g., Metropolis–Hastings) is required to recover the exact posterior.

Methodology
The paper treats the Maximum Likelihood Estimator (MLE) not merely as an optimization result, but as a deterministic map $\hat{\theta}(d)$ from data space to parameter space. The MC replica distribution is identified as the pushforward of the Gaussian noise $\eta$ under the map $\eta \mapsto \hat{\theta}(D_0 + \eta)$, where $D_0$ is the observed data.

To derive the distribution analytically, the authors perform a first-order Taylor expansion of the MLE map $\hat{\theta}(\cdot)$ around the observed data $D_0$. This reduces the pushforward to a linear-Gaussian map, the covariance of which is determined by the Jacobian of the MLE map, $\partial \hat{\theta} / \partial d$, evaluated at $D_0$.

The core technical innovation is the derivation of this Jacobian in closed form using implicit differentiation of the first-order optimality condition (the vanishing gradient of the loss function). By defining the loss $L(\theta; d) = \frac{1}{2}\|d - T(\theta)\|^2$ and its gradient $F(\theta, d) = J(\theta)^\top(T(\theta) - d) = 0$, the authors differentiate the identity $F(\hat{\theta}(d), d) \equiv 0$ with respect to $d$. This yields the Jacobian:
$$\frac{\partial \hat{\theta}}{\partial d}(D_0) = A_\star^{-1} J_\star^\top$$
where $J_\star$ is the theory Jacobian at the best fit, and $A_\star$ is the Hessian of the loss function at the best fit. Crucially, $A_\star$ is decomposed into the Gauss–Newton matrix ($J_\star^\top J_\star$) and a "residual-weighted Hessian" term involving the second derivatives of the theory map weighted by the residuals.

Key Contributions and Results

1. Closed-Form Leading-Order Distribution: The paper derives a closed-form expression for the leading-order MC replica covariance matrix:
   $$\Sigma_{\text{MC}} = \sigma^2 A_\star^{-1} J_\star^\top J_\star A_\star^{-1}$$
   (assuming isotropic noise $C = \sigma^2 I$).

2. Comparison with Bayesian Posterior: The authors compare this MC covariance with the Laplace approximation of the Bayesian posterior, which is Gaussian with covariance $\Sigma_{\text{Bayes}} = \sigma^2 A_\star^{-1}$.

3. The Discrepancy Matrix: The difference between the two distributions is shown to be a single computable matrix:
   $$\Sigma_{\text{MC}} - \Sigma_{\text{Bayes}} = -\sigma^2 A_\star^{-1} \left( \sum_a r_a H_{a,\cdot\cdot}(\hat{\theta}_\star) \right) A_\star^{-1}$$
   where $r_a$ are the residuals and $H_{a,\cdot\cdot}$ are the components of the theory Hessian. This term, the residual-weighted Hessian, is the sole source of deviation between the MC method and Bayesian inference at leading order.

4. Regimes of Agreement and Disagreement:

   • Linear Theory: The Hessian term vanishes, and the MC distribution coincides exactly with the Bayesian posterior.
   • Nonlinear Theory, Good Fit: If residuals are small, the discrepancy is negligible, and the MC method is a faithful approximation.
   • Nonlinear Theory, Poor Fit: If residuals are significant and the theory has curvature, the MC method may over- or under-estimate uncertainties depending on the sign of the residual-weighted Hessian.

5. Toy Examples: The paper validates the derivation using two single-parameter examples:

   • Quadratic Theory: Demonstrates that in the "unreachable" regime (where data lies outside the theory's range), the MC distribution collapses to a delta function (zero variance) because the Jacobian vanishes, whereas the Bayesian posterior remains a non-degenerate Gaussian.
   • Circular Theory: Illustrates how the MC method can over-cover (when data is inside the circle) or under-cover (when data is outside) relative to the Bayesian posterior, consistent with previous numerical findings.

Significance and Claims
The paper claims to provide a "complementary derivation" that makes the nonlinear departure of the MC replica method from Bayesian inference explicit and computable. The primary significance lies in the derivation of the residual-weighted Hessian as a diagnostic tool.

The authors argue that evaluating this matrix is computationally cheap compared to running corrective sampling schemes. Consequently, this expression allows practitioners to determine a priori whether the bare MC replica distribution is sufficient for their required accuracy or if the discrepancy is large enough to necessitate more rigorous Bayesian sampling. The paper positions this as a step toward understanding the faithfulness of MC methods in complex fits (like PDFs and SMEFT) without necessarily replacing them, but rather providing a metric for their validity. The authors explicitly state that evaluating this matrix in full PDF and SMEFT fits is a "natural next step," but do not present such a full-scale application in this work.