  We present a detailed mathematical study of the Monte Carlo replica method as
applied in the global fitting literature from the high-energy physics theory
community. For the first time, we provide a rigorous derivation of the
parameter distributions implied by the method, and show that, whilst they agree
with Bayesian posteriors for linear models, they disagree otherwise. We proceed
to numerically quantify the disagreement between the Monte Carlo replica method
and the Bayesian method in the context of two phenomenologically relevant
scenarios: fits of the SMEFT Wilson coefficients, and fits of PDFs (albeit in a
toy scenario). In both scenarios, we find that uncertainty estimates of the
quantities of interest are discrepant between the two approaches when
non-linearity is relevant. Our findings motivate future investigation of
Bayesian methodologies for global PDF fits, especially in the context of
simultaneous determination of PDFs and SMEFT Wilson coefficients.
