#  A generalized statistical model for fits to parton distributions

Mengshi Yan  Department of Physics and State Key Laboratory of Nuclear Physics
and Technology, Peking University, Beijing 100871, China<EMAIL_ADDRESS>Tie-Jiun Hou  School of Nuclear Science and Technology, University of South
China, Hengyang 421001, China<EMAIL_ADDRESS>Zhao Li  Institute of High
Energy Physics, Chinese Academy of Sciences, Beijing 100049, China;
School of Physics Sciences, University of Chinese Academy of Sciences, Beijing
100039, China  Center for High Energy Physics, Peking University, Beijing
100871, China<EMAIL_ADDRESS>Kirtimaan Mohan<EMAIL_ADDRESS>C.–P. Yuan
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, U.S.A<EMAIL_ADDRESS>
###### Abstract

Parton distribution functions (PDFs) form an essential part of particle
physics calculations. Currently, the most precise predictions for these non-
perturbative functions are generated through fits to global data. A problem
that several PDF fitting groups encounter is the presence of tension in data
sets that appear to pull the fits in different directions. In other words, the
best fit depends on the choice of data set. Several methods to capture the
uncertainty in PDFs in presence of seemingly inconsistent fits have been
proposed and are currently in use. These methods are important to ensure that
uncertainty in PDFs are not underestimated. Here we propose a novel method for
estimating the uncertainty by introducing a generalized statistical model
inspired by unsupervised machine learning techniques, namely the Gaussian
Mixture Model (GMM). Using a toy model of PDFs, we demonstrate how the GMM can
be used to faithfully reconstruct the likelihood associated with PDF fits,
which can in turn be used to accurately determine the uncertainty on PDFs,
especially in presence of tension in the fitted data sets. We further show how
this statistical model reduces to the usual chi-squared likelihood function
for a consistent data set and provide measures to optimize the number of
Gaussians in the GMM.

††preprint: MSUHEP-24-002

###### Contents

  1. I Introduction
  2. II Tension in a toy model
    1. II.1 Two pseudo-data sets with tension
    2. II.2 Least Squares Fits and Tension
  3. III Gaussian mixture model in PDF fit
    1. III.1 Performing parameter fits with the GMM
    2. III.2 Mean and Uncertainty in the GMM
  4. IV Application of the GMM in a toy model of PDFs
    1. IV.1 Applying GMM to different scenarios
    2. IV.2 The GMM method compared to different tolerance criteria
  5. V Consistency checks for the GMM
    1. V.1 Consistent results from consistent data
    2. V.2 Information Criteria for an Optimal number of base Gaussians
  6. VI Conclusion

##  I Introduction

Parton distribution functions (PDFs) are non-perturbative functions
quantifying the probability of finding a quark or gluon with certain momentum
fraction $x$ from a nucleon. Taking advantage of a wide variety of
experimental data and advanced computating packages for perturbative QCD (and
electroweak) calculations, various groups [1, 2, 3, 4, 5] provide a
sophisticated determination of PDFs. The recently released NNPDF [6] and
xFitter [7] fitting codes open opportunities for general users to also conduct
global QCD analyses.

These PDFs form a critical component of the theory community’s infrastructure
[8] and developing the next generation of PDFs are crucial for achieving the
physics goals of the Large Hadron Collider (LHC) and the Electron-Ion Collider
(EIC). An unbiased determination of PDFs are essential for precision
phenomenology and the search for new physics, such as the recent measurement
of $W$ boson mass by the CDF collaboration [9]. A robust understanding of the
the uncertainties arising from PDFs is necessary to provide precise
predictions. PDF uncertainties constitute one of the dominant systematic
errors both for direct searches of new physics [10] and for precision tests of
the Standard Model (SM), such as determining the properties of the Higgs boson
[11]. A faithful determination of PDFs and its uncertainties are thus
necessary in order to make precise theory predictions at hadron-hadron and
lepton-hadron colliders. Improving PDFs will require new experimental data,
new theoretical (higher-order) calculations and a determination of fits to
global data in a statistically robust way [8].

Many issues could be encountered when trying to perform global fits to data
sets from a variety of experiments. One of these issues, which is the primary
focus of this paper, is in dealing with data that appear to be in tension with
each other, leading to inconsistencies within the global fit. For example,
suppose we have two data sets that pull the fit parameters to different
values. In this scenario, the data are in tension with each other. The source
of such tension could arise due to deficiencies in theoretical models and
calculations or even from underestimation of experimental uncertainties. To
keep the discussion simple, we work with an example of the latter, although,
the method presented here is agnostic to the source of tension.

An example of tension between data sets would be the measurement of the ATLAS
8 TeV $W$ and $Z$ rapidity distribution [12], which was found to be
inconsistent with other data included in the PDF global analysis and favored
an enhanced strange quark distribution function [1, 2, 13]. Although a non-
vanishing strangeness asymmetry $(s-\bar{s})(x)$ can partially release this
tension [2, 13], it is still inadequate to accommodate the ATLAS 8 TeV $W/Z$
data [12] consistently in a PDF global analysis.

Another example is the inconsistency between the E866 NuSea and E906 SeaQuest
measurements of the cross-section ratio $\sigma(pd)/2\sigma(pp)$, which
approximates the ratio of anti-quark PDFs
$\sigma(pd)/2\sigma(pp)\approx\big{(}1+\bar{d}_{p}(x_{2})/\bar{u}_{p}(x_{2})\big{)}/2$.
The E866 NuSea and E906 SeaQuest measurements provide valuable information on
the PDF ratio $(\bar{d}/\bar{u})(x)$ in the large-$x$ region. However, they
are found to be in tension at large-$x$ [14, 15, 13]. The effect of tension
between certain data sets generally becomes negligible for large data sets,
when outliers appear with small probability. However, in this case, E866 NuSea
and E906 SeaQuest, the measurements probe a specific kinematic region that is
otherwise weakly constrained and hence the tension between them cannot be
ignored. It is thus necessary to develop a prescription to accommodate such
inconsistencies in a PDF global analysis in order to properly estimate the PDF
uncertainty.

Several methods have been proposed and are currently in use to provide an
estimate of uncertainties especially in the presence of tension[16]. One naive
approach would be to simply omit one of the troublesome data sets or perhaps
attach a significantly lower weight to this data set. This approach introduces
biases in the fitting procedure that needs to be justified with some prior
knowledge. In some cases, such a justification may be possible. For example,
if perhaps the methodology of the omitted experiment is analyzed and found to
be unsatisfactory. In such a scenario, the prior bias in omitting the data set
from this experiment is justified, but this is not always the case, since it
is not always straightforward to reanalyze particle physics experiments and
their data.

Modern PDFs perform fits to data by defining a likelihood or loss function.
For PDFs produced by CTEQ-TEA (CT) [1] and MSHT [2], the likelihood is the
$\Delta\chi^{2}$. Under normal circumstances with a consistent set, the
uncertainty can be defined by the $68\%$ confidence level limit of the
$\chi^{2}$ distribution, which amounts to varying the fit and determining
boundaries where $\Delta\chi^{2}=1$. The advanced methodology adopted by both
CTEQ-TEA and MSHT is to introduce a tolerance $T$ such that the fits should be
varied within $\Delta\chi^{2}=T^{2}$. Values of $T$ are chosen either
dynamically [2] or set to a fixed value [1]. Typical values of $T$ are much
larger than $1$, which increases the uncertainties in the PDFs. NNPDF uses the
Monte-Carlo method for determining uncertainties, where they produce
representative sample of PDFs by sampling their likelihood or loss function
with the help of pseudo-data,which might be one of the reasons for the
uncertainties from NNPDF’s analysis being notably different [17]. Although, as
pointed out in Refs. [18, 19, 20], the level of inconsistency with data sets
is found to be negligible in the NNPDF3.1 analyses, further study is needed to
understand differences between NNPDF and the other fitting groups. One such
task was performed in Ref. [21], which details benchmark comparisons of three
global PDF sets - CT18, MSHT20 and NNPDF3.1 \- and their similarities and
differences that have been observed. The end result of that study was a new
PDF4LHC21 ensemble of combined PDFs suitable for a wide range of LHC
applications.

In this paper we describe how to apply a more general statistical framework to
accommodate potential inconsistencies in PDF data. Specifically, we show how
the Gaussian Mixture Model (GMM), which is commonly used as an unsupervised
machine learning technique, can be adapted to produce PDF uncertainties that
faithfully represent the likelihood of PDF fits. Our major findings are
summarized below for the benefit of the reader.

  * •

The Least Squares(LS) method is based on the assumption that the likelihood
function looks like a single Gaussian distribution. Thus the likelihood is not
faithfully represented when dealing with inconsistent data sets. The resulting
uncertainties can be highly overestimated in some regions whereas they can
simultaneously be highly underestimated in other regions.

  * •

The GMM provides a way to faithfully represent the likelihood of a combination
of data sets with tension, which in turn gives a more accurate estimate of
uncertainties. Therefore, the implementation of GMM in the PDF global fitting
can be expected to accommodate the data sets with varying degrees of tension.

  * •

It is necessary to assess the level of tension in the data sets and to
determine the optimal hyperparameters of the GMM. We use the Akaike
Information Criterion (AIC) and Bayesian Information Criterion (BIC) to make
such a determination and avoid either over-fitting or under-fitting.

The rest of the paper is organized as follows. In Sec. II, we set up a toy PDF
model and describe different pseudo-data sets, some of which have obvious
tension between them. Section III introduces the Gaussian Mixture Model as the
statistical model to be used in order to deal with inconsistencies in data
sets. In Sec. IV, the results of the GMM fit are compared with those of
conventional fits using the usual least squares fit. Consistency between the
GMM and the commonly-used Least Squares (LS) method is evaluated using a data
set with no tension in Sec. V. Finally, we discuss our results and conclude in
Sec. VI.

##  II Tension in a toy model

In this section we introduce a toy model of PDFs and perform LS fits in order
to illustrate issues that arise when there is tension between data sets.

###  II.1 Two pseudo-data sets with tension

We begin by utilizing the toy model of pseudo-PDF suggested in Ref. [22] and
choose the functional form of the underlying truth to be,

$g(x)=a_{0}\ x^{a_{1}}(1-x)^{a_{2}}e^{xa_{3}}(1+xe^{a_{4}})^{a_{5}}.$ (1)

Here $\left\\{a_{0},a_{1},a_{2},a_{3},a_{4},a_{5}\right\\}$ are parameters of
the function. To keep the discussion clear we will only perform fits to
subsets of these parameters, while keeping others fixed. We would like to
emphasize that although we choose an explicit functional form, the starting
point could have also been a deep neural network or a polynomial
parametrization.111The use of a deep neural network is supported by the
Universal Approximation Theorem [23] whereas a polynomial is supported by the
Stone-Weierstrass Theorem [24, 25]. These theorems ensure that any smooth
function can be well approximated either by neural networks or by polynomials.
Here, we restrict ourselves to use only the functional form in Eq. (1) for
clarity of presentation.

We generate pseudo-data sets by setting values of the parameters $a_{i}$ to
certain fixed values. Then, for each fixed value of $x$ we generate a pseudo-
data point by assuming the data follows a normal distribution about the truth
with uncertainty given by

$\Delta g(x)=\frac{\alpha}{\sqrt{g(x)}}.$ (2)

Here, the parameter $\alpha$ in Eq. (2) is used to scale the size of
uncertainty. In the region where the magnitude of pseudo-PDF is small, the
size of uncertainty is large, mimicking the case of real world data which has
larger statistical uncertainty in regions where event numbers are small. We
generate pseudo-data sets using

$g_{D}(x)=\Big{(}1+r\times\Delta g(x)\Big{)}g(x),$ (3)

where the random variable $r$ follows a Gaussian distribution with unit
standard deviation. Although, it is possible to include correlated systematic
uncertainties or even full likelihoods from experimental data, for simplicity,
we treat the pseudo-PDF uncertainty to be uncorrelated and purely statistical.

For the following analysis, we create six pseudo-data sets that are summarized
later in Table 4. In this section we describe only two of those pseudo-data
sets, which have different underlying truth values and thus are in tension
with each other. We refer to these two data sets as pseudo-data $\\#1$ and
pseudo-data $\\#2$. These pseudo-data sets are generated on the same array of
$x$ values but with different values of parameters $a_{i}$ as shown in Table
1. The comparison of these two pseudo-data sets, along with that of their
respective pseudo-PDF truths, is shown in Fig. 1. The tension between the
pseudo-data sets is strong in the range $10^{-2}<x<0.1$, and less significant
in other ranges of $x$. Generating the pseudo-data sets for the same values of
$x$ ensures that the tension is explicit without any wiggle-room to smooth
over tension with a flexible parameterization of pseudo-PDFs or due to a
different theoretical model. It is important to note that, in general,
inconsistencies in data might be a result of an inadequate theoretical model,
such as underestimation of theoretical uncertainties or a less flexible
parametrization of non-perturbative PDFs, or from an under-estimation of the
experimental data error. In either case, the GMM method we propose here, will
be able to capture uncertainties from both theoretical and experimental.
Hence, we do not concern ourselves with the sources of data tension, but
rather, we investigate methods to estimate uncertainty when fitting
inconsistent data with a fixed theoretical prescription.

| $N_{\text{pt}}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$
---|---|---|---|---|---|---|---
pseudo-data #1 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.4 | -3.0
pseudo-data #2 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.6 | -2.8
Table 1: The choice of parameters ($a_{i}$) of the function $g(x)$ defined in
Eq.(1) used to generate the two pseudo-data sets #1 and #2. $N_{pt}$ refers to
the number of data points in the set.

Figure 1:  Two pseudo-data sets, $g_{D}(x;\\#1)$ (black) and $g_{D}(x;\\#2)$
(red) are in tension with each other. The underlying truth $g(x)$ for
individual pseudo-data sets are also shown by the solid curves which are
referred to as truth $g(x;\\#1)$ and truth $g(x;\\#2)$. The plot on the left
shows the pseudo-PDFs whereas the plot on the right shows their ratio with
respect to the underlying truth $g(x;\\#1)$.

###  II.2 Least Squares Fits and Tension

Before introducing the GMM method we first illustrate the method of dealing
with inconsistent data in the Least Squares fitting framework. Fits are
performed by optimizing the log likelihood, which in our case is the
$\chi^{2}$ function defined as

$\chi^{2}(\theta)=\sum_{j=1}^{N_{\text{pt}}}\bigg{(}\frac{y_{i}-T(\theta)}{\Delta
y_{i}}\bigg{)}^{2}.$ (4)

Given the known $N_{\text{pt}}$ measurements $y_{i}$ associated with
uncertainties $\Delta y_{i}$, the $\chi^{2}$ provides a normalized distance
from theoretical predictions $T(\theta)$ to the measurements. Here $\theta$
are the fit parameters, which are some subset of the parameters $a_{i}$, and
are varied in order to optimize the $\chi^{2}$ function.

To illustrate the problem with the LS fit when tension exists, three parallel
LS fits are performed, which we refer to as SG-A (performed only with pseudo-
data #1), SG-B (performed only with pseudo-data #2) and SG-C (performed by
combining both pseudo-data #1 and pseudo-data #2). We label these fits as “SG”
to refer to the single-Gaussian nature of the likelihood in LS method, in
contrast to the multi-modal distribution that we will introduce later. Results
of these fits are shown in Table 2. In these fits, we use the true functional
form in Eq. (1), for the theory calculation, and fix parameters
$\\{a_{0},a_{1},a_{2},a_{3}\\}$ to their true values which can be found in
Table 1. The parameters $a_{4}$ and $a_{5}$ are varied and the values that
minimize the $\chi^{2}$ function are shown in Table 2. The values of
$\chi^{2}_{\\#1}/N_{\text{pt}}$, $\chi^{2}_{\\#2}/N_{\text{pt}}$, and
$\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$ indicate the goodness of fit to
individual pseudo-data sets # 1, # 2 and the combined pseudo-data set
respectively. A good fit will have values of a $\chi^{2}/N_{\text{pt}}\sim 1$,
whereas a poor fit has $\chi^{2}/N_{\text{pt}}$ values significantly larger
than 1. The best-fit of SG-A, where only pseudo-data set # 1 is included, is
found at the minimum of $\chi^{2}_{\\#1}/N_{\text{pt}}=0.88$, with values of
$a_{4}=2.32$ and $a_{5}=-3.22$. Also shown in the same row for SG-A are values
of $\chi^{2}_{\\#2}/N_{\text{pt}}$ and $\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$,
these indicate the values of each of the likelihood functions
$\chi^{2}_{\\#2}/N_{\text{pt}}$ and $\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$
determined at the best fit values of $a_{4}$ and $a_{5}$ for fits SG-A. We see
that for SG-A the values of $\chi^{2}_{\\#2}/N_{\text{pt}}=6.55$ and
$\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}=3.72$ are significantly larger than one,
indicating that there is significant tension between both the combined data
set and the data set #2.

Likewise, the best-fit of SG-B, where only pseudo-data set # 2 is included, is
found at the minimum of $\chi^{2}_{\\#2}/N_{\text{pt}}=1.02$, indicating a
reasonably good fit. For SG-B, the values $\chi^{2}_{\\#1}/N_{\text{pt}}$ and
$\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$ are determined by using the best fit
value of $a_{4}=2.63$ and $a_{5}=2.72$. Once again large values of
$\chi^{2}_{\\#1}/N_{\text{pt}}=7$ and
$\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}=4.01$ indicate tension with both the
combined data set and data set #1. Recall, that tension between the two data-
sets, pseudo-data #1 and pseudo-data #2 are generated by starting with
different choices of $a_{4}$ and $a_{5}$ as given in Table 1. It can be seen,
as expected, that the best-fit values of $a_{4}$ and $a_{5}$ in the two fits
without tension, namely SG-A and SG-B, are indeed close to their underlying
truth values.

The fit, SG-C, to the the combined data set where both pseudo-data sets # 1
and # 2 are included, is found at the minimum of
$\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}=2.42$, indicating, as expected, a poor fit
to the combined data set. In this case, the values of
$\chi^{2}_{\\#1}/N_{\text{pt}}=2.27$, $\chi^{2}_{\\#2}/N_{\text{pt}}=2.42$
indicate that the combined fit, still has tension with both individual data
sets.

fits | pseudo-data | $N_{\text{pt}}$ | best-fit $a_{4}$ | best-fit $a_{5}$ | $\chi^{2}_{\\#1}/N_{\text{pt}}$ | $\chi^{2}_{\\#2}/N_{\text{pt}}$ | $\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$
---|---|---|---|---|---|---|---
SG-A | # 1 | 50 | 2.32 | -3.22 | 0.88 | 6.55 | 3.72
SG-B | # 2 | 50 | 2.63 | -2.73 | 7.00 | 1.02 | 4.01
SG-C | # 1 and # 2 | 100 | 2.48 | -2.94 | 2.27 | 2.56 | 2.42
truth | # 1 | - | 2.4 | -3.0 | - | - | -
truth | # 2 | - | 2.6 | -2.8 | - | - | -
Table 2:  The table shows the best-fit values of $a_{4}$ and $a_{5}$, as well
as a measure of the goodness of fit $\chi^{2}/N_{\text{pt}}$ for three Least
Square fits SG-A, SG-B, and SG-C, performed using pseudo-data sets # 1 and #
2. The columns $\chi^{2}_{\\#1}/N_{\text{pt}}$,
$\chi^{2}_{\\#2}/N_{\text{pt}}$ and $\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$
indicate the values of these log-likelihood functions evaluated at values of
$a_{4}$ and $a_{5}$ in the given rows. Also shown, for comparison, in the last
two rows are the true values of $a_{4}$ and $a_{5}$ that were used to generate
the two psuedo-data sets.

Although the fit to the combined data is poor, the main problem arises when
attempting to determine the relevant uncertainty associated with such a fit.
This is illustrated in Fig. 2, which presents Results of the LS fits SG-A,
SG-B and SG-C. For all the three fits, results are presented in the form of
replica sets, where the probability distribution of the pseudo-PDF is sampled
with the Monte Carlo method as described in [6]. For the purpose of
comparison, we also estimate the uncertainty of SG-C fit using the Hessian
method [1], which is consistent with the Monte Carlo method. Here we only
consider the Hessian uncertainty generated by the variation of $\chi^{2}$ and
determine the uncertainty with the criterion that $\Delta\chi^{2}=1$. Here, we
do not include any additional tolerance criterion, such as the two-tier
penalty for the CT uncertainty [1] and the dynamic tolerance for the MSHT
uncertainty [2]. Later, we will include such tolerance criteria for comparison
with the GMM.

In panels 2 and 2, pseudo-PDFs from fits SG-A, SG-B and SG-C are found to be
in disagreement in the range $10^{-3}<x<0.2$. Even the $1\sigma$ uncertainty
bands do not overlap in this region. For $x>0.2$, the tension in data is weak
due to larger uncertainty in the data and this can be seen from the overlap of
the uncertainty bands in this region. Even though the uncertainty bands
overlap at large $x$, the impact of significant tension in the range
$10^{-3}<x<0.2$ still affects pseudo-PDF at large $x$ as can be observed by
the different shapes of the uncertainty bands for the different LS fits. For
the LS fit SG-C, the Hessian uncertainty with $\Delta\chi^{2}=1$ is also
provided, which deviates from the Monte Carlo uncertainty by a negligible
amount. The tension between the different fits, SG-A, SG-B and SG-C, can also
be observed by looking at the result of the fits to the parameters $a_{4}$ and
$a_{5}$, shown in Fig. 2. As expected, the true values of the parameters
$a_{4}$ and $a_{5}$ fall within the regions covered by the Monte Carlo sample
of 1000 replicas for each of the fits SG-A and SG-B. However, the true values
of the fit parameters are far from any of the replicas generated with SG-C.
Furthermore, none of the replicas generated for either of the fits overlap and
to make matters worse, neither the $1\sigma$ nor $3\sigma$ uncertainty bands
of SG-C, generated using the Hessian method, overlap with the two truth values
for $a_{4}$ and $a_{5}$ shown in Table 1. The tolerance criteria used by the
PDF fitting groups would increase the uncertainty, but as to be shown in
Section IV, this approach will not be sufficient to incorporate both data sets
into a single fit.

One might wonder why the fit SG-C, which combines both data sets, is unable to
reproduce the desired uncertainty bands that should look like some combination
of the uncertainty bands of both SG-A and SG-B. Such a combination of
uncertainty bands would arise only if the likelihood function were multi-
modal. On the other hand, fit SG-C assumes that there is only a single mode
for the likelihood. In other words, we have ignored the two individual
likelihoods generated by the two data sets and tried to define a new
likelihood in SG-C that we have assumed to be uni-modal. Another way of
stating this is to say that a naive combination of Monte-Carlo samples from
SG-A and SG-B, seen in Fig. 2, is a better representation of the combined
likelihood of both data sets. This combined likelihood is multi-modal, unlike
that of SG-C. Of course, just by adding the Monte Carlo samples does not
suffice and we need a proper method to combine the samples. In order to
combine likelihoods in a statistically robust way, we introduce the GMM method
in the next section which provides a way to combine likelihoods into a multi-
modal distribution.

Before we conclude this section, we comment on alternate methods to deal with
such inconsistencies in fits. For example, the Particle Data Group (PDG) [26],
recommends increasing the uncertainty on data sets, either all data sets if
they are of similar precision or the data set from the experiment with the
highest precision, by a scale factor $S=\sqrt{\chi^{2}/N_{\text{pt}}}$. For
data sets with tension, $\sqrt{\chi^{2}/N_{\text{pt}}}>1$ and hence the scale
factor $S>1$, implying an increase in uncertainty. For a large data set
($N_{\text{pt}}\gg 1$), increasing the uncertainty on all data sets implies
that $\chi^{2}/N_{\text{pt}}\simeq 1$, thus ensuring a good fit. This idea has
been further explored and elucidated in both the Bayesian and Frequentist
approach [27, 28, 29]. The GMM presented here, not only provides an alternate
way to combine inconsistent data sets without increasing their uncertainty,
but also, if desired, provides an alternate estimate for the scale factor $S$
required to make data sets consistent. We describe some of these ideas briefly
in Sec. V.2.

Figure 2:  Panel (a): the comparison of pseudo-data and optimal pseudo-PDFs.
The uncertainties in replica sets correspond to the 1 standard deviation
interval. For the LS fit $C$, the Hessian uncertainty with $\Delta\chi^{2}=1$
is also provided, which deviates from the Monte Carlo uncertainty by a
negligible amount. The sizes of uncertainties are small in view of the
magnitudes of pseudo-PDFs. Panel (b): similar to panel (a), but for ratios to
the central value of the Hessian set of the SG-C fit. Panel (c): clouds of
fitted values of parameters $a_{4}$ and $a_{5}$ in replica sets, as compared
to true values.

##  III Gaussian mixture model in PDF fit

As mentioned in the earlier section, the likelihood from each of the data sets
(pseudo-data #1 and pseudo-data #2) are significantly different and needs to
be combined in a meaningful way. We propose the use of the Gaussian Mixture
Model to construct the combined likelihood for data with tension.

The GMM was first used by Karl Pearson in 1894 [30]. It is now used
extensively as an unsupervised machine learning technique to classify data
[31]. In general, the GMM is a mixture of $K$ base Gaussian distributions.
Given the data set $Y=\\{(y_{1},\Delta y_{1}),(y_{2},\Delta
y_{2}),\cdots,(y_{N_{\text{pt}}},\Delta y_{N_{\text{pt}}})\\}$, where $\Delta
y_{j}$ is the uncertainty associated with the measured value $y_{j}$222 For
clarity of presentation, we use an uncorrelated data set and present formulas
for this scenario. These can easily be generalized when correlated
uncertainties are present., we can define the GMM likelihood as

$\displaystyle L(Y|\vec{\theta})=\prod_{j=1}^{N_{\text{pt}}}\pi(y_{j},\Delta
y_{j}|\vec{\theta})=\prod_{j=1}^{N_{\text{pt}}}\sum_{i=1}^{K}\omega_{i}\mathcal{N}(y_{j},\Delta
y_{j}|\theta_{i}),$ (5) $\displaystyle 0\leq\omega_{i}\leq
1\quad\text{and}\quad\sum_{i}\omega_{i}=1,$

where the $i$’th base distribution $\mathcal{N}(y_{j},\Delta
y_{j}|\theta_{i})$ is associated with a normalized weight $\omega_{i}$. Here
$\theta_{i}$ corresponds to the fit parameters. For example, for the fits
performed in the previous section
$\theta_{i}=\left\\{a_{4_{i}},a_{5_{i}}\right\\}$ where $i=1,2,...K$.
Therefore, performing a two parameter fit with K Gaussians, requires $2\times
K$ parameters. To simplify the following expressions, we will abbreviate
$\mathcal{N}$ and $\pi$ for the $i$’th base distribution and the $j$’th
measurement as follows

$\mathcal{N}(y_{j},\Delta
y_{j}|\theta_{i})=\mathcal{N}_{ij},\quad\pi(y_{j},\Delta
y_{j}|\vec{\theta})=\pi_{j}\ .$ (6)

Notice that the number of original parameters $N_{\text{parm}}$ has increased
in the GMM by a factor $K$. This means that we need to determine not only the
central value of $K\times N_{\text{parm}}$ parameters but also the uncertainty
on each of them. Later we will describe how this fit to $K\times
N_{\text{parm}}$ parameters can be combined in a meaningful way to produced a
fit and uncertainty on only $N_{\text{parm}}$ parameters.333 For example,
there were two fit parameters ($a_{4}$ and $a_{5}$) in the LS fits in Sec. II,
whereas in the GMM there would be $2\times K$ parameters. Ultimately, we want
to be able determine the uncertainty on two parameters and not $2\times K$
parameters. In the GMM, the logarithm of the likelihood function is

$\displaystyle-\ln L$ $\displaystyle=$
$\displaystyle-\ln\bigg{(}\prod_{j=1}^{N_{\text{pt}}}L(y_{j},\Delta
y_{j}|\vec{\theta})\bigg{)}=-\sum_{j=1}^{N_{\text{pt}}}\ln\bigg{(}\sum_{i=1}^{K}\omega_{i}\mathcal{N}_{ij}\bigg{)}$
(7) $\displaystyle=$
$\displaystyle-\sum_{j=1}^{N_{\text{pt}}}\ln\bigg{(}\sum_{i=1}^{K}\frac{\omega_{i}}{\sqrt{2\pi}\Delta
y_{j}}\text{exp}\Big{[}-\frac{1}{2}\Big{(}\frac{y_{j}-T(\theta_{i})}{\Delta
y_{j}}\Big{)}^{2}\Big{]}\bigg{)}.$

For a single Gaussian distribution with $K=1$, the log-likelihood in Eq. (7)
is proportional to $\chi^{2}$, which reproduces the Least Square method as
shown in the equation below

$-\ln L=-\sum_{j=1}^{N_{\text{pt}}}\ln\bigg{(}\mathcal{N}(y_{j},\Delta
y_{j}|\theta)\bigg{)}=\frac{1}{2}\chi^{2}+\sum_{j=1}^{N_{\text{pt}}}\ln\Big{(}\sqrt{2\pi}\Delta
y_{j}\Big{)}\ .$ (8)

Notice, in Eq. (7), we have chosen to define the likelihood, such that each
Gaussian has information about all experimental data points, i.e. for each
Gaussian we sum over all values of $j=1$ to $j=N_{\text{pt}}$. An alternate
way of defining the likelihood would be to distribute the experimental data
among the various Gaussians, i.e. each data set would have its own Gaussian.
These two definitions have different interpretations. In the latter, we are
trying to assign different theoretical predictions to different experiments,
whereas in the former we are trying to find the theoretical values that
different combinations of all experiments predict. Since we are trying to
combine all experiments into a single global fit, we focus on using the
definition as given in Eq. (7).

The GMM method allows us to extend the statistical model assumed in the
optimization procedure, from a single Gaussian distribution to a more
generalized configuration, where a multi-modal feature is allowed.
Importantly, for a consistent data set we will see that the multiple Gaussian
modes will overlap and yield a result that is equivalent to a single Gaussian,
i.e. the Least Squares fit. We will discuss the consistency between the GMM
and LS methods in Sec. V.

###  III.1 Performing parameter fits with the GMM

In order to determine the best fit parameters, we minimize the the log-
likelihood defined in Eq. (7). We remark here that, conventinally the
Expectation Maximization algorithm [31] is used to learn the best fit
parameters of the mixture model.444We remark here that it is possible to use
the Expectation Maximization algorithm if we produce a large pseudo-data set
and apply the Gaussian mixture model directly to the theory parameter space.
The advantage of this algorithm is that it provides an iterative method that
does not require information about the derivative of the likelihood. However,
in our setup, a straightforward implementation of the Expectation Maximization
algorithm is not possible since we are trying to fit a nonlinear function. For
the purpose of demonstration, we use methods that rely on information about
the derivative of the likelihood. Specifically, we use TensorFlow [32] to
determine derivatives of the likelihood and minimize the log-likelihood. As
described in the previous section, the number of parameters have increased and
for $K$ Gaussians we have

$\displaystyle N_{\text{parm}}^{\text{GMM}}=K\times N_{\text{parm}}+K-1\ ,$
(9)

where, the last two terms, $K-1$, come from the constraint on the weights
$\omega_{i}$ that are required to sum to unity. We carry out minimization of
the likelihood by varying all $N_{\text{parm}}^{\text{GMM}}$ parameters.555It
is possible to determine the value of the weights $\omega_{i}$ in other more
complicated ways since they indicate the degree of importance to each
Gaussian. Here we take the approach that we learn the values of $\omega_{i}$
directly from the data and hence it is also varied during the minimization
process. The gradient of the log-likelihood with respect to the fit parameters
that belong to $i$-th Gaussian is

$\displaystyle\frac{\partial}{\partial\theta_{i}}\ln L$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N_{\text{pt}}}\omega_{i}\frac{\mathcal{N}_{ij}}{\pi_{j}}\frac{y_{j}-T(\theta_{i})}{\Delta
y^{2}_{j}}\frac{\partial T(\theta_{i})}{\partial\theta_{i}},$ (10)

and that with respective to weight $\omega_{i}$ is

$\displaystyle\frac{\partial}{\partial\omega_{i}}\ln L$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N_{\text{pt}}}\frac{\mathcal{N}_{ij}}{\big{(}\sum_{i=1}^{K}\omega_{i}\mathcal{N}_{ij}\big{)}}=\sum_{j=1}^{N_{\text{pt}}}\frac{\mathcal{N}_{ij}}{\pi_{j}}\,.$
(11)

Optimal values, $\hat{\theta}_{i}$, for the parameters ${\theta}_{i}$ are
found by extremizing the log-likelihood by solving [33]

$\displaystyle\frac{\partial}{\partial\theta_{i}}\ln
L\bigg{|}_{\theta_{i}=\hat{\theta}_{i}}$ $\displaystyle=$ $\displaystyle 0$
(12)

along with a similar equation involving the weights $\omega_{i}$,

$\displaystyle\frac{\partial}{\partial w_{i}}\ln
L\bigg{|}_{w_{i}=\hat{w}_{i}}$ $\displaystyle=$ $\displaystyle 0.$ (13)

In this work, we extremize the log-likelihood function by using the Adam
optimizer with learning rate being $10^{-6}$ provided by the TensorFlow API
[32]. Note that, in Eq. (10), the best-fit value of the prediction
$y_{j}(\theta_{i})$ is weighted not only by the variance in measurements
($\Delta y_{i}$), but also by the weights $\omega_{i}$ and the factor
$\mathcal{N}(y_{j}|\theta_{i})/\pi(y_{j},\Delta y_{j}|\vec{\theta})$ of the
$i$-the base distribution, which becomes trivial for a single Gaussian. The
Jacobian $\partial y_{j}(\theta_{i})/\partial\theta_{i}$, transforms
information from the data space to the PDF space. The presence of this
transformation makes it difficult to use the Expectation Maximum algorithm.

In this work we use the above equation as a constraint that determines the
value of $\omega_{i}$. This is generally the approach taken for other
unsupervised machine learning applications. In general, this strategy tends to
assign higher weights to experiments with higher precision. If desired, it is
possible to control the weights by increasing or decreasing the value of the
uncertainties associated with a certain experimental data set.666One method of
trying to deal with inconsistent data sets is to artificially inflate
uncertainties until we arrive at a consistent data set [16]. Although not
described in this work, we remark that the GMM in conjunction with certain
information criteria, provides a way to determine the extent to which it is
necessary to inflate the uncertainty on data until the fits become consistent.
Thus along with $\theta_{i}$, the weights $\omega_{i}$ are also varied in
order to optimize the log-likelihood.

###  III.2 Mean and Uncertainty in the GMM

It is straightforward to show that the expectation value of parameters is the
weighted sum of the optimal parameters for individual base distributions,

$\mathbb{E}[\theta]=\sum_{i=1}^{K}\omega_{i}\hat{\theta}_{i}.$ (14)

Consequently, an experiment with higher precision would be associated with a
higher weight and, as a result, the central value of the parameters is pulled
towards those values that better agree with higher precision experiments.

There are several ways to determine the uncertainty with the modified log-
likelihood. Below, we describe the use of the observed Fisher Information
Matrix to estimate uncertainty [33] on the fit parameters. We note that this
is only one way to estimate the uncertainty and there are other possible
methods which will be discussed in detail elsewhere.

For the case of single Gaussian or when all the Gaussians have the same mean
and covariance, the observed Fisher information matrix reduces to the usual
Hessian matrix used to determine the covariance of the fit parameters in the
LS method.777 For the GMM with multiple Gaussians, the likelihood function as
defined in Eq. (7) includes all experimental points for each Gaussian. Thus,
if the means of each Gaussian are identical, then this implies that the
covariance, which is related to the second order derivative about the
minimimum of the likelihood, must also be identical, c.f. Eq. (22). In
general, the covariance matrices for Gaussians with identical means may be
different in the GMM if the likelihoods were to be defined differently.  Like
the Hessian matrix, the information matrix can also be found with numerical
methods. Here we use algorithmic differentiation [34] within TensorFlow to
determine derivatives of the likelihood function. We note that the information
matrix is only an estimate of the covariance of the GMM and it is possible to
determine the covariance matrix as well, which will be described in detail
elsewhere.

Elements of the information matrix ($\mathcal{I}$) are determined by taking a
further derivative on Eq. (12) as follows,

$\displaystyle\left(\mathcal{I}\right)_{mk}$
$\displaystyle=-\frac{\partial^{2}}{\partial\theta_{m}\partial\theta_{k}}\ln
L=\frac{\partial}{\partial\theta_{m}}\left(-\sum_{j=1}^{N_{\text{pt}}}\frac{\omega_{k}}{\pi_{j}}\frac{\partial\mathcal{N}_{kj}}{\partial\theta_{k}}\right)$
$\displaystyle=\sum_{j=1}^{N_{\text{pt}}}\left[-\frac{\omega_{k}}{\pi_{j}}\frac{\partial^{2}\mathcal{N}_{mj}}{\partial\theta_{m}\partial\theta_{m}}\delta_{mk}+\frac{\omega_{k}\omega_{m}}{\pi_{j}^{2}}\frac{\partial\mathcal{N}_{mj}}{\partial\theta_{m}}\frac{\partial\mathcal{N}_{kj}}{\partial\theta_{k}}\right]$
(15)

Notice that for a fit with $N_{\text{parm}}$ parameters per Gaussian, the
information matrix is a $(KN_{\text{parm}})\times(KN_{\text{parm}})$ symmetric
matrix, i.e. it spans a space that is a direct product of an $N_{\text{parm}}$
parameter and a $K$ parameter space. The second term contributes to the
elements of the information matrix ($\mathcal{I}$) that are off-diagonal in
$K$-space but not in $N_{\text{parm}}$ space. In contrast, the off-diagonal
terms in $N_{\text{parm}}$ space contribute to the covariance matrix. On the
other hand, the $K$-space off-diagonal terms are small since they must add up
to the left hand side of Eq. (12), which vanishes at
$\theta_{i}=\hat{\theta}_{i}$ by definition. This is also verified numerically
for all of the fits performed here. For simplicity, we will neglect these
particular $K$-space off-diagonal terms and only consider the
$N_{\text{parm}}\times N_{\text{parm}}$ block diagonal terms of the
information matrix. Although, it is not necessary to neglect these off
diagonal terms, this allows us to quickly approximate the probability
distribution function in parameter space ($\theta$) as a mixture of Gaussians.
In other words, we are trying to use the Likelihood function defined in Eq.
(7), to approximate the probability distribution function of $\theta$. Thus,
we assume that, when evaluated at the mean ($\theta_{i}=\hat{\theta}_{i}$) the
information matrix takes on an approximate block diagonal form, with each
$N_{\text{parm}}\times N_{\text{parm}}$ block belonging to the subspace of the
$i$’th Gaussian $\mathcal{N}(y_{j},\Delta y_{j}|\theta_{i})$. We denote the
information matrix for individual base distributions of the GMM as
$\mathcal{I}_{\text{GMM},i}=-(\partial^{2}\ln
L/\partial\theta_{i}\partial\theta_{i})/\omega_{i}$, which are the unweighted
$N_{\text{parm}}\times N_{\text{parm}}$ diagonal blocks of the full
information matrix $\mathcal{I}$. An estimate of the covariance of individual
base Gaussian distributions of the GMM is then,

$\text{cov}_{\text{GMM},i}\simeq\bigg{(}\mathcal{I}_{\text{GMM},i}\bigg{)}^{-1}=\Bigg{(}\sum_{j=1}^{N_{\text{pt}}}\left[-\frac{1}{\pi_{j}}\frac{\partial^{2}\mathcal{N}_{ij}}{\partial\theta_{i}^{2}}+\frac{\omega_{i}}{\pi_{j}^{2}}\left(\frac{\partial\mathcal{N}_{ij}}{\partial\theta_{i}}\right)^{2}\right]\Bigg{)}^{-1}$
(16)

We can now define the probability distribution of the parameters of the fit as
mixture of Gaussians as follows:

$\displaystyle\mathcal{P}(x|\theta_{i},\text{cov}_{\text{GMM,i}})$
$\displaystyle=\sum_{i=1}^{K}\frac{\omega_{i}}{(2\pi)^{N_{\text{parm}}/2}\det(\text{cov}_{\text{GMM,i}})^{1/2}}$
$\displaystyle\times\text{exp}\Big{[}-\frac{1}{2}\Big{(}(x-\theta_{i})^{T}\cdot\left(\text{cov}_{\text{GMM,i}}\right)^{-1}\cdot(x-\theta_{i})\Big{)}\Big{]}\
.$ (17)

With this definition, we can now determine estimates of the uncertainty on the
parameters $\theta$. We first note that when $K=1$, we have a single Gaussian
and Eq. (15) reduces to the Hessian matrix of the LS method,

$-\frac{\partial^{2}}{\partial\theta\partial\theta}\ln
L\bigg{|}_{\theta=\hat{\theta}}=\frac{1}{2}\frac{\partial^{2}}{\partial\theta\partial\theta}\chi^{2}\bigg{|}_{\theta=\hat{\theta}}=\sum_{j=1}^{N_{\text{pt}}}\frac{1}{\Delta
y^{2}_{j}}\frac{\partial^{2}T(\theta)}{\partial\theta^{2}}\bigg{|}_{\theta=\hat{\theta}},$
(18)

and the covariance is

$\text{cov}_{\chi^{2}}=\bigg{(}\sum_{j=1}^{N_{\text{pt}}}\frac{1}{\Delta
y^{2}_{j}}\frac{\partial^{2}T(\theta)}{\partial\theta^{2}}\bigg{|}_{\theta=\hat{\theta}}\bigg{)}^{-1}.$
(19)

For the GMM, in order to provide a useful estimate of uncertainty, we need to
reduce the $(KN_{\text{parm}})\times(KN_{\text{parm}})$ space to a
$(N_{\text{parm}})\times(N_{\text{parm}})$ space in a meaningful way. To that
end, let us consider the covariance of the GMM which is defined as

$\displaystyle\text{cov}_{\text{GMM}}$ $\displaystyle=$
$\displaystyle\mathbb{E}[\theta^{2}]-\mathbb{E}[\theta]^{2}$ (20)
$\displaystyle=$
$\displaystyle\sum_{i=1}^{K}\omega_{i}\hat{\theta}_{i}^{(2)}-\big{(}\sum_{i=1}^{K}\omega_{i}\hat{\theta}_{i}\big{)}^{2},$

where $\hat{\theta}_{i}^{(2)}$ denotes the estimated second moment of
parameters for the $i$-th base distribution. It can be re-expressed using the
covariance of $i$-th base Gaussian distribution,

$\hat{\theta}_{i}^{(2)}=\text{cov}_{\text{GMM},i}+\hat{\theta}^{2}_{i}.$ (21)

Combining Eqs. (20) and (21), an estimate of the covariance of the GMM is
given by

$\displaystyle\text{cov}_{\text{GMM}}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{K}\omega_{i}\
\text{cov}_{\text{GMM},i}+\sum_{i=1}^{K}\omega_{i}(\mathbb{E}[\theta]-\hat{\theta}_{i})^{2}$
(22) $\displaystyle\simeq$
$\displaystyle\sum_{i=1}^{K}\omega_{i}\Bigg{(}\sum_{j=1}^{N_{\text{pt}}}\left[-\frac{1}{\pi_{j}}\frac{\partial^{2}\mathcal{N}_{ij}}{\partial\theta_{i}^{2}}+\frac{\omega_{i}}{\pi_{j}^{2}}\left(\frac{\partial\mathcal{N}_{ij}}{\partial\theta_{i}}\right)^{2}\right]\Bigg{)}^{-1}+\sum_{i=1}^{K}\omega_{i}(\mathbb{E}[\theta]-\hat{\theta}_{i})^{2}.$

The covariance of the GMM shown in the above equation receives two
contributions: the first is the weighted sum of covariances for the individual
base distributions, and the second term accounts for differences in the mean
values of each of the base distributions. When the means of each Gaussian
overlap, the second term vanishes and the first term reproduces the covariance
in LS method, Eq. (19).

##  IV Application of the GMM in a toy model of PDFs

In this section, we compare results obtained by the GMM method, defined in Eq.
(7), to those obtained using the LS method. When estimating uncertainty using
the LS method, we first diagonalize the Hessian and determine eigenvalue
directions. We then determine values of the parameters that satisfy the
$\Delta\chi^{2}=1$ condition and define the uncertainty on the fit parameters,
which corresponds to the 68% confidence level interval. Here, we also
implement the CT [1] and MSHT [2] tolerance criteria, which provide a more
conservative estimate of uncertainty under the potential imperfections of
fits. The comparison is presented in three cases of pseudo-data, including the
two pseudo-data sets introduced in Sec. II.1.

###  IV.1 Applying GMM to different scenarios

In this section we consider three cases. In each case we combine certain
pseudo-data sets. We list and describe each of these cases below:

  * •

Case-1: Here we use pseudo-data sets $\\#~{}1$ and $\\#~{}2$ which were
introduced in Sec. II.1. This case represents the scenario when there is
strong tension between data sets. With this case, we will demonstrate how the
GMM can provide an estimate of uncertainty in a statistically robust way.

  * •

Case-2: Here we combine pseudo-data sets $\\#~{}3$ and $\\#~{}4$. The central
values used to generate these data sets are the same as those used in
$\\#~{}1$ and $\\#~{}2$, however, their uncertainties are tripled by
increasing the value of $\alpha$ in Eq. (2) by a factor of three, as shown in
Table 3. The central values of data sets $\\#~{}1$ and $\\#~{}3$ as well as
$\\#~{}2$ and $\\#~{}4$ do not overlap, because different random numbers
($r$), cf. Eq. (3), are used to generate them. With this combination we will
be able to show how the GMM behaves when the errors of the data sets are large
enough so that there is no tension. We can connect this to the alternative
suggestions of dealing with uncertain data by artificially increasing the
errors of the data sets [16].

  * •

Case-3: This case contains pseudo-data sets $\\#~{}5$ and $\\#~{}6$. The
tension in this case is generated only via variation of the parameter $a_{2}$
. The parameter $a_{2}$ controls the behaviour of the pseudo-PDF in the
large-$x$ region, where the uncertainty is large. We fit this case with
parameters $a_{2}$ and $a_{4}$. With this combination we demonstrate
overfitting with the GMM method, and later in Sec. V.2 we will demonstrate how
to deal with overfitting.

The true values of the parameters for cases 1 through 3 are presented in Table
3 and Fig. 3.

case | pseudo-data | $\alpha$ | $N_{\text{pt}}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$
---|---|---|---|---|---|---|---|---|---
Case-1 | #1 | 0.05 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.4 | -3.0
| #2 | 0.05 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.6 | -2.8
Case-2 | #3 | 0.15 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.4 | -3.0
| #4 | 0.15 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.6 | -2.8
Case-3 | #5 | 0.05 | 50 | 30 | 0.5 | 2.4 | 4.3 | 2.4 | -3.0
| #6 | 0.05 | 50 | 30 | 0.5 | 1.4 | 4.3 | 2.4 | -3.0
Table 3:  Similar to Table 1, the summary of true values of the pseudo-PDF
parameters and values of $\alpha$ in Eq. (2) for three cases. For comparison,
we also include the true values of the pseudo-PDF parameters for Case-1 in
this table.

(a) Case-2

(b) Case-2

(c) Case-3

(d) Case-3

Figure 3:  Similar to Fig. 1, but for Case-2 and Case-3.  Case-1 | fits | pseudo-data | $a_{4}$ | $a_{5}$ | $\chi^{2}_{\\#1}/N_{\text{pt}}$ | $\chi^{2}_{\\#2}/N_{\text{pt}}$ | $\chi^{2}_{\\#1+\\#2}/N_{\text{pt}}$ | weight
---|---|---|---|---|---|---|---|---
| SG-A | # 1 | 2.32 | -3.22 | 0.88 | 6.55 | 3.72 | -
| SG-B | # 2 | 2.63 | -2.73 | 7.00 | 1.02 | 4.01 | -
| SG-C | # 1 and # 2 | 2.48 | -2.94 | 2.27 | 2.56 | 2.42 | -
| GMM Gauss.-1 | # 1 and # 2 | 2.31 | -3.25 | 0.88 | 6.55 | 3.72 | 0.52
| GMM Gauss.-2 | # 1 and # 2 | 2.64 | -2.72 | 7.00 | 1.02 | 4.01 | 0.48
Case-2 | fits | pseudo-data | $a_{4}$ | $a_{5}$ | $\chi^{2}_{\\#3}/N_{\text{pt}}$ | $\chi^{2}_{\\#4}/N_{\text{pt}}$ | $\chi^{2}_{\\#3+\\#4}/N_{\text{pt}}$ | weight
| SG-A | # 3 | 2.20 | -3.50 | 0.76 | 1.68 | 1.22 | -
| SG-B | # 4 | 2.74 | -2.54 | 1.61 | 0.89 | 1.25 | -
| SG-C | # 3 and # 4 | 2.48 | -2.98 | 0.96 | 1.11 | 1.04 | -
| GMM Gauss.-1 | # 3 and # 4 | 2.40 | -3.10 | 0.81 | 1.36 | 1.09 | 0.69
| GMM Gauss.-2 | # 3 and # 4 | 2.66 | -2.73 | 1.68 | 0.90 | 1.29 | 0.31
Case-3 | fits | pseudo-data | $a_{2}$ | $a_{4}$ | $\chi^{2}_{\\#5}/N_{\text{pt}}$ | $\chi^{2}_{\\#6}/N_{\text{pt}}$ | $\chi^{2}_{\\#5+\\#6}/N_{\text{pt}}$ | weight
| SG-A | # 5 | 2.27 | 2.40 | 1.16 | 2.33 | 1.75 | -
| SG-B | # 6 | 1.25 | 2.41 | 2.46 | 0.82 | 1.64 | -
| SG-C | # 5 and # 6 | 1.60 | 2.41 | 1.54 | 1.14 | 1.34 | -
| GMM Gauss.-1 | # 5 and # 6 | 2.06 | 2.40 | 1.19 | 1.96 | 1.58 | 0.58
| GMM Gauss.-2 | # 5 and # 6 | 1.13 | 2.41 | 2.58 | 0.82 | 1.70 | 0.42
Table 4:  The best-fit values of parameters, and $\chi^{2}/N_{\text{pt}}$
values for individual and combined pseudo-data sets obtained by the LS method
and the GMM fit.

The best-fit obtained by the GMM method for the Case-1 is shown in Figs. 4 and
4 and Table 4. For the GMM likelihood function, we choose to fit with a
mixture of two Gaussian distributions. As a result, the mean values of each
Gaussian, can describe the two distinct pseudo-data sets well and in fact are
very similar to what we obtain for SG-A and SG-B. We find the weights
$\omega_{i}$ to be roughly equal, which is a result of the errors in each of
the data sets $\\#1$ and $\\#2$ being almost equal and also having the same
number of data points. In Figs. 4and 4, the pseudo-PDF uncertainty estimated
by the GMM is compared to uncertainties predicted using the LS method. For
$10^{-3}<x<0.1$, the (small) error band of SG-C pseudo-PDF exclude those of
SG-A and SG-B fits, the error band estimated by the GMM is consistent and
includes areas covered by both SG-A and SG-B pseudo-PDFs. For $x>0.1$, where
the error in the data sets are large, uncertainties from SG-C and the GMM
overlap. Fig. 4 shows the results of the different fits in the fit parameter
space. The 1-$\sigma$ correlation ellipse estimated by SG-C is small, as has
been mentioned before, and does not overlap with SG-A and SG-B replicas. On
the other hand, the GMM correlation ellipse is large enough and, importantly,
also possesses a different axes orientations, that enables the GMM covariance
to overlap with some of the replicas of both SG-A and SG-B. Note, that we have
represented the GMM by a single covariance matrix but this should not be
confused with the covariance arising from a single Gaussian, since the shape
of the GMM likelihood is multi-modal and not uni-modal.

For Case-2, shown in Fig. 5, the tension in the data set has been reduced by
increasing the error on the data, so that SG-C describes the combined pseudo-
data set with $\chi^{2}_{\\#3+\\#4}/N_{\text{pt}}=1.04$. Even for such a mild
tension, the fits SG-A and SG-B do not overlap with SG-C. On the other hand,
the GMM fit overlaps with both SG-A and SG-B. The difference between the axes
orientations of the correlation ellipses between SG-C and the GMM is small.
Although, the GMM method produces a larger uncertainty which is more
consistent with SG-A and SG-B, the differences with SG-C are small and we need
a way to determine whether we should stick with SG-C or with the GMM method.
We will return to a discussion of this in Sec. V.2.

The performances for Case-3 of the GMM and LS methods are similar to Case-1.
The LS method cannot fit both pseudo-data sets # 5 and # 6 well, c.f. the
bottom rows of Table 4, even though the uncertainties in both data sets are
large. The GMM is able to learn about the tension in these data sets and shows
a strong inclination to two modes. In terms of the fitted pseudo-PDFs, the GMM
in Fig. 6 systematically enlarges the uncertainty and envelopes both data sets
# 5 and # 6 simultaneously.

Figure 4:  Pseudo-PDFs obtained by the GMM method and the LS method for Case-1
are compared. Panel (a): The pseudo-data sets are compared to best-fit pseudo-
PDFs. Although shown separately, the two best-fit pseudo-PDFs from the GMM
method are eventually combined according to weights shown in Table 4. Panel
(b): Similar to panel (a), but presented as ratios of the best-fit obtained
from the SG-C fit. Panel (c): The absolute values of pseudo-PDFs for all SG-A,
-$B$, -$C$, and GMM fits. Panel (d): Similar to panel (c), but presented as
ratios to the best-fit obtained from the SG-C Hessian set. Panel (e): The
correlation ellipse from the GMM method is compared to SG-C correlation
ellipses, and to replicas from SG-A, -$B$, and -$C$ fits. The best-fit points
are located at their corresponding values in Table 4, while the GMM best-fit
corresponds to the weighted average of best-fit of two individual Gaussian
distributions as given by Eq. (14).

Figure 5:  Similar to Fig. 4, but for the Case-2.

Figure 6:  Similar to Fig. 4, but for the Case-3.

###  IV.2 The GMM method compared to different tolerance criteria

Estimations of PDF uncertainties are needed for inference from QCD experiments
at the LHC and other facilities. Various factors that determine the accuracy
of PDFs include the theoretical uncertainty, the insufficiently understood
source of uncertainty due to the prescription for experimental systematic
errors, the tension between experimental data, and the sampling of PDF
functional forms and fitting methodologies, as detailed in Ref. [35]. Global
PDF analysis groups are aware of the shortcomings of using a fixed tolerance
$\Delta\chi^{2}=1$ and have proposed alternate methods. Specifically, CTEQ-TEA
[1] and MSHT [2] global PDF analyses, both adopt certain tolerance criteria,
known as the two-tier tolerance and the dynamic tolerance criteria,
respectively, on top of the basic Hessian method.

In Sec. IV.1, we compared results of the GMM with the usual LS method of
determining 1-$\sigma$ contours by setting $\Delta\chi^{2}=1$. In this
subsection, we compare the GMM results presented in Sec. IV.1 to the Hessian
results with the CT and MSHT tolerance criteria.

We remind the reader that CTEQ-TEA Hessian eigenvector PDFs are constructed
using a two-tier prescription [36, 37] that prevents unexpectedly large
increments to either the global $\chi^{2}$ or to the $\chi^{2}_{E}$ values of
individual experiments. Hence, the CT PDF error sets are determined using a
two-tier (global+dynamic) tolerance ($T$) that accounts for theory,
experimental, parametrization, and methodological uncertainties and results in
error bands that are wider than the $\Delta\chi^{2}=1$ criterion. The CT
Hessian uncertainty, at the 68% confidence level, nominally corresponds to
$T^{2}\equiv\Delta\chi^{2}_{total}\approx 37$, where $\chi^{2}_{total}$ is the
sum of the global $\chi^{2}$ and the so-called Tier-2 (dynamic) penalty [36,
37]. The MSHT PDF fits by default apply the “dynamic tolerance” procedure, as
detailed in [38, 39, 2]. This enlarges the uncertainties beyond the
$\Delta\chi^{2}=1$ definition to account for data sets with tension, as well
as potential mismatch of data and theory due to imprecision in the theory or
parametrization and experimental measurements. The impact of different
tolerance criteria on the determination of PDF induced uncertainties has been
examined in Ref. [40]. In this study, we compare the results of applying the
CT and MSHT tolerance criteria with $T^{2}=37$ to the ones with
$\Delta\chi^{2}=1$.

The Hessian uncertainties based on the single Gaussian LS method with a
tolerance criteria of $T^{2}=37$ for all three cases are compared to the GMM
method in Figs. 7 to 9.

  * •

Case-1: Both Hessian sets with CT and MSHT tolerance criteria estimate a
smaller pseudo-PDF uncertainty for $10^{-3}<x<0.1$, compared to the GMM.
Recall for Case-1, the GMM produced a more faithful representation of
uncertainties and overlapped with SG-A and SG-B. We see that the alternate
tolerance criteria is not able to account for the larger uncertainties in this
$x$ range. Further, the uncertainties given by the modified tolerance criteria
are much larger than required when $x>0.1$. The same can also be observed in
the fit parameter space. Here, the SG-C correlation ellipses with tolerance
cover larger regions than the default Hessian set, but still cannot predict a
correlation that spans regions of parameter space that cover portions of SG-A
and SG-B.

  * •

Case-2: As shown in Fig. 8, the larger tolerance leads to an overestimation of
pseudo-PDF uncertainty, particularly for $x>0.1$. In fact, the correlation
ellipse for the larger tolerance, seen in the fit parameter space, has regions
that are not populated by any of the replica sets.

  * •

Case-3: As shown in Fig. 9(b), the GMM uncertainty (gray shaded regions of the
pseudo-PDF ratios) shows better agreement with the uncertainty bands of the
higher tolerance criteria. However, when comparing results in the fit
parameter space, shown in Figs. 9(a), there appears to be significant mismatch
between the higher tolerance ellipse and the GMM ellipse.

The tolerance criteria, currently in use by PDF collaborations, was introduced
to capture larger uncertainties when there is tension in the fits. However,
these methods are still based on the assumption that the likelihood behaves
like a single Gaussian and is uni-modal. The tolerance criteria are not able
to faithfully capture the effect of tension in the data. In some regions, the
higher tolerance underestimates the uncertainty, and in others it
overestimates the uncertainty. The GMM, on the other hand, is able to provide
a more faithful representation of the combined likelihoods of two data sets
with tension.

(a) GMM compared to tolerance criteria

(b) pseudo-PDF ratio

Figure 7:  The pseudo-PDF uncertainty estimated by the GMM method in Case-1 is
compared to pseudo-PDF uncertainties predicted by the LS method along with the
CT and MSHT tolerance criteria.

(a) GMM compared to tolerance criteria

(b) pseudo-PDF ratio

Figure 8:  Similar to Fig. 7, but for the Case-2.

(a) GMM compared to tolerance criteria

(b) pseudo-PDF ratio

Figure 9:  Similar to Fig. 7, but for the Case-3.

##  V Consistency checks for the GMM

In Sec. IV, we demonstrated how the GMM with two base Gaussian distributions
is able to successfully describe the combined likelihood of two pseudo-data
sets with tension. In this section, we show how the GMM reduces to the usual
single Gaussian LS method when data are consistent. Further, as with any
present day machine learning technique, when implementing the GMM, one must
deal with over-fitting. In this section, we also remark on the consequences of
over-fitting with the GMM.

###  V.1 Consistent results from consistent data

In Sec. III, we have shown how the GMM likelihood function reduces to the LS
method when the number of Gaussians used in the GMM is $K=1$. Another
possibility for reducing to the usual $\chi^{2}$ fit is when all base
Gaussians end up with the same mean and covariance, i.e.

$\mathbb{E}[\theta]\sim\hat{\theta}_{i}\sim\hat{\theta}_{k},\quad\text{cov}_{\text{GMM},i}\sim\text{cov}_{\text{GMM},k},\quad\text{and}\quad\mathcal{N}(y_{j},\Delta
y_{j}|\theta_{i})\sim\mathcal{N}(y_{j},\Delta y_{j}|\theta_{k}).$ (23)

Consequently, the second term in Eq. (22), which accounts for the difference
of means, vanishes. In this scenario, the covariance of the GMM reduces to

$\displaystyle\text{cov}_{\text{GMM}}$ $\displaystyle\simeq$
$\displaystyle\sum_{i=1}^{K}\omega_{i}\Bigg{(}\sum_{j=1}^{N_{\text{pt}}}\left[-\frac{1}{\pi_{j}}\frac{\partial^{2}\mathcal{N}_{ij}}{\partial\theta_{i}^{2}}+\frac{\omega_{i}}{\pi_{j}^{2}}\left(\frac{\partial\mathcal{N}_{ij}}{\partial\theta_{i}}\right)^{2}\right]\Bigg{)}^{-1}+\sum_{i=1}^{K}\omega_{i}(\mathbb{E}[\theta]-\hat{\theta}_{i})^{2}$
(24) $\displaystyle\simeq$
$\displaystyle\bigg{(}\sum_{j=1}^{N_{\text{pt}}}\frac{1}{\Delta
y^{2}_{j}}\bigg{)}^{-1},$

which is the covariance found by the LS method in Eq. (19). When fitting a
consistent data set, consistency between the GMM method and LS method are
expected as shown analytically in Eq. (24). In practice, due to overfitting,
we will see that the GMM reduces to the LS method only if no local minima
exist in the log-likelihood, i.e. there is only one minimum of the log-
likelihood.

Here, we verify Eq. (24) by performing two more fits on pseudo-data sets, as
shown in Fig. 10:

  * •

Case-4: This case concerns a fit to the pseudo-data $\\#~{}7$, which has the
same truth as pseudo-data $\\#~{}1$, but without any statistical fluctuations.
In this case, only a single extremum of the log-likelihood of the GMM exists,
i.e. $\theta_{i}=\theta_{k}$, and we shall see that numerically as well, the
GMM will reduce to the LS method fit.

  * •

Case-5: This case concerns a fit to the pseudo-data $\\#~{}1$. Here, due to
overfitting, we will see that although $\theta_{i}\simeq\theta_{k}$, the GMM
tends to find multiple possible values of $\theta_{k}$. This in turn implies
that the results from the GMM method are different from the LS method.

(a) Truth for Case-4

(b) Truth for Case-5

(c) best-fits for Case-4

(d) best-fits for Case-5

Figure 10:  Upper panels: the truths are compared to pseudo-data for Case-4
and Case-5. Lower panels: the best-fit pseudo-PDFs are compared to pseudo-data
for both cases.

We begin by looking at results from the fit in Case-4. The uncertainties
predicted by both methods are shown in Figs. 11(a), 11(c), 12(a), and 12(c).
The figures indicate, as expected, that the central value of the fit for both
the GMM and the LS methods are in good agreement. In Fig. 11(b) and Fig.
11(c), we see that the ellipses from the covariance matrix estimated using
both the GMM and the LS methods overlap and are identical. On the other hand,
for Case-5, we notice in Fig. 11(b) and Fig. 11(d), that the covariance
predicted by the GMM method $\text{cov}_{\text{GMM}}$ is marginally larger
than the covariance given by the LS method $\text{cov}_{\chi^{2}}$. Fig. 11(d)
shows that the covariance of individual base Gaussian distributions of the GMM
$\text{cov}_{\text{GMM},i}$, which is larger than that obtained by the LS
method $\text{cov}_{\chi^{2}}$, and that the optimal values of parameters
$\hat{\theta}_{i}$ for both base Gaussian distributions of the GMM differ by a
small amount. The overestimation of uncertainty is a result of overfitting in
GMM, where the GMM ends up learning statistical fluctuations rather than real
tension between data sets.

The GMM method looks for tension and tries to categorize data into two or more
different Gaussians. In general, for the GMM, overfitting gives rise to larger
uncertainties and in that sense the GMM always provides a conservative
estimate of the uncertainty.888This is in contrast to overfitting by
increasing the number of fit parameters in a LS fit, where overfitting may
lead to a reduction in $\chi^{2}/N_{\text{pt}}<1$. In this sense, the GMM
method provides conservative estimates of the uncertainty. However, one would
like to provide a precise PDF which is essential for the physics goals of the
community. In order to determine if there is overfitting in the GMM, we
describe certain information criteria in the next section that can be used to
determine the optimal number of Gaussians, i.e., the value of $K$. The result
of Case-5 suggest that applying the GMM to a consistent data set will lead to
larger uncertainties, which although conservative, should also be avoided.

(a) Correlations for Case-4

(b) Correlations for Case-5

(c) Correlations for Case-4

(d) Correlations for Case-5

Figure 11:  Upper panels: the comparison of covariances obtained by the GMM
method and the LS method, for the Case-4 and Case-5. Lower panels: similar to
upper panels, but include the covariances of individual base Gaussian
distributions of the GMM.

(a) Pseudo-PDF ratio for Case-4

(b) Pseudo-PDF ratio for Case-5

(c) Size of pseudo-PDF uncertainty for Case-4

(d) Size of pseudo-PDF uncertainty for Case-5

Figure 12:  Upper panels: the comparison of pseudo-PDF ratio obtained by both
methods for Case-4 and Case-5, given the reference being the central value
obtained by the LS method. Lower panels: the comparison of the sizes of
pseudo-PDF uncertainty obtained by the GMM method and the LS method, for
Case-4 and Case-5.

###  V.2 Information Criteria for an Optimal number of base Gaussians

The results shown in Sec. V.1 indicate that the GMM method might over-estimate
the uncertainty, especially if one includes more base Gaussian distributions
than needed. Therefore, a procedure to determine the optimal value of $K$, for
a mixture of $K$ base Gaussian distributions of the GMM, is required. Here, we
propose the use of the Akaike information criterion (AIC) [41] and Bayesian
information criterion (BIC) [42] as indicators that will help determine
optimal values for $K$, i.e., the number of Gaussians. The values of AIC and
BIC, defined in Eq. (25), increase as more parameters are introduced. They
also take into account the change in goodness of fit and reduce as the fit
improves. Both AIC and BIC are thus able to account for the trade-off between
introducing additional parameters and an improvement to the fit. The minimum
values for the AIC and BIC indicate to us the best choice of parameters, in
this case, the best choice for the value of $K$. The AIC and BIC are defined
as

AIC $\displaystyle=$ $\displaystyle N_{p}\ \ln N_{\text{pt}}-2\ln
L\big{|}_{\theta=\hat{\theta}},$ BIC $\displaystyle=$ $\displaystyle
2N_{p}-2\ln L\big{|}_{\theta=\hat{\theta}}.$ (25)

We count the number of independent parameters $N_{p}$ for the GMM as
following,

$N_{p}=K\times N_{\text{parm}}+(K-1),$ (26)

where $N_{\text{parm}}$ is the number of independent parameters for each base
Gaussian distribution. The first term in Eq. (26) is the total number of
fitting parameters over $K$ copies of base distributions. The second term
corresponds to the number of independent weights. In our pseudo-data study,
there are always two free variables in the pseudo-PDF parametrization. Hence
the total number of independent parameters for the GMM method in our current
fits is

$N_{\text{parm}}=2K+(K-1).$ (27)

We scan the AIC and BIC of the GMM up to $K=4$ for all 5 cases considered in
this work, as summarized in Table 5. For Case-1 and Case-3, both AIC and BIC
suggest two base Gaussian distributions would be optimal in order to
faithfully extract information from their pseudo-data sets. This confirms the
results shown in Sec. IV, particularly in Figs. 4 and 6, where two base
Gaussian distributions can estimate a correlation that covers the two pseudo-
data sets in tension. The AIC and BIC suggest for Case-2 a sinlge Gaussian
distribution is enough to fit this pseudo-data set. In Fig. 5, the LS method
predicts a correlation, whose slope is close to the correlation estimated with
two base Gaussian distributions. This is consistent with the fact that Case-2
has minimal tension due to the large values of data errors. At this point we
would like to point out that the GMM in conjunction with AIC and BIC, can be
used to determine how much the errors in the data should be inflated in order
to make them consistent so that we can use the LS method, as suggested in Ref.
[16, 27, 26, 28, 29].999 Specifically, one possible measure is to increase the
errors in the data sets until the AIC and BIC indicate that $K=1$ is the
optimal value of $K$. This provides an alternate way to determine the scale
factor used in Ref. [26]. For Case-4 and Case-5, the AIC and BIC of course
suggest that there is no need of using the GMM method, since by doing so one
is over-fitting the data. We reiterate that the consequence of over-fitting is
to increase our uncertainty and therefore leads to a conservative estimate of
uncertainty, i.e., using the optimum value of $K$ provides the most precise
PDFs, whereas larger values of $K$ will lead to PDFs with larger uncertainty.
We also notice in Table 5 that a higher number of base distributions $K$
sometimes would lead to equal or slightly lower value of $-\ln L$ than the $K$
suggested by AIC and BIC, indicating that one can equally or even better fit
the data with more base distributions. However this is likely due to over-
fitting the statistical fluctuations in the data. To conclude, although the
GMM method is expected to reduce to the LS method for a consistent data set,
one should avoid using it without detailed assessment of the level of tension.
The combination of AIC and BIC indicators with the generalized statistical
model of the GMM method could effectively reproduce the correct correlation in
the presence of inconsistency in data. At the same time we can determine
optimal values of $K$ and provide a formalism to determine more precise PDFs.

|  | $K=1$ | $K=2$ | $K=3$ | $K=4$
---|---|---|---|---|---
Case-1 | AIC | -102.2 | -203.6 | -194.9 | -187.9
| BIC | -106.1 | -211.2 | -206.4 | -203.2
$N_{\text{pt}}$=100 | $-\ln L$ | -55.0 | -109.6 | -109.2 | -109.6
Case-2 | AIC | -21.2 | -15.4 | -7.9 | -0.2
| BIC | -25.0 | -23.0 | -19.3 | -15.5
$N_{\text{pt}}$=100 | $-\ln L$ | -14.5 | -15.5 | -15.7 | -15.7
Case-3 | AIC | -219.3 | -220.2 | -212.8 | -205.0
| BIC | -223.2 | -227.8 | -224.3 | -220.3
$N_{\text{pt}}$=100 | $-\ln L$ | -113.6 | -117.9 | -117.9 | -118.1
Case-4 | AIC | -117.8 | -109.9 | -102.1 | -94.3
| BIC | -121.6 | -117.6 | -113.6 | -109.6
$N_{\text{pt}}$=50 | $-\ln L$ | -62.8 | -62.8 | -62.8 | -62.8
Case-5 | AIC | -169.3 | -161.5 | -153.6 | -145.8
| BIC | -173.1 | -169.1 | -165.1 | -161.1
$N_{\text{pt}}$=50 | $-\ln L$ | -88.6 | -88.6 | -88.6 | -88.6
Table 5:  The AIC, BIC and $-\ln L$ for Case-1 to Case-5 obtained by the GMM
with various numbers of base distributions. The lowest values of AIC, BIC and
$-\ln L$ for each cases are written in bold font.

##  VI Conclusion

An accurate PDF uncertainty estimation is a theoretical necessity for
precision measurements as well as new physics searches at the LHC and future
colliders. The sources of PDF uncertainty include the theoretical uncertainty,
the experimental systematic errors, and the sampling of PDF functional forms
and fitting methodologies, etc. [35]. The existence of tension among different
data sets can also reduce the precision of PDFs and necessitates the use of
alternate statistical models that provide more accurate estimates of
uncertainty. In this work, we described a formalism that generalizes the
commonly used Least-Squares (LS) method and provides a way to combine
likelihoods from different experiments in a statistically robust and
meaningful way.

The proposed generalized likelihood modeled via the Gaussian mixture model
(GMM) includes the possibility of the likelihood being either multi-modal or
uni-modal. We demonstrated the shortcomings of the LS model when dealing with
data sets that are in tension, using a toy model of PDFs. We found that, with
a single Gaussian distribution assumed in the LS method, neither the Monte
Carlo method nor the Hessian method can provide an accurate sampling of the
probability distribution of the pseudo PDFs. Although the Monte Carlo method
directly samples the probability distribution in the data space, the
assumption that there is only one Gaussian distribution in the PDF space leads
to a LS fit that would be stuck in a narrow valley of $\chi^{2}$. The same
situation applies to the Hessian method. Inflating the Hessian uncertainty
with certain tolerance criteria is a conservative prescription for the
uncertainty estimation currently used by PDF fitting collaborations. We find
that doing so is not a guarantee of a good approximation of the combined
likelihood of two data sets that are in tension, particularly in the case
where the tension is complicated. Specifically, the value of the tolerance
chosen is not based on statistically robust arguments. Further, we showed how
increasing the tolerance does not accurately capture tension and leads to
inadequate uncertainty in some regions of the PDF but at the same time an
overestimation of uncertainty in other regions of the PDF, c.f. Fig. 4. This
can also be seen in the correlation ellipses of the fit parameter space. What
is required here is not just an extension of the original $\chi^{2}$ ellipse,
but also a rotation of the axes of the correlation ellipse. The latter can
neither be achieved with a modified tolerance criteria nor by using a Monte-
Carlo sampling technique.

The novel GMM pseudo-PDF fit overcomes these shortcomings of the LS method by
parametrizing the likelihood as a mixture of Gaussian distributions, by
allowing for a multi-modal feature of the likelihood. We have described how
this unsupervised machine learning technique can be leveraged to learn about
tension in the fits as well as be used to determine precisely the
uncertainties on the PDFs. Specifically, it is possible to use the GMM for
performing the unsupervised machine learning task of clustering the data in
order to identify and learn about the tension in the data sets. We have
demonstrated how the GMM is simply generalization of the $\chi^{2}$ log-
likelihood and reduces to it when data is consistent. Finally, It is possible
to over-fit in the GMM, where the model might end up learning statistical
fluctuations rather than real tension between data sets. We have described the
use of the Akaike information criterion (AIC) and Bayesian information
criterion (BIC) to avoid over-fitting and can be used to determine optimal
value of the number of base distributions for the GMM method.

## Acknowledgment

MY thanks Yandong Liu and Pavel Nadolsky for the discussions. MY has left PKU
at the time of this paper being posted. This work was supported in part by the
National Science Foundation of China under grants No. 12075251, and the
Natural Science Foundation of Hunan province of China under Grant No.
2023JJ30496. It is also supported in part by the U.S. National Science
Foundation under Grant No. PHY-2310291 and PHY-2310497.

## References

  * [1] T.-J. Hou et al., “New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC,” Phys. Rev. D 103 no. 1, (2021) 014013, arXiv:1912.10053 [hep-ph].
  * [2] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne, “Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs,” Eur. Phys. J. C 81 no. 4, (2021) 341, arXiv:2012.04684 [hep-ph].
  * [3] NNPDF Collaboration, R. D. Ball et al., “The path to proton structure at 1% accuracy,” Eur. Phys. J. C 82 no. 5, (2022) 428, arXiv:2109.02653 [hep-ph].
  * [4] S. Alekhin, J. Blümlein, S. Moch, and R. Placakyte, “Parton distribution functions, $\alpha_{s}$, and heavy-quark masses for LHC Run II,” Phys. Rev. D 96 no. 1, (2017) 014011, arXiv:1701.05838 [hep-ph].
  * [5] ATLAS Collaboration, G. Aad et al., “Determination of the parton distribution functions of the proton using diverse ATLAS data from $pp$ collisions at $\sqrt{s}=7$, 8 and 13 TeV,” Eur. Phys. J. C 82 no. 5, (2022) 438, arXiv:2112.11266 [hep-ex].
  * [6] NNPDF Collaboration, R. D. Ball et al., “An open-source machine learning framework for global analyses of parton distributions,” Eur. Phys. J. C 81 no. 10, (2021) 958, arXiv:2109.02671 [hep-ph].
  * [7] S. Alekhin et al., “HERAFitter,” Eur. Phys. J. C 75 no. 7, (2015) 304, arXiv:1410.4412 [hep-ph].
  * [8] S. Amoroso et al., “Snowmass 2021 whitepaper: Proton structure at the precision frontier,” Acta Phys. Polon. B 53 no. 12, (2022) A1, arXiv:2203.13923 [hep-ph].
  * [9] CDF Collaboration, T. Aaltonen et al., “High-precision measurement of the W boson mass with the CDF II detector,” Science 376 no. 6589, (2022) 170–176.
  * [10] M. Cepeda et al., “Report from Working Group 2: Higgs Physics at the HL-LHC and HE-LHC,” CERN Yellow Rep. Monogr. 7 (2019) 221–584, arXiv:1902.00134 [hep-ph].
  * [11] LHC Higgs Cross Section Working Group Collaboration, D. de Florian et al., “Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector,” arXiv:1610.07922 [hep-ph].
  * [12] ATLAS Collaboration, M. Aaboud et al., “Precision measurement and interpretation of inclusive $W^{+}$ , $W^{-}$ and $Z/\gamma^{*}$ production cross sections with the ATLAS detector,” Eur. Phys. J. C 77 no. 6, (2017) 367, arXiv:1612.03016 [hep-ex].
  * [13] T.-J. Hou, H.-W. Lin, M. Yan, and C. P. Yuan, “Impact of lattice $s(x)-\bar{s}(x)$ data in the CTEQ-TEA global analysis,” arXiv:2204.07944 [hep-ph].
  * [14] M. Guzzi et al., “NNLO constraints on proton PDFs from the SeaQuest and STAR experiments and other developments in the CTEQ-TEA global analysis,” SciPost Phys. Proc. 8 (2022) 005, arXiv:2108.06596 [hep-ph].
  * [15] T.-J. Hou, M. Yan, J. Liang, K.-F. Liu, and C. P. Yuan, “Connected and disconnected sea partons from the CT18 parametrization of PDFs,” Phys. Rev. D 106 no. 9, (2022) 096008, arXiv:2206.02431 [hep-ph].
  * [16] K. Kovařík, P. M. Nadolsky, and D. E. Soper, “Hadronic structure in high-energy collisions,” Rev. Mod. Phys. 92 no. 4, (2020) 045003, arXiv:1905.06957 [hep-ph].
  * [17] NNPDF Collaboration, R. D. Ball, A. Candido, J. Cruz-Martinez, S. Forte, T. Giani, F. Hekhorn, K. Kudashkin, G. Magni, and J. Rojo, “Evidence for intrinsic charm quarks in the proton,” Nature 608 no. 7923, (2022) 483–487, arXiv:2208.08372 [hep-ph].
  * [18] R. Stegeman, “The negligible impact of experimental inconsistencies in the NNPDF4.0 global dataset,” PoS ICHEP2022 (2022) 787, arXiv:2212.07703 [hep-ph].
  * [19] N. Lambri, “Optimized regression models for parton distribution functions determination using deep learning models,” 2020.  http://nnpdf.mi.infn.it/wp-content/uploads/2021/01/thesis_lambri.pdf.
  * [20] L. Talon, “Optimization of parton density uncertainties,” 2019.  http://nnpdf.mi.infn.it/wp-content/uploads/2019/06/Talon_MS_thesis.pdf.
  * [21] PDF4LHC21 combination group Collaboration, T. Cridge, “PDF4LHC21: Update on the benchmarking of the CT, MSHT and NNPDF global PDF fits,” SciPost Phys. Proc. 8 (2022) 101, arXiv:2108.09099 [hep-ph].
  * [22] H. Paukkunen and P. Zurita, “PDF reweighting in the Hessian matrix approach,” JHEP 12 (2014) 100, arXiv:1402.6623 [hep-ph].
  * [23] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” Mathematics of control, signals, and systems 2 no. 4, (1989) 303–314.
  * [24] M. H. Stone, “Applications of the theory of boolean rings to general topology,” Transactions of the American Mathematical Society 41 no. 3, (1937) 375–481.
  * [25] M. H. Stone, “The generalized weierstrass approximation theorem,” 1948.
  * [26] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [27] J. Erler and R. Ferro-Hernández, “Alternative to the application of PDG scale factors,” Eur. Phys. J. C 80 no. 6, (2020) 541, arXiv:2004.01219 [physics.data-an].
  * [28] G. Cowan, “Statistical Models with Uncertain Error Parameters,” Eur. Phys. J. C 79 no. 2, (2019) 133, arXiv:1809.05778 [physics.data-an].
  * [29] G. D’Agostini, “Sceptical combination of experimental results: General considerations and application to epsilon-prime / epsilon,” arXiv:hep-ex/9910036.
  * [30] K. Pearson, “Contributions to the mathematical theory of evolution,” Philosophical Transactions of the Royal Society of London. A 185 (1894) 71–110.
  * [31] K. P. Murphy, Machine learning: a probabilistic perspective.  MIT Press, Cambridge, MA, 2012.
  * [32] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems,” arXiv:1603.04467.
  * [33] K. Basford, D. Greenway, G. Mclachlan, and D. Peel, “Standard errors of fitted component means of normal mixtures,” Computational Statistics 12 (01, 1997) 1.
  * [34] H. M. Bücker and G. F. Corliss, “A bibliography on automatic differentiation,” in Automatic Differentiation: Applications, Theory, and Implementations, H. M. Bücker, G. F. Corliss, P. D. Hovland, U. Naumann, and B. Norris, eds., vol. 50 of Lecture Notes in Computational Science and Engineering, pp. 321–322.  Springer, New York, NY, 2005.
  * [35] A. Courtoy, J. Huston, P. Nadolsky, K. Xie, M. Yan, and C. P. Yuan, “Parton distributions need representative sampling,” arXiv:2205.10444 [hep-ph].
  * [36] H.-L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin, and C. P. Yuan, “New parton distributions for collider physics,” Phys. Rev. D 82 (2010) 074024, arXiv:1007.2241 [hep-ph].
  * [37] J. Gao, M. Guzzi, J. Huston, H.-L. Lai, Z. Li, P. Nadolsky, J. Pumplin, D. Stump, and C. P. Yuan, “CT10 next-to-next-to-leading order global analysis of QCD,” Phys. Rev. D 89 no. 3, (2014) 033009, arXiv:1302.6246 [hep-ph].
  * [38] A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt, “Parton distributions for the LHC,” Eur. Phys. J. C 63 (2009) 189–285, arXiv:0901.0002 [hep-ph].
  * [39] L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, “Parton distributions in the LHC era: MMHT 2014 PDFs,” Eur. Phys. J. C 75 no. 5, (2015) 204, arXiv:1412.3989 [hep-ph].
  * [40] X. Jing et al., “Quantifying the interplay of experimental constraints in analyses of parton distributions,” Phys. Rev. D 108 no. 3, (2023) 034029, arXiv:2306.03918 [hep-ph].
  * [41] H. Akaike, “A new look at the statistical model identification,” IEEE Transactions on Automatic Control 19 no. 6, (1974) 716–723.
  * [42] G. Schwarz, “Estimating the Dimension of a Model,” The Annals of Statistics 6 no. 2, (1978) 461 – 464. https://doi.org/10.1214/aos/1176344136.