# Bayesian Approach to Inverse Problems: an Application to NNPDF Closure
Testing

Luigi Del Debbio  Higgs Centre for Theoretical Physics, School of Physics and
Astronomy, Peter Guthrie Tait Road, Edinburgh EH9 3 FD, United Kingdom.
Tommaso Giani  Department of Physics and Astronomy, Vrije Universiteit, NL
1081 HV Amsterdam  Nikhef Theory Group, Science Park 105, 1098 XG Amsterdam,
The Netherlands  Michael Wilson  Higgs Centre for Theoretical Physics, School
of Physics and Astronomy, Peter Guthrie Tait Road, Edinburgh EH9 3 FD, United
Kingdom.

###### Abstract

We discuss the Bayesian approach to the solution of inverse problems and apply
the formalism to analyse the closure tests performed by the NNPDF
collaboration. Starting from a comparison with the approach that is currently
used for the determination of parton distributions (PDFs) by the NNPDF
collaboration, we discuss some analytical results that can be obtained for
linear problems and use these results as a guidance for the more complicated
non-linear problems. We show that, in the case of Gaussian distributions, the
posterior probability density of the parametrized PDFs is fully determined by
the results of the NNPDF fitting procedure. In the particular case that we
consider, the fitting procedure and the Bayesian analysis yield exactly the
same result. Building on the insight that we obtain from the analytical
results, we introduce new estimators to assess the statistical faithfulness of
the fit results in closure tests. These estimators are defined in data space,
and can be studied analytically using the Bayesian formalism in a linear model
in order to clarify their meaning. Finally we present numerical results from a
number of closure tests performed with current NNPDF methodologies. These
further tests allow us to validate the NNPDF4.0 methodology and provide a
quantitative comparison of the NNPDF4.0 and NNPDF3.1 methodologies. As PDFs
determinations move into precision territory, the need for a careful
validation of the methodology becomes increasingly important: the error bar
has become the focal point of contemporary PDFs determinations. In this
perspective, theoretical assumptions and other sources of error are best
formulated and analysed in the Bayesian framework, which provides an ideal
language to address the precision and the accuracy of current fits.

##  1 Introduction

Inverse problems are the typical example of inference where a model is sought
starting from a finite-dimensional set of experimental observations. These
problems are notoriously difficult, and often require trying to guess a
continuous function, i.e. an element of an infinite dimensional space, from a
finite amount of data. As emphasised by Hadamard a long time ago, it is easy
to end up in a situation where we deal with ill-posed problems, in the sense
that the solution may not exist, may not be unique, or may be unstable under
small variations of the input data. The determination of parton distributions
from experimental data, or the reconstruction of spectral densities from
lattice QCD simulations, are just two examples where these problems arise in
particle physics. In all these cases, finding a robust solution becomes very
challenging, if not impossible, making these questions all the more urgent,
especially at a time when precision studies are the ultimate challenge in
order to highlight divergences from the Standard Model.

A Bayesian approach provides an apter tool for addressing inverse problems.
Starting from a prior distribution for the model, which encodes our
theoretical knowledge or prejudice, basic statistical tools allow us to
determine the posterior distribution of the solution, after taking into
account a set of experimental observations. The prior and posterior
probabilities encode our knowledge about the model before and after
incorporating the experimental results. There are multiple advantages to this
formulation: the inverse problem is well defined, and the prior distribution
ensures that all the assumptions are explicit, while regulating the likelihood
distribution.

Note also that probability measures can be defined in infinite-dimensional
spaces. In cases where we are looking to determine a continuous function, the
Bayesian approach provides, at least in principle, a formulation of the
problem directly in terms of the posterior probability measure of the model
function. It is often convenient for practical reasons to define a
parametrization of the model function in terms of a finite, albeit large,
number of parameters and reduce the problem to a finite dimensional one. The
price to pay for this simplification is the introduction of some bias, which
needs to be understood and possibly quantified. An expressive parametrization
clearly helps in this case.

A Bayesian approach to inverse problems has been actively developed by
mathematicians for a long time, and this development has dramatically
accelerated in the last decade. In this paper we aim at summarising the
existing framework and adapt it to analyse the fits of parton distribution
functions obtained by the NNPDF collaboration. We review the Bayesian
formalism in Sec. 2, where we define the notation that will be used in the
rest of the paper. We report some known examples for illustrative purposes.
Even though these are well-known results, we find it useful to summarise them
in the context of our specific problem. In Sec. 3, we try to connect the
Bayesian approach with the NNPDF fits based on a Monte Carlo methodology,
where the distribution of the PDFs is encoded in a set of fits to artificially
generated data, called replicas. We can anticipate here that, under the
hypotheses that the data are Gaussian and the model is linear, the NNPDF
procedure does characterise completely the posterior proabillity density. When
the model is non-linear two modifications need to be taken into account. First
of all the analytical calculation that we present in Sec. 2.1 is no longer
possible, and one needs to rely on the fact that a linearization of the model
yields an accurate prediction. Even though linear models are known to provide
good approximations, the systematic errors introduced by this approximation
are not easy to quantify. There is also a more subtle effect that needs to be
taken into account. When working with linear models, the minimization
procedure is know to have a unique minimum, which can be computed exactly.
Non-linear models can be plagued by multiple minima, and more importantly by
inefficiencies of the algorithm in finding them. While it would be foolish to
ignore the limitations of the analytical calculation, it is nonetheless very
useful to have an explicit result as a template to guide the analysis of our
numerical investigations.

In the NNPDF fitting framework, the posterior probability of the parton
distribution functions, is encoded in an ensemble of fits to replicas of the
data, where the data replicas have been generated in order to reproduce the
fluctuations described by the experimental central values and uncertainties.
This bootstrap procedure propagates the data fluctuations in the space of
fitted PDFs. A successful fit should yield robust confidence intervals for
observables, in particular for those that are more relevant for phenomenology.

The idea of closure tests is to test this procedure in a fit to artificial
data that have been generated from a known set of input PDFs. In this case the
underlying law is known and we can check how the posterior distribution
compares to the underlying law. This is the basis of a closure test, which is
summarised at the end of Sec. 3. Closure tests have already been used to test
to validity of previous iterations of the NNPDF methodology. Here we aim to
refine some of the pre-existing closure test estimators and with the help of
fast fitting methodology perform a more extensive study of how faithful our
uncertainties are. For this purpose we introduce in Sec. 4 new estimators that
allow us to quantify the quality of our fits. These estimators are defined in
the space of data and need to be understood as stochastic variables that are
characterised by a probability density. Where possible, we use the Bayesian
formalism in order to compute analytically these probability densities and
compare with numerical results. In order to perform analytical calculations we
often need to make some simplifying assumptions. While it is incorrect to use
the analytical results for making quantitative predictions for the results of
the realistic case of a non-linear fit, the analytical results provide a guide
to the interpretation of the observed patterns.

The results of our numerical studies are summarised in Sec. 5, where we also
compare the NNPDF4.0 methodology introduced in the latest NNPDF fit [1] to the
methodology used by the collaboration in previous fits.

The Bayesian formalism, by providing posterior probability distributions,
paves the way to explore a number of issues. We highlight in our conclusions
some possible questions that we defer to future studies.

##  2 Inverse Problems

The problem of determining PDFs from a set of experimental data falls under
the general category of inverse problems, i.e. the problem of finding the
input to a given model knowing a set of observations, which are often finite
and noisy. In this section we are going to review the Bayesian formulation of
inverse problems. It is impossible to do justice to this vast subject here.
Instead we try to emphasise the aspects that are relevant for quantifying
uncertainties on PDF determinations.

###  2.1 Statement of the problem

The space of inputs is denoted by $X$, while $R$ denotes the space of
responses. The model is specified by a forward map

$\displaystyle G:~{}$ $\displaystyle X\to R$ $\displaystyle u\mapsto
r=G(u)\,,$ (1)

which associates a response $r\in R$ to the input $u\in X$, where we assume
that $X$ and $R$ are Banach spaces. 111Banach spaces are complete normed
vector spaces. We do not need to get into a more detailed discussion here, but
it is important to note that working in Banach spaces allows us to generalise
the results to infinite-dimensional spaces of functions. As an example we can
think of $u$ as being a PDF, i.e. a function defined on the interval $[0,1]$,
and $r$ a DIS structure function. The structure function is related to the PDF
by a factorization formula involving perturbative coefficient functions:

$\displaystyle r(x,Q^{2})=\int_{x}^{1}\frac{dz}{z}\,C(z,Q^{2})u(x/z,Q^{2})\,.$
(2)

Note that in this example the forward map maps one real function into another
real function. In this case the space of models $X$ is the space of continuous
functions defined on the interval $[0,1]$, which satisfy integrability
conditions. Even though this is an infinite-dimensional space, it is possible
to define a probability measure on such a space and construct a Bayesian
solution to the inverse problem. In current determinations of PDFs, the
functional form of the PDF is dictated by some kind of parametrization, with
different parametrizations being used by different collaborations. In all
cases, the space $X$ is a finite-dimensional space, ${}^{N_{\rm model}}$,
where $N_{\rm model}$ is the number of parameters. In the case of the NNPDF
fits discussed below, the weights of the neural networks are the parameters
that determine the functional form of the PDFs. Alternatively, one could think
of a finite-dimensional representation defined by the value of the PDF at
selected values of $x$, i.e. $u_{i}=u(x_{i})$ for $i=1,\ldots,N_{\rm model}$.
Depending on the context we will denote by $u$ either the function $u(x)$, or
the vector of real parameters that are used to determine the function $u(x)$.
Disambiguation should hopefully be straightforward.

Experiments will not have access to the full function $r$ but only to a subset
of $N_{\rm data}$ observations. In order to have a formal mathematical
expression that takes into account the fact that we have a finite number of
measurements, we introduce an observation operator

$\displaystyle O:~{}$ $\displaystyle R\to Y$ $\displaystyle r\mapsto y\,,$ (3)

where $y\in Y$ is a vector in a finite-dimensional space $Y$ of experimental
results, e.g. the value of the structure function for some values of the
kinematic variables $x$ and $Q^{2}$. In general we will assume that
$y\in^{N_{\rm data}}$, i.e. we have a finite number $N_{\rm data}$ of real
experimental values. The quantity of interest is the composed operator

$\displaystyle\mathcal{G}:~{}$ $\displaystyle X\to^{N_{\rm data}}$
$\displaystyle\mathcal{G}=O\circ G\,,$ (4)

which maps the input $u$ to the set of data. Taking into account the fact that
experimental data are subject to noise, we can write

$\displaystyle y=\mathcal{G}(u)+\eta\,,$ (5)

where $\eta$ is a random variable defined over ${}^{N_{\rm data}}$ with
probability density $\rho(\eta)$. We will refer to $\eta$ as the observational
noise. In this setting, the inverse problem becomes finding $u$ given $y$. It
is often the case that inverse problems are ill-defined in the sense that the
solution may not exist, may not be unique, or may be unstable under small
variations of the problem.

In solving the inverse problem, we are going to adopt a Bayesian point of
view, as summarised e.g. in Ref. [2]: our prior knowledge about $u$ is encoded
in a prior probability measure $\mu_{X}^{0}$, where the suffix $X$ indicates
that the measure is defined in the space of models, and the suffix 0 refers to
the fact that this is a prior distribution. We use Bayes’ theorem to compute
the posterior probability measure of $u$ given the data $y$, which we denote
as $\mu_{X}^{\mathcal{G}}$. When the probability measure can be described by a
probability density, we denote the probability densities associated to
$\mu_{X}^{0}$ and $\mu_{X}^{\mathcal{G}}$, by $\pi_{X}^{0}$ and
$\pi_{X}^{\mathcal{G}}$ respectively. Then, using Eq. (5), we can write the
data likelihood, i.e. the probability density of $y$ given $u$,

$\displaystyle\pi_{Y}(y|u)=\rho(y-\mathcal{G}(u))\,,$ (6)

and Bayes’ theorem yields

$\displaystyle\pi_{X}^{\mathcal{G}}(u)=\pi_{X}(u|y)\propto\pi_{X}^{0}(u)\rho(y-\mathcal{G}(u))\,.$
(7)

Even though the concepts that we have introduced so far should sound familiar,
it is worthwhile clarifying some ideas and present an explicit case where all
the probability densities are carefully defined. This is best exemplified by
considering the case where both the observational noise and the model prior
are Gaussian. We assume that we are given a set of central values
$y_{0}\in^{N_{\rm data}}$ and their covariance matrix $C_{Y}$. Then the prior
probability density of the observable $y$ is

$\pi_{Y}^{0}(y|y_{0},C_{Y})\propto\exp\left(-\frac{1}{2}\left|y-y_{0}\right|_{C_{Y}}^{2}\right)\,,$
(8)

where, similarly to the convention used above, the suffix $Y$ emphasises the
fact that this is a probability density in data space, and the notation
explicitly reminds us that this is the probability density given the central
values $y_{0}$ (and the covariance matrix $C_{Y}$). Similarly we can choose a
Gaussian distribution for the prior distribution of the input model,
characterized by a central value $u_{0}$ and a covariance $C_{X}$:

$\displaystyle\pi_{X}^{0}(u|u_{0},C_{X})$
$\displaystyle\propto\exp\left(-\frac{1}{2}\left|u-u_{0}\right|_{C_{X}}^{2}\right)\,.$
(9)

Following the convention above, we use a suffix $X$ here to remind the reader
that we are looking at a probability density in the space of models. Note that
in the expressions above we used the norms in $X$ and ${}^{N_{\rm data}}$
respectively, and introduced the short-hand notation

$\displaystyle\left|a\right|_{M}^{2}=\left|M^{-1/2}a\right|^{2}\,,$ (10)

where $a$ denotes a generic element of $X$, $R$ or ${}^{N_{\rm data}}$. For
the case where $a\in^{N_{\rm data}}$, we use the Euclidean norm and

$\displaystyle\left|a\right|_{M}^{2}=\sum_{i,j}a_{i}M^{-1}_{ij}a_{j}\,,$ (11)

where the indices $i,j$ run from 1 to $N_{\rm data}$, which eventually yields
the usual expression for the $\chi^{2}$ of correlated data. Up to this point
data and models are completely independent, and the joint distribution is
simply the product of $\pi_{Y}^{0}$ and $\pi_{X}^{0}$.

The forward map induces a correlation between the input model and the
observables, so we introduce a probability density $\theta$ that describes
these correlations due to the underlying theory,

$\theta(y,u|\mathcal{G})=\delta\left(y-\mathcal{G}(u)\right)\,,$ (12)

where the Dirac delta corresponds to the case where there are no theoretical
uncertainities. Theoretical uncertainties can be introduced by broadening the
distribution of $y$ away from the exact prediction of the forward map, e.g.
using a Gaussian with covariance $C_{T}$,

$\theta(y,u|\mathcal{G})=\exp\left(-\frac{1}{2}\left|y-\mathcal{G}(u)\right|_{C_{T}}^{2}\right)\,.$
(13)

In the context of PDF fitting a similar recipe to take into account
theoretical errors has recently been advocated in Refs. [3, 4]. Note that
there are no rigorous arguments favouring the assumption that theoretical
errors are normally distributed; it is nonetheless a useful working
assumption, and a definite improvement compared to ignoring the theoretical
errors altogether. The net effect of the theory errors is a redefinition of
the covariance of the data, which has no major impact in our discussion, and
therefore will be ignored. Taking the correlation $\theta(y,u|\mathcal{G})$
into account, the joint distribution of $y$ and $u$ is

$\displaystyle\pi^{\mathcal{G}}(y,u|y_{0},C_{Y},u_{0},C_{X})\propto\pi_{X}^{0}(u|u_{0},C_{X})\pi_{Y}^{0}(y|y_{0},C_{Y})\theta(y,u|\mathcal{G})\,.$
(14)

We can now marginalize with respect to y, neglecting theory errors,

$\displaystyle\pi^{\mathcal{G}}_{X}(u|y_{0},C_{Y},u_{0},C_{X})$
$\displaystyle\propto\int
dy\,\pi_{X}^{0}(u|u_{0},C_{X})\pi_{Y}^{0}(y|y_{0},C_{Y})\theta(y,u|\mathcal{G})$
(15) $\displaystyle\propto\pi_{X}^{0}(u|u_{0},C_{X})\int
dy\,\pi_{Y}^{0}(y|y_{0},C_{Y})\delta\left(y-\mathcal{G}(u)\right)$ (16)
$\displaystyle\propto\pi_{X}^{0}(u|u_{0},C_{X})\,\pi_{Y}^{0}(\mathcal{G}(u)|y_{0},C_{Y})\,.$
(17)

We see that we have recovered Eq. 7. The log-likelihood in the Gaussian case
is simply the $\chi^{2}$ of the data, $y_{0}$, to the theory prediction,
$\mathcal{G}(u)$:

$-\log\pi_{Y}^{0}(\mathcal{G}(u)|y_{0},C_{Y})=\frac{1}{2}\sum_{i,j=1}^{N_{\rm
data}}\left(\mathcal{G}(u)-y_{0}\right)_{i}\left(C_{Y}^{-1}\right)_{ij}\left(\mathcal{G}(u)-y_{0}\right)_{j}\,.$
(18)

In the notation of Eq. 7

$\pi_{Y}^{0}(\mathcal{G}(u)|y_{0},C_{Y})=\rho\left(\mathcal{G}(u)-y_{0}\right)\,,$
(19)

where in this case

$\displaystyle\rho(\eta)$
$\displaystyle\propto\exp\left(-\frac{1}{2}\left|\eta\right|_{C_{Y}}^{2}\right)\,.$
(20)

The probability density $\pi^{\mathcal{G}}_{X}(u|y_{0},C_{Y},u_{0},C_{X})$ was
called $\pi_{X}^{\mathcal{G}}(u)$ in Eq. 7, where the suffix $\mathcal{G}$ is
a short-hand to denote the posterior probability in model space, taking into
account all the conditional variables. Hence, for the Gaussian case, the
result from Bayes’ theorem reduces to

$\displaystyle\pi_{X}^{\mathcal{G}}(u)$
$\displaystyle\propto\exp\left[-\frac{1}{2}\left|y_{0}-\mathcal{G}(u)\right|_{C_{Y}}^{2}-\frac{1}{2}\left|u-u_{0}\right|_{C_{X}}^{2}\right]$
(21) $\displaystyle\propto\exp\left[-S(u)\right]\,.$ (22)

Note that in the argument of the likelihood function we have the central
values of the data points $y_{0}$ as given by the experiments. Eq. (21) is the
Bayesian answer to the inverse problem, our knowledge of the model $u$ is
encoded in the probability measure $\mu_{X}^{\mathcal{G}}$, which is fully
specified by the density $\pi_{X}^{\mathcal{G}}$. There are several ways to
characterise a probability distribution, a task that becomes increasingly
difficult in high-dimensional spaces. As discussed later in this study, the
NNPDF approach is focused on the determination of the Maximum A Posteriori
(MAP) estimator, i.e. the element $u_{*}\in X$ that maximises
$\pi_{X}^{\mathcal{G}}(u)$:

$\displaystyle u_{*}=\arg\min_{u\in
X}\left(\frac{1}{2}\left|y_{0}-\mathcal{G}(u)\right|_{C_{Y}}^{2}+\frac{1}{2}\left|u-u_{0}\right|_{C_{X}}^{2}\right)\,.$
(23)

For every instance of the data $y_{0}$, the MAP estimator is computed by
minimising a regulated $\chi^{2}$, where the regularization is determined by
the prior that is assumed on the model $u$. We will refer to this procedure as
the classical fit of experimental data to a model. Note that in the Bayesian
approach, the regulator appears naturally after having specified carefully all
the assumptions that enter in the prior. In this specific example the
regulator arises from the Gaussian prior for the model input $u$, which is
normally distributed around a solution $u_{0}$. The MAP estimator provides the
explicit connection between the Bayesian approach and the classical fit.

###  2.2 Comparison with classical fitting

Analytical results make the connection between the two approaches more
quantitative, and therefore more transparent. We are going to summarise these
results here without proofs, referring the reader to the mathematical
literature for the missing details. Working in the finite-dimensional case, we
assume

$\displaystyle u$ $\displaystyle\in^{N_{\rm model}}\,,$ $\displaystyle y$
$\displaystyle\in^{N_{\rm data}}\,,$

and we are going to review in detail two examples from Ref. [2], which
illustrate the role of the priors in the Bayesian setting. We report here the
results in Ref. [2] because their simplicity provides a striking example of
the role of priors, which is sometimes underestimated. It is particularly
useful to distinguish the case of an underdetermined system from the case of
an overdetermined one.

#### Underdetermined system

The first case that we are going to analyse is the case of a linear system
that is underdetermined by the data. The linear model is completely specified
by a vector of coefficients $g\in^{N_{\rm model}}$,

$\mathcal{G}(u)=\left(g^{T}u\right)\,.$ (24)

Assuming that we have only one datapoint, i.e. $N_{\rm data}=1$,

$y=(g^{T}u)+\eta\,,$ (25)

where $\eta\sim\mathcal{N}(0,\gamma^{2})$ is one Gaussian number, whose
probability density is centred at $0$ and has variance $\gamma^{2}$.

$X_{\parallel}$$X_{\perp}$$u$$g$$y/|g|$ Figure 1: A linear model in a two-
dimensional space is constrained by a single data point $y$. Once the vector
$g$ is given, the model is fully specified by the vector $u$, viz. $y=g^{T}u$.
All models $u$ along the vertical red line reproduce exactly the data point.
In the absence of prior knowledge, every model along that one-dimensional
subspace is a legitimate solution of the inverse problem, which is clearly
underdetermined.

For simplicity we are going to assume that the prior on $u$ is also a multi-
dimensional Gaussian, centred at $0$ with covariance matrix $C_{X}$. In this
case the posterior distribution can be written as

$\pi_{X}^{\mathcal{G}}(u)\propto\exp\left[-\frac{1}{2\gamma^{2}}\left|y-(g^{T}u)\right|^{2}-\frac{1}{2}\left|u\right|_{C_{X}}^{2}\right]\,,$
(26)

which is still a Gaussian distribution for $u$. The mean and covariance are
respectively

$\displaystyle m$
$\displaystyle=\frac{(C_{X}g)y}{\gamma^{2}+(g^{T}C_{X}g)}\,,$ (27)
$\displaystyle\Sigma$
$\displaystyle=C_{X}-\frac{(C_{X}g)(C_{X}g)^{T}}{\gamma^{2}+(g^{T}C_{X}g)}\,.$
(28)

Because the argument of the exponential is a quadratic form, the mean of the
distribution coincides with the MAP estimator. Hence in this case a fit that
minimises the $\chi^{2}$ of the data to the theory prediction yields exactly
the mean of the posterior distribution. It is instructive to look at these
quantities in the limit of infinitely precise data, i.e. in the limit
$\gamma\to 0$:

$\displaystyle m_{\star}$ $\displaystyle=\lim_{\gamma\to
0}m=\frac{(C_{X}g)y}{(g^{T}C_{X}g)}\,,$ (29) $\displaystyle\Sigma_{\star}$
$\displaystyle=\lim_{\gamma\to
0}\Sigma=C_{X}-\frac{(C_{X}g)(C_{X}g)^{T}}{(g^{T}C_{X}g)}\,.$ (30)

These values satisfy

$\displaystyle(g^{T}m_{\star})=y\,,$ (31)
$\displaystyle(\Sigma_{\star}g)=0\,,$ (32)

which shows that the mean of the distribution is such that the data point is
exactly reproduced by the model, and that the uncertainty in the direction
defined by $g$ vanishes. It should be noted that the uncertainty in directions
perpendicular to $g$ does not vanish and is determined by a combination of the
prior and the model, viz. $C_{X}$ and $g$ in our example. This is a particular
example of a more general feature: for underdetermined systems the information
from the prior still shapes the probability distribution of the solution even
in the limit of vanishing statistical noise.

It is interesting to analyse what happens when the prior on the model is
removed. For this purpose we can start from a covariance in model space that
is proportional to the identity matrix,

$C_{X}=\sigma\mathbb{1}\,,$ (33)

and take the limit $\sigma\to\infty$. In this limit

$m_{\star}=y\frac{g}{(g^{T}g)}\,,$ (34)

while the posterior covariance becomes

$\Sigma_{\star}=\sigma\left(\mathbb{1}-\frac{gg^{T}}{(g^{T}g)}\right)\,.$ (35)

This covariance is the projector on the subspace orthogonal to $g$ multiplied
by $\sigma$. In the simple two-dimensional example depicted in Fig. 1,
choosing the direction of $g$ as one of the basis vector, we obtain

$\Sigma_{\star}=\begin{pmatrix}0&\\\ &\sigma\end{pmatrix}\,,$ (36)

which shows that the error in the $g$ direction vanishes, while the error in
the orthogonal subspace diverges as we remove the prior.

#### Overdetermined system

We are now going to consider an example of an overdetermined system and
discuss again the case of small observational noise. We consider $N_{\rm
data}\geq 2$ and $N_{\rm model}=1$, with a linear forward map such that

$y=gu+\eta\,,$ (37)

where $\eta$ is an $N_{\rm data}$-dimensional Gaussian variable with a
diagonal covariance $\gamma^{2}\mathbb{1}$, and $\mathbb{1}$ denotes the
identity matrix. For simplicity we are going to assume a Gaussian prior with
unit variance for $u$, which yields for the posterior distribution:

$\pi_{X}^{\mathcal{G}}(u)\propto\exp\left(-\frac{1}{2\gamma^{2}}\left|y-gu\right|^{2}-\frac{1}{2}u^{2}\right)\,.$
(38)

The posterior is Gaussian and we can easily compute its mean and variance:

$\displaystyle m$ $\displaystyle=\frac{(g^{T}y)}{\gamma^{2}+|g|^{2}}\,,$ (39)
$\displaystyle\sigma^{2}$
$\displaystyle=\frac{\gamma^{2}}{\gamma^{2}+|g|^{2}}\,.$ (40)

In this case, in the limit of vanishing observational noise, we obtain

$\displaystyle m_{\star}$ $\displaystyle=\frac{(g^{T}y)}{|g|^{2}}\,,$ (41)
$\displaystyle\sigma_{\star}^{2}$ $\displaystyle=0\,.$ (42)

The mean is given by the weighted average of the datapoints, which is also the
solution of the $\chi^{2}$ minimization

$m_{\star}=\arg\min_{u\in}\left|y-gu\right|^{2}\,.$ (43)

Note that in this case the variance $\sigma_{\star}$ vanishes independently of
the prior. In the limit of small noise, the distribution tends to a Dirac
delta around the value of the MAP estimator.

###  2.3 Linear Problems

Linear problems in finite-dimensional spaces are characterized by a simple
forward law,

$y=\mathcal{G}u\,,$ (44)

where $\mathcal{G}$ is a matrix. In this framework one can readily derive
analytical solutions that are useful to understand the main features of the
Bayesian approach. Assuming that the priors are Gaussian again, the cost
function $S(u)$ is a quadratic function of $u$,

$\displaystyle S(u)$
$\displaystyle=\left(\mathcal{G}u-y_{0}\right)^{T}C_{Y}^{-1}\left(\mathcal{G}u-y_{0}\right)+\left(u-u_{0}\right)^{T}C_{X}^{-1}\left(u-u_{0}\right)$
(45)
$\displaystyle=\left(u-\tilde{u}\right)\tilde{C}_{X}^{-1}\left(u-\tilde{u}\right)+\mathrm{const}\,,$
(46)

where

$\displaystyle\tilde{C}_{X}^{-1}$
$\displaystyle=\left(\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G}+C_{X}^{-1}\right)\,,$
(47) $\displaystyle\tilde{u}$
$\displaystyle=\tilde{C}_{X}\left(\mathcal{G}^{T}C_{Y}^{-1}y_{0}+C_{X}^{-1}u_{0}\right)\,.$
(48)

The case where we have no prior information on the model is recovered by
taking the limit $C_{X}^{-1}\to 0$, which yields

$\displaystyle\tilde{C}_{X}^{-1}$
$\displaystyle=\left(\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G}\right)\,,$ (49)
$\displaystyle\tilde{u}$
$\displaystyle=\tilde{C}_{X}\left(\mathcal{G}^{T}C_{Y}^{-1}y_{0}\right)\,.$
(50)

The action of $C_{Y}^{-1}$ on the vector of observed data $y_{0}$ is best
visualised using a spectral decomposition

$C_{Y}^{-1}=\sum_{n}\frac{1}{\sigma_{n}^{2}}P_{n}\,,$ (51)

where $P_{n}$ denotes the projector on the $n$-th eigenvector of $C_{Y}$, and
$\sigma_{n}^{2}$ is the corresponding eigenvalue. The action of $C_{Y}^{-1}$
is to perform a weighted average of the components of $y_{0}$ in the
directions of the eigenvectors of $C_{Y}$.

An explicit expression for the posterior distribution of data can be obtained
from the joint distribution by marginalising over the model input $u$:

$\displaystyle\pi_{Y}^{\mathcal{G}}(y|y_{0},C_{Y},u_{0},C_{X})$
$\displaystyle=\int du\,\pi^{\mathcal{G}}(y,u|y_{0},C_{Y},u_{0},C_{X})$ (52)
$\displaystyle\propto\exp\left(-\frac{1}{2}\left(y-\tilde{y}\right)^{T}\tilde{C}_{Y}^{-1}\left(y-\tilde{y}\right)\right)\,,$
(53)

where

$\displaystyle\tilde{y}$ $\displaystyle=\mathcal{G}\tilde{u}\,,$ (54)
$\displaystyle\tilde{C}_{Y}$
$\displaystyle=\mathcal{G}\tilde{C}_{X}\mathcal{G}^{T}\,.$ (55)

Note that this is the naive error propagation from the covariance of the
model, $\tilde{C}_{X}$, to the space of data.

#### Posterior distribution of unseen data

In real-life cases we are also interested in the posterior distribution of a
set of data that have not been included in the fit. Because different datasets
are described by the same theory, the knowledge of one dataset will inform our
knowledge of the underlying theory – i.e. we will determine a posterior
distribution for the model. That new knowledge about the model will then
propagate to any other – unseen – set of data, even if the experiments are
completely unrelated. In the Bayesian framework that we have developed, this
situation can be modeled by having two independent sets of data $y$ and
$y^{\prime}$, for which we have a prior distribution

$\displaystyle\pi_{Y}^{0}$
$\displaystyle\left(y,y^{\prime}|y_{0},C_{Y},y^{\prime}_{0},C^{\prime}_{Y}\right)=\pi_{Y}^{0}\left(y^{\prime}|y^{\prime}_{0},C^{\prime}_{Y}\right)\pi_{Y}^{0}\left(y|y_{0},C_{Y}\right)$
(56)
$\displaystyle\propto\exp\left[-\frac{1}{2}\left(y^{\prime}-y^{\prime}_{0}\right)^{T}(C^{\prime}_{Y})^{-1}\left(y^{\prime}-y^{\prime}_{0}\right)\right]\,\exp\left[-\frac{1}{2}\left(y-y_{0}\right)^{T}(C_{Y})^{-1}\left(y-y_{0}\right)\right]\,.$
(57)

Note that the prior distribution is factorised as the product of individual
distributions for $y$ and $y^{\prime}$ since the datasets are assumed to be
independent. Following the derivation above, we can write the joint
distribution for the data and the model

$\pi^{\mathcal{G}}(y,y^{\prime},u)\propto\pi_{Y}^{0}(y,y^{\prime}|y_{0},C_{Y},y^{\prime}_{0},C^{\prime}_{Y})\pi_{X}^{0}(u)\delta\left(y-\mathcal{G}u\right)\delta\left(y^{\prime}-\mathcal{G}^{\prime}u\right)\,.$
(58)

Because both sets of data are derived from the same model $u$, the joint
distribution above introduces a correlation between the data sets. The same
model $u$ appears in both delta functions in the equation above. We can now
marginalise with respect to the dataset $y$,

$\begin{split}\pi(y^{\prime},u)\propto&\exp\left[-\frac{1}{2}\left(y^{\prime}-y^{\prime}_{0}\right)^{T}(C^{\prime}_{Y})^{-1}\left(y^{\prime}-y^{\prime}_{0}\right)\right]\,\exp\left[-\frac{1}{2}\left(u-\tilde{u}\right)^{T}(\tilde{C}_{X})^{-1}\left(u-\tilde{u}\right)\right]\\\
&\quad\times\delta\left(y^{\prime}-\mathcal{G}^{\prime}u\right)\,.\end{split}$
(59)

where $\tilde{C}_{X}$ and $\tilde{u}$ are given respectively in Eqs. 47 and
48. By marginalising again, this time with respect to the model, we derive the
posterior distribution of the unseen data,

$\pi^{y}_{Y}(y^{\prime})\propto\exp\left[\frac{1}{2}\left(y^{\prime}-\tilde{y}^{\prime}\right)^{T}(\tilde{C}^{\prime}_{Y})^{-1}\left(y^{\prime}-\tilde{y}^{\prime}\right)\right]\,,$
(60)

where

$\displaystyle\tilde{C}^{\prime}_{Y}$
$\displaystyle=\mathcal{G}^{\prime}\tilde{C}^{\prime}_{X}\mathcal{G}^{\prime
T}$ (61) $\displaystyle\tilde{y}^{\prime}$
$\displaystyle=\mathcal{G}^{\prime}\tilde{u}^{\prime}\,,$ (62)

and we have introduced the variables $\tilde{u}^{\prime}$ and
$\tilde{C}^{\prime}_{X}$,

$\displaystyle\tilde{C}_{X}^{\prime-1}$ $\displaystyle=\mathcal{G}^{\prime
T}C_{Y}^{\prime-1}\mathcal{G}^{\prime}+\tilde{C}_{X}^{-1}$ (63)
$\displaystyle\tilde{u}^{\prime}$
$\displaystyle=\tilde{C}^{\prime}_{X}\left(\mathcal{G}^{\prime
T}C_{Y}^{\prime-1}y_{0}^{\prime}+\tilde{C}_{X}^{-1}\tilde{u}\right)\,.$ (64)

The variables $\tilde{C}_{X}$ and $\tilde{u}$ have been defined above. We
repeat their definition here in order to have all the necessary equations
collected together:

$\displaystyle\tilde{C}_{X}^{-1}$
$\displaystyle=\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G}+C_{X}^{-1}$ (65)
$\displaystyle\tilde{u}$
$\displaystyle=\tilde{C}_{X}\left(\mathcal{G}^{T}C_{Y}^{-1}y_{0}+C_{X}^{-1}u_{0}\right)\,.$
(66)

Eqs. 61, 62, 63, 64, 65 and 66 yield the central value and the variance of the
posterior distribution of unseen data as a function of $y_{0}$,
$y_{0}^{\prime}$, $C_{Y}$, $C_{Y}^{\prime}$, $u_{0}$ and $C_{X}$.

#### A comment on non-linear models

The linear models that we have discussed so far may look over-simplified at
first sight. In practice, it turns out that non-linear models can often be
linearised around the central value of the prior distribution,

$\mathcal{G}(u)=\mathcal{G}(u_{0})+G\left(u-u_{0}\right)+\ldots\,,$ (67)

where

$G^{i}_{\alpha}=\left.\frac{\partial\mathcal{G}^{i}}{\partial
u_{\alpha}}\right|_{u_{0}}\,,$ (68)

and we have neglected higher-order terms in the expansion of $\mathcal{G}(u)$.

If these terms are not negligible, another option is to find the MAP
estimator, and then expand the the forward map around it, which yields
equations very similar to the previous ones, with $u_{0}$ replaced by $u_{*}$.
If the posterior distribution of $u$ is sufficiently peaked around the MAP
estimator, then the linear approximation can be sufficiently accurate.

###  2.4 The infinite-dimensional case

In the finite-dimensional case, where the probability measures are specified
by their densities with respect to the Lebesgue measure, Eq. (7) can be
rephrased by saying that $\rho$ is the Radon-Nikodym derivative of the
probability measure $\mu^{\mathcal{G}}$ with respect to $\mu_{0}$, viz.

$\displaystyle\frac{d\mu^{\mathcal{G}}}{d\mu^{0}}(u)\propto\rho(y-\mathcal{G}(u))\,.$
(69)

Using the fact that the density $\rho$ is a positive function, we can rewrite

$\displaystyle\rho(y-\mathcal{G}(u))=\exp\left(-\Phi(u;y)\right)\,,$ (70)

and therefore

$\displaystyle\frac{d\mu^{\mathcal{G}}}{d\mu^{0}}(u)\propto\exp\left(-\Phi(u;y)\right)\,.$
(71)

In finite-dimensional spaces, the three equations above are just definitions
that do not add anything to the above discussion in terms of probability
densities. Their interest resides in the fact that the last expression, Eq.
(71), can be properly defined when $X$ is infinite-dimensional, allowing a
rigorous extension of the Bayesian formulation of inverse problems to the case
of infinite-dimensional spaces.

Summarising the details of probability measure in infinite-dimensional spaces,
is beyond the scope of this work. Adopting instead a heuristic approach, we
can say that a function $f$ is a random function if $f(x)$ is a random
variable for all values of $x$. Since the values of the function at different
values of $x$ can be correlated, a random function is fully characterised by
specifying the joint probability densities

$\pi\left(f_{1},\ldots,f_{n};x_{1},\ldots x_{n}\right)\,,$ (72)

where $f_{i}=f(x_{i})$, for all values of $n$, and all values of
$x_{1},\ldots,x_{n}$. This infinite set of finite-dimensional densities allows
the definition of a probability measure.

For a Gaussian random function, these densities only depend on a mean value
function $m(x)$ and a covariance $C(x,x^{\prime})$. The probability densities
for the variables $f_{i}$, for any value of $n$ is

$\pi\left(f_{1},\ldots,f_{n};x_{1},\ldots,x_{n}\right)\propto\exp\left[-\frac{1}{2}\sum_{ij}\left(f_{i}-m_{i}\right)C^{-1}(x_{i},x_{j})\left(f_{j}-m_{j}\right)\right]\,.$
(73)

The covariance $C$ is such that

$C(x,x^{\prime})=\int
df\,df^{\prime}\,\left(f-m(x)\right)\left(f^{\prime}-m(x^{\prime})\right)\pi\left(f,f^{\prime};x,x^{\prime}\right)\,,$
(74)

which shows that the two-point probability density determines all the other
distributions.

#### Functional linear problems

This formalism allows us to formulate a Bayesian solution of the inverse
problem

$y^{i}=\int dx\,G^{i}(x)u(x)\,,$ (75)

where $y^{i}$ is a discrete set of observables and $u(x)$ is a random
function, with a Gaussian prior with mean $u_{0}(x)$ and covariance
$C_{X}(x,x^{\prime})$. The vector of observed values is denoted $y_{0}$, and
we assume that the prior distribution of $y$ is a Gaussian centred at $y_{0}$
with covariance $C_{Y}$.

Similarly to the finite-dimensional case, the Bayesian solution yields a
Gaussian random function for the posterior distribution of the solution
$u(x)$. In order to characterise the posterior Gaussian distribution we need
explicit expressions for its mean and its covariance. Introducing the matrix

$S^{ij}=\int
dxdx^{\prime}\,G^{i}(x)C_{X}(x,x^{\prime})G^{j}(x^{\prime})+C_{Y}^{ij}\,,$
(76)

and its inverse $T^{ij}=\left(S^{-1}\right)^{ij}$, the posterior Gaussian
field is centred at

$\displaystyle\tilde{u}(x)=u_{0}(x)+\int
dx^{\prime}\,C_{X}(x,x^{\prime})G^{i}(x^{\prime})T^{ij}\left(y_{0}^{j}-\int
dx^{\prime\prime}\,G^{j}(x^{\prime\prime})u_{0}(x^{\prime\prime})\right)\,,$
(77)

which is the expected generalization of the finite-dimensional example
discussed above. Interestingly, this can be rewritten as

$\tilde{u}(x)=u_{0}(x)+\int
dx^{\prime}\,C_{X}(x,x^{\prime})\psi(x^{\prime})\,,$ (78)

where

$\displaystyle\psi(x)=G^{i}(x)\delta y^{i}\,,$ (79)

and the weighted residuals are given by

$\delta y^{i}=T^{ij}\left(y_{0}^{j}-\int
dx^{\prime}\,G^{j}(x^{\prime})u_{0}(x^{\prime})\right)\,.$ (80)

Defining

$\Psi^{i}(x)=\int dx^{\prime}\,C_{X}(x,x^{\prime})G^{i}(x^{\prime})\,,$ (81)

the posterior covariance can be written as

$\tilde{C}_{X}(x,x^{\prime})=C_{X}(x,x^{\prime})-\Psi^{i}(x)T^{ij}\Psi^{j}(x^{\prime})\,.$
(82)

It is instructive to compare the Bayesian result summarised above with the
method proposed by Backus and Gilbert [5] to solve the same inverse problem.
Assuming that there exists an unknown ’true’ model $u_{\mathrm{true}}$, such
that the observed data are

$y_{0}^{i}=\int dx\,G^{i}(x)u_{\mathrm{true}}(x)\,,$ (83)

we look for an estimate $u_{\mathrm{est}}$ of the true solution in the form

$u_{\mathrm{est}}(x)=Q^{i}(x)y^{i}_{0}\,,$ (84)

so that the problem is now recast as finding the functions $Q^{i}(x)$. Using
Eq. 83 we obtain

$u_{\mathrm{est}}(x)=\int
dx^{\prime}R(x,x^{\prime})u_{\mathrm{true}}(x^{\prime})\,,$ (85)

which states that with a finite amount of data we can only resolve a filtered
version of the true solution. The kernel $R$ is given by

$R(x,x^{\prime})=Q^{i}(x)G^{i}(x^{\prime})\,.$ (86)

The coefficient functions $Q^{i}(x)$ can be chosen so that the kernel is as
close as possible to a delta function,

$R(x,x^{\prime})\simeq\delta(x,x^{\prime})~{}~{}\Longrightarrow~{}~{}u_{\mathrm{est}}\simeq
u_{\mathrm{true}}\,.$ (87)

Approximating the delta function can be achieved by minimising

$\int dx^{\prime}\,\left(R(x,x^{\prime})-\delta(x-x^{\prime})\right)^{2}\,,$
(88)

which yields

$Q^{i}(x)=\left(S^{-1}\right)^{ij}G^{j}(x)\,,$ (89)

where

$S^{ij}=\int dx\,G^{i}(x)G^{j}(x)\,.$ (90)

The interesting observation is that the central value of the Bayesian solution
presented above reduces to the Backus-Gilbert $u_{\mathrm{est}}$ in the case
where $u_{0}$ is just white noise and therefore

$C_{X}(x,x^{\prime})=\delta(x-x^{\prime})\,.$ (91)

The connections between the Bayesian treatment and the Backus-Gilbert solution
and its regulated variations, deserves further investigations, which we defer
to future studies. Note that the Bayesian solution allows a variety of priors
to be explicitly declared and compared to the Backus-Gilbert solution.

##  3 NNPDF Monte Carlo approach to inverse problems

In this section we discuss the NNPDF approach to inverse problems, and make
contact explicitly with the formalism laid out in Sec. 2. In the Bayesian
formulation, Eq. 21 gives a quantitative description of how the information
contained in the experimental data propagates into our knowledge of the space
of models. In practice, it should be possible to sample directly from the
posterior distribution, or to find a solution using some kind of Backus-
Gilbert methodology [5]. These new types of approach are not straightforward
and we defer their investigation to further, dedicated studies. Here we focus
instead on the standard NNPDF fitting procedure and investigate its relation
with the Bayesian result. The NNPDF approach generates an ensemble of fit
results, which are supposed to describe the posterior probability distribution
for the model (i.e. in the space of PDFs) given the experimental data. In the
case of a linear map, we show here that this is exactly the case: the NNPDF
replicas are distributed exactly accroding to the posterior density that was
obtained in the previous section.

###  3.1 Fitting replicas

The approach for generating a sample in model space utilised by NNPDF can
broadly be described as fitting model replicas to pseudo-data replicas. As
discussed in Eq. 5 the experimental values are subject to observational noise.
If we assume this observational noise to be multigaussian then the
experimental central values, $y_{0}$, are given explicitly by

$y_{0}=f+\eta,$ (92)

where $f$ is the vector of true observable values, and the observational noise
is drawn from a Gaussian centred on zero such as in Eq. 20, i.e.
$\eta\sim\mathcal{N}(0,C_{Y})$ where $C_{Y}$ is the experimental covariance
matrix. In Eq. 92, each basis vector corresponds to a separate data point, and
the vector of shifts $\eta$ permits correlations between data points according
to the covariance matrix provided by the experiments. Given the data, the
NNPDF approach is to compute a MAP estimator along the lines discussed in the
previous section, i.e. finding the model that minimises the $\chi^{2}$ to the
data. The key difference between the NNPDF approach and the classical MAP
estimator is that instead of fitting the observational data given by Eq. 92,
an ensemble of model replicas are fitted each to an independently sampled
instance of pseudo-data, which is generated by augmenting $y_{0}$ with some
noise, $\epsilon^{(k)}$,

$\mu^{(k)}=y_{0}+\epsilon^{(k)}=f+\eta+\epsilon^{(k)},$ (93)

where $k$ is the replica index and each instance of the noise, $\epsilon$, is
drawn independently from the same Gaussian from which the observational noise
is drawn from, i.e. $\epsilon\sim\mathcal{N}(0,C_{Y})$. For each replica $k$,
$\mu^{(k)}$ is a vector in ${}^{N_{\rm data}}$. The distribution of pseudo-
data in a simple one-dimensional example is shown in Fig. 2. Note that, if we
were to repeat this construction multiple times, the true value $f$ would be
within a 1$\sigma$ interval centred at $y_{0}$ with a 68% probability.

Figure 2: Histogram showing the distribution of $10^{4}$ replicas generated
around an experimental value $y_{0}$ with unit variance. The central value
$y_{0}$, which is represented by the solid dot at the centre of the replica
distribution, is drawn from a Gaussian distribution with unit variance centred
at the true value $f$, which is assumed to be the origin in this plot. The
central value $y_{0}$ is kept fixed for the whole set of replicas.

The parameters for each model replica maximise the likelihood evaluated on the
corresponding pseudo-data. We can think of this approach as a special case of
MAP estimation, as described in Eq. 23, where there is no model prior that
regulates the likelihood. Another way of viewing this is to take
$C_{X}^{-1}\to 0$ in Eq. 23, as was done to obtain the result in Eq. 49.
Either way, there is no prior information about the model. The
parameterisation of the model is fixed, so the model space is the space of
parameters $u\in^{N_{\rm model}}$. In ${}^{N_{\rm model}}$, we find the
parameters which minimise the $\chi^{2}$ between the predictions from the
model and the corresponding pseudo-data $\mu^{(k)}$

$\begin{split}u_{*}^{(k)}&=\arg\min_{u^{(k)}}{\chi^{2}}^{(k)}\\\
&=\arg\min_{u^{(k)}}\sum_{ij}\left(\mathcal{G}(u^{(k)})-\mu^{(k)}\right)^{T}C_{Y}^{-1}\left(\mathcal{G}(u^{(k)})-\mu^{(k)}\right)\,,\end{split}$
(94)

where, as usual, minimising the $\chi^{2}$ is equivalent to maximising the
likelihood, $\mathcal{L}$, since $\chi^{2}\equiv-\log{\mathcal{L}}$.

As a final note: since we do not include the model prior, overall
normalisations can be omitted in Eq. 94. It is clear however that if we were
including a model prior in our MAP, it is important that the relative
normalisation between the likelihood function and the model prior is clearly
specified.

###  3.2 Fluctuations of fitted values

It is not immediately obvious that the MC methodology, maximising the
likelihood on an ensemble of pseudo-data replicas, guarantees that the model
replicas are indeed sampled from the posterior distribution of parameters
given data as described e.g. in Eq. 21. In order to investigate this issue, we
will again consider a model whose predictions are linear in the model
parameters, so that the posterior distribution of model parameters can be
computed explicitly. DIS observables can be seen as an example of a linear
model. Let us assume that the PDF are parametrized by their values at selected
values of $x_{i}$, so that $u$ is a finite-dimensional vector

$u_{i}=u(x_{i})\,,\quad i=1,\ldots,N_{\rm model}\,.$ (95)

The observables are then computed by taking a matrix-vector multiplication of
the vector $u$ by a $N_{\rm data}\times N_{\rm model}$ matrix, $F_{ij}$, which
is called FK table in the NNPDF jargon,

$y_{i}=\sum_{j=1}^{N_{\rm model}}F_{ij}u_{j}\,,$ (96)

where $i=1,\ldots,N_{\rm data}$. In this case the forward map coincides
exactly with the FK table, $\mathcal{G}_{ij}=F_{ij}$ and the sum in Eq. 96
approximates the convolution in Eq. 2. However it should be clear that the
discussion below applies to any model whose forward map can be approximated as
Eq. 96, like a linear approximation of neural networks [6]. Let us consider
for instance the actual NNPDF parametrization, where the value of the PDF is
given by a neural network whose state is determined by a set of weights
$\theta$. In this case $X$ is the space of weights, the model will be
specified by a vector $\theta$ in the space of weights, and we are going to
continue using $u$ to denote the PDF, so that

$u_{i}=u(x_{i};\theta)\,.$ (97)

If the parameters have a sufficiently narrow distribution around some value
$\bar{\theta}$, then we can expand the expression for the observables:

$\displaystyle y_{i}$
$\displaystyle=\sum_{j}F_{ij}\left[u(x_{j};\bar{\theta})+\sum_{k}\left.\frac{\partial
u}{\partial\theta_{k}}\right|_{x_{j},\bar{\theta}}\left(\theta-\bar{\theta}\right)_{k}\right]$
(98)

and therefore

$\displaystyle y_{i}-\bar{y}_{i}$
$\displaystyle=\sum_{k}\mathcal{G}_{ik}\Delta\theta_{k}\,,$ (99)

where

$\displaystyle\mathcal{G}_{ik}$
$\displaystyle=\sum_{j}F_{ij}\left.\frac{\partial
u}{\partial\theta_{k}}\right|_{x_{j},\bar{\theta}}\,,$ (100)
$\displaystyle\bar{y}_{i}$
$\displaystyle=\sum_{j}F_{ij}u(x_{j};\bar{\theta})\,,$ (101)
$\displaystyle\Delta\theta_{k}$ $\displaystyle=\theta_{k}-\bar{\theta}_{k}\,.$
(102)

In order to get an exact analytical solution for the linear model, we
additionally require $\mathcal{G}$ to have linearly independent rows, and
therefore $\mathcal{G}C_{Y}\mathcal{G}^{T}$ is invertible. As discussed in
Sec. 2, in the absence of prior information on the model, the posterior
distribution of model parameters is a Gaussian with mean and covariance given
by Eqs. 49 and 50.

Let us now stick to the parametrization in Eq. 95 and let us deploy the NNPDF
Monte Carlo method to fitting model replicas, then in the case under study
$\arg\min_{u^{(k)}}{\chi^{2}}^{(k)}$ is found analytically by imposing that
the derivative of ${\chi^{2}}^{(k)}$ with respect to the model parameters is
zero, i.e.

$\begin{split}u_{*}^{(k)}&=(\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G})^{-1}\left(\mathcal{G}^{T}C_{Y}^{-1}y_{0}+\mathcal{G}^{T}C_{Y}^{-1}\epsilon^{(k)}\right)\,.\end{split}$
(103)

Eq. 103 shows that $u_{*}$ is a linear combination of the Gaussian variables
$\epsilon$, and therefore is also a Gaussian variable. Its probability density
is then completely specified by the mean and covariance, which can be
calculated explicitly, given that the probability density for $\epsilon$ is
known:

$\displaystyle\mathbf{E}_{\\{u_{*}\\}}\left[u_{*}\right]$
$\displaystyle=\tilde{u}=(\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G})^{-1}\mathcal{G}^{T}C_{Y}^{-1}y_{0}\,,$
(104) $\displaystyle{\rm cov}(u_{*})$
$\displaystyle=\tilde{C}_{X}=(\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G})^{-1}\,.$
(105)

These two equations show that, under the assumptions specified above,
$u_{*}\sim\mathcal{N}(\tilde{u},\,\tilde{C}_{X})$. In other words, when the
model predictions are linear in the model parameters, the NNPDF MC method is
shown to produce a sample of models exactly distributed according to the
expected posterior distribution of model parameters given the data. When we
fit PDFs, parameterised as deep fully connected neural networks, to data which
includes hadronic observables, it is clear that the forward map is non-linear,
and therefore this proof does not strictly apply. As discussed above, even for
non-linear models we can make a linear approximation of the forward map
provided that we are expanding around the MAP estimator. This means the NNPDF
MC methodology should reproduce the posterior distribution of the model given
the data, at least close to $\tilde{u}$, the central value of the fitted
replicas. Furthermore, by fluctuating the data and fitting the replicas, the
fluctuations in data space are propagated to model space non-linearly. So even
for non-linear problems, the NNPDF MC methodology can be expected to produce a
sample of models which are at least approximately distributed according to the
posterior model distribution. It remains to be shown, however, that further
away from the MAP estimator the approximation holds despite the non-linear
dependence of the model replicas on the data uncertainities.

###  3.3 Closure test

The concept of the closure test, which was first introduced in Ref. [7], is to
construct artificial data by using a known pre-existing function to generate
the true observable values, $f$. One way of achieving this is by choosing a
model $w$ and then compute $f=\mathcal{G}(w)$. Then the experimental central
values are artificially generated according to Eq. 92, where the observational
noise is pseudo-randomly generated from the assumed distribution. In Ref. [7],
$f$ is referred to as level 0 (L0) data and $y_{0}$ is referred to as level 1
(L1) data. Finally, if we use the NNPDF MC method to fit artificially
generated closure data, the pseudo-data replicas that are fitted by the model
replicas are referred to as level 2 (L2) data.

The outcome of the fit is then compared with the known true value to check the
precision and accuracy of the fitting procedure. Finding quantitative
estimators that allow us to characterise the quality of the closure fit is one
of the main problems that need to be addressed. We will discuss a new class of
estimators in the next section.

Note that in a closure test, the assumed prior of the data is fully consistent
with the particular instance of observed central values, $y_{0}$: by
construction, there are no inconsistent data sets. In the original closure
test of Ref. [7] there was also no modelisation uncertainty, the true
observable values were assumed to be obtained by applying the forward map
$\mathcal{G}$ to a vector in model space $w$. It is worth noting that the
assumption of zero modelisation uncertainties is quite strong and likely
unjustified in many areas of physics. In the context of fitting parton
distribution functions there are potentially missing higher order
uncertainties (MHOUs) from using fixed order perturbative calculations as part
of the forward map. MHOUs have been included in parton distribution fits [3,
4] and in the future these should be included in the closure test, however
this is beyond the scope of the study presented here, since MHOUs are still
not included in the NNPDF methodology by default. In the results presented in
the rest of this paper we do include nuclear and deuteron uncertainties, as
presented in [8, 9, 10], since they are to be included in NNPDF fits by
default. Extensive details for including theoretical uncertainties, modelled
as theoretical covariance matrices can be found in those references. For the
purpose of this study the modelisation uncertainty is absorbed into the prior
of the data, since

$y_{0}=\mathcal{G}(u)+\eta+\delta$ (106)

where $\delta\sim\mathcal{N}(0,C^{\rm theory})$. As long as the modelisation
uncertainty is independent of the data uncertainty, we can absorb $\delta$
into $\eta$ by modifying the data prior: $\eta\sim\mathcal{N}(0,C+C^{\rm
theory})$. In doing that, we must also update the likelihood of the data given
the model to use the total covariance $(C+C^{\rm theory})$. From now onwards
we will omit $C^{\rm theory}$ because it is implicit that we always sample and
fit data using the total covariance matrix which includes any modelisation
uncertainty we currently take into account as part of our methodology.

Mapping the closure test procedure to the quantities used in the Bayesian
treatment presented in the previous section will allow us to derive a number
of analytical results in Sec. 4.4.

##  4 Data space estimators

In order to perform a quantitative analysis of the results obtained in the
closure tests, we discuss several estimators, which are computed from the
outcome of the closure test fits. These results depend on the pseudo-data that
have been generated and therefore are stochastic variables which can
fluctuate. The values of the estimators on a single replica will not tell
anything about the quality of our fits: we need to understand their
probability distributions in order to validate our fitting procedure. We begin
this section by defining estimators in data space, i.e. estimators that are
computed from the model predictions for a set of experimental data points.
Having defined the estimators, we define criteria to characterise faithful
uncertainties. We conclude this section with a discussion of the predictions
that can be obtained for these estimators in the case of a linear model, where
analytical calculations can be performed. As we already saw in the previous
section, the analytical results cannot be applied directly to the NNPDF fits,
but they are useful examples that illustrate the expected behaviour of these
quantities.

###  4.1 Deriving the data space estimators

For a given model $u_{*}^{(k)}$, obtained from fitting the $k$-th replica, we
start by defining the model error as the $\chi^{2}$ between the model
predictions and some data central values ${{y_{0}}^{\prime}}$, normalised by
the number of data points

$\frac{1}{N_{\rm
data}}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{y_{0}}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{y_{0}}^{\prime}}\right)\,,$
(107)

and we purposely denoted the data which the model error is evaluated on as
${{y_{0}}^{\prime}}$, as opposed to the training data $y_{0}$, which is used
to determine the model parameters. The corresponding covariance is denoted
${{C_{Y}}^{\prime}}$. Note that in Eq. 107, ${{y_{0}}^{\prime}}$ is a
stochastic variable, but also $u_{*}^{(k)}$ is a stochastic variable, with its
pattern of fluctuations, since the fitted model depends on the data
$\mu^{(k)}$ that enter the fit. We define the model error $\mathcal{E}^{\rm
out}$ on the set of data ${{y_{0}}^{\prime}}$ by taking the average over the
models,

$\mathcal{E}^{\rm out}=\frac{1}{N_{\rm
data}}\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{y_{0}}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{y_{0}}^{\prime}}\right)\right]\,,$
(108)

where we defined the expectation value over the ensemble of model replicas as

$\mathbf{E}_{\\{u_{*}\\}}\left[x\right]\equiv\frac{1}{N_{\rm
replicas}}\sum_{k=1}^{N_{\rm replicas}}x^{(k)}\,.$ (109)

We could of course set ${{y_{0}}^{\prime}}=y_{0}$ and evaluate the model
performance on the fitted data however, as is common in machine learning
literature, we intend to use a separate set of test data. Ideally we would
choose ${{y_{0}}^{\prime}}$ such that ${{y_{0}}^{\prime}}$ and $y_{0}$ are
statistically independent, as in Eq. 56. This is achieved by choosing the
split such that the experimental covariance matrix is block diagonal:

$C_{X}^{\rm total}=\begin{bmatrix}C_{X}&0\\\ 0&{{C_{X}}^{\prime}}\\\
\end{bmatrix}\,.$ (110)

It is useful to perform a decomposition of Eq. 108, following usual
manipulations of the likelihood function associated with least-squares
regression in [11]. Least-squares regression is a special case of minimum
likelihood estimation, where the uncertainty on each data point is equal in
magnitude and uncorrelated. Here we review the decomposition in the more
general framework of data whose uncertainty is multigaussian. Starting with
Eq. 108 (evaluated on the ideal test data), we can complete the square

$\begin{split}&\mathcal{E}^{\rm out}=\frac{1}{N_{\rm
data}}\biggl{\\{}\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{f}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{f}^{\prime}}\right)\right]+\\\
&+\mathbf{E}_{\\{u_{*}\\}}\left[\left({{f}^{\prime}}-{{y_{0}}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{f}^{\prime}}-{{y_{0}}^{\prime}}\right)\right]+\\\
&+2\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{f}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{f}^{\prime}}-{{y_{0}}^{\prime}}\right)\right]\biggr{\\}}\,.\end{split}$
(111)

Let us now discuss these terms one by one, starting from the last two. The
second term is the shift associated with evaluating the model error on noisey
test data and the final term is a cross term which we will deal with later. We
therefore focus next on further decomposing the first term,

$\begin{split}&\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{f}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-{{f}^{\prime}}\right)\right]=\\\
&=\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)\right]+\\\
&+\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)\,,\end{split}$
(112)

where we have used the fact that the second term is constant across replicas
and the cross term that arises in this decomposition is zero when the
expectation value across replicas is taken. The first term in this expression
we call the variance and the second term is the bias.

As previously mentioned $\mathcal{E}^{\rm out}$ should be considered a
stochastic estimator, in theory we could take the expectation value across
training data $y_{0}$ and test data ${{y_{0}}^{\prime}}$, the latter of which
cancels the cross term in Eq. LABEL:eq:EoutDecomposition. The final result of
that would be

$\mathbf{E}_{y_{0},{{y_{0}}^{\prime}}}[\mathcal{E}^{\rm
out}]=\mathbf{E}_{y_{0}}[{\rm bias}]+\mathbf{E}_{y_{0}}[{\rm
variance}]+\mathbf{E}_{{{y_{0}}^{\prime}}}[{\rm noise}]\,.$ (113)

We are not interested in the observational noise term, since it is independent
of the model and in the limit of infinite test data
$\mathbf{E}_{{{y_{0}}^{\prime}}}[{\rm noise}]\to 1$. The two estimators of
interest are independent of the test data, and therefore we only need to take
the expectation value over the training data. In practical terms, this can be
achieved by running multiple closure fits, each with a different observational
noise vector $\eta$, and taking the average i.e.

$\mathbf{E}_{y_{0}}[x]=\frac{1}{N_{\rm fits}}\sum_{j=1}^{N_{\rm fits}}x.$
(114)

Clearly this is resource intensive, and requires us to perform many fits. In
NNPDF3.0 [7], single replica proxy fits were used to perform a study of the
uncertainties. Here we have expanded the data-space estimators used in the
closure fits and also will be using multiple full replica fits to calculate
various expectation values - made possible by our next generation fitting
code.

###  4.2 Geometric Interpretation

It is possible to interpret the relevant data space estimators geometrically,
by considering a coordinate system where each basis vector corresponds to an
eigenvector of the experimental covariance matrix normalised by the square
root of the corresponding eigenvalue. An example of this is given in Fig. 3,
where for simplicity we have considered a system with just two data points,
i.e. a two-dimensional data space, with a diagonal covariance. The origin of
the coordinate system is the true value of the observable. The observational
noise in these coordinates corresponds to a unit circle centred in the origin
as shown in Fig. 3. If the experimental covariance is faithful, there is a 68%
probability that the experimental value $y_{0}$ is within this unit circle.
Fig. 3 shows one possible instance of $y_{0}$. Repeating the entire fit
procedure multiple times requires generating new sets of experimental data
$y_{0}$. The average over $y_{0}$ mentioned above, is precisely the average
over multiple fits, restarting the procedure each time from a new instance of
$y_{0}$.

For a given $y_{0}$ the replicas are generated as a set of points Gaussianly
distributed around it and therefore, in the limit of a large number of
replicas, 68% of them will fall within a unit circle centred in $y_{0}$. This
is the dashed circe in the figure. Clearly there is also a 68% probability
that the true value (i.e. the origin in our plot) is inside this second
circle. The model predictions, one for each replica, are then a set of points,
whose mean is $E_{\epsilon}[g]$. The mean squared radius of those points is
what we call the variance. The bias is the l2-norm of the vector between the
origin and the mean of the model predictions.

A faithful representation of the errors requires that the true value, i.e. the
origin of the coordinate system in our figure, has 68% probability of being
within 1$\sigma$ from the central value of the fit, which is given by
$E_{\epsilon}[g]$. Looking at the figure again, the probability for the origin
to be inside the shaded circle must be 68%. We will discuss faithful errors in
more detail in the next subsection.

$1/\sigma_{1}$$1/\sigma_{2}$$y_{0}$$E_{\epsilon}[g]$ Figure 3: Example of
geometric interpretation of closure test estimators. The origin is the true
observable values for each data point. The level one data (or experimental
central values) are shifted away from this by $\eta$. In this example the
covariance matrix is diagonal, so the eigenvectors correspond to the two data
points, the square root of the eigenvalues are simply the standard deviation
of those points. This is without loss of generality because any multivariate
distribution can be rotated into a basis which diagonalises the covariance
matrix. The 1-sigma observational noise confidence interval is a unit circle
centered on the origin. Some closure estimators can be understood as l2-norms
of the vectors connecting points, i.e the bias is the l2-norm of the vector
from the origin to the central value of the predictions.

###  4.3 Faithful uncertainties in data space

The two closure estimators of interest, bias and variance, can be used to
understand faithful uncertainties in a practical sense. If we return to Eq.
113 we can examine both estimators in a more detail.

#### Variance

The variance in the above decomposition refers to the variance of the model
predictions in units of the covariance

$\begin{split}{\rm variance}&=\frac{1}{N_{\rm
data}}\mathbf{E}_{\\{u_{*}\\}}\Big{[}\\\
&\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)\Big{]}\,,\end{split}$
(115)

which can be interpreted as the model uncertainty in the space of the test
data. It is instructive to rephrase Eq. 115 as

${\rm variance}=\frac{1}{N_{\rm data}}{\rm Tr}\left[{{C}^{\prime}}^{(\rm
replica)}{{C_{Y}}^{\prime}}^{-1}\right],$ (116)

where

${{C}^{\prime}}^{(\rm
replica)}=\mathbf{E}_{\\{u_{*}\\}}\left[\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)\left({{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]\right)^{T}\right]$
(117)

is the covariance matrix of the predictions from the model replicas. Note that
we can rotate to a basis where ${{C_{Y}}^{\prime}}$ is diagonal,

$\left({{C_{Y}}^{\prime}}^{-1}\right)_{ij}=\frac{1}{\left({{\sigma}^{\prime}}_{i}\right)^{2}}\delta_{ij}\,,$
(118)

then we can rewrite Eq. 116 as

${\rm variance}=\frac{1}{N_{\rm data}}\,\sum_{i}\frac{{{C}^{\prime}}^{(\rm
replica)}_{ii}}{\left({{\sigma}^{\prime}}_{i}\right)^{2}}\,.$ (119)

The numerator in the right-hand side of the equation above is the variance of
the theoretical prediction obtained from the fitted replicas, while the
denominator is the experimental variance, the average is now taken over
eigenvectors of the experimental covariance matrix. Note that this ratio does
not need to be equal to one, the model error on a given point can be smaller
than the experimental one, since all theoretical estimates are correlated by
the underlying theoretical law.

#### Bias

Similarly, the bias is defined as the difference between the expectation value
of the model predictions and the true observable values in units of the
covariance, i.e.

${\rm bias}=\frac{1}{N_{\rm
data}}\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)^{T}{{C_{Y}}^{\prime}}^{-1}\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)\,.$
(120)

The smaller the bias, the closer the central value of the predictions is to
the underlying law. In Eq. 113, the expectation value is taken across the
prior distribution of the training data, which yields

$\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{1}{N_{\rm data}}{\rm
Tr}\left[{{C}^{\prime}}^{(\rm central)}{{C_{Y}}^{\prime}}^{-1}\right]\,,$
(121)

where we have introduced ${{C}^{\prime}}^{(\rm central)}$ as the covariance of
the expectation value of the model predictions,

${{C}^{\prime}}^{(\rm
central)}=\mathbf{E}_{y_{0}}\left[\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}\right)^{T}\right]\,.$
(122)

The point is that the bias on the test data is a stochastic variable which
depends on the central value of the training data through $u_{*}^{(k)}$. The
matrix ${{C}^{\prime}}^{(\rm central)}$ describes the fluctuations of the
central value of the model prediction around the true observable values as we
scan different realisations of the training data.

It is important to stress the difference between variance and bias. In the
case of the variance, we are looking at the fluctuations of the replicas
around their central value for fixed $y_{0}$. This is related to the ensemble
of model replicas we provide as the end product of a fit and can be calculated
when we have one instance of $y_{0}$, provided by the experiments. In the case
of the bias we consider the fluctuations of the central value over replicas
around the true theoretical prediction as the values of $y_{0}$ fluctuate
around $f$. This latter procedure is only possible in a closure test, where
the underlying true observable is known. The bias as defined here yields an
estimate of the fluctuations of the MAP estimator if we could do multiple
independent instances of their measurements from each experiment.

#### Bias-variance ratio

Finally, the bias-variance ratio is defined as

$\mathcal{R}_{bv}\equiv\sqrt{\frac{\mathbf{E}_{y_{0}}[{\rm
bias}]}{\mathbf{E}_{y_{0}}[{\rm variance}]}}\,,$ (123)

where we have taken the square root, since bias and variance are both mean
squared quantities. The value of $\mathcal{R}_{bv}$ yields a measurement of
how much uncertainties are over or under estimated. If the uncertainties are
completely faithful, then $\mathcal{R}_{bv}=1$. We note that
$\mathcal{R}_{bv}$ is not completely general: it is not a measure defined in
model space and depends on the choice of test data. Therefore it only gives
local information on the model uncertainties. If the distribution of the
expectation value of model predictions is gaussian centered on the true
observable values, with covariance ${{C}^{\prime}}^{(\rm central)}$ and the
distribution of the model replicas is also gaussian, with covariance
${{C}^{\prime}}^{(\rm replica)}$ then model uncertainties are faithful if

${{C}^{\prime}}^{(\rm central)}{{{C}^{\prime}}^{(\rm replica)}}^{-1}=1.$ (124)

The difficulty with calculating Eq. 124 comes from the fact that
${{C}^{\prime}}^{(\rm replica)}$ is likely to have large correlations which
would lead it to be singular or ill-conditioned. As a result, any error
estimating ${{C}^{\prime}}^{(\rm replica)}$ from finite number of replicas
could lead to unstable results. $\mathcal{R}_{bv}$ overcomes this instability
by taking the ratio of the average across test data of these matrices, in
units of the experimental covariance matrix. There may still be large relative
errors for smaller eigenvalues of ${{C}^{\prime}}^{(\rm replica)}$, but these
should not lead to instabilities in $\mathcal{R}_{bv}$ unless they correspond
to directions with very low experimental uncertainty. As an extra precaution,
we shall estimate an uncertainty on $\mathcal{R}_{bv}$ by performing a
bootstrap sample on fits and replicas.

#### Quantile statistics

When the closure test was first presented in [7], there was an estimator
introduced in the space of PDFs which also aimed to estimate faithfulness of
PDF uncertainties using the combined assumption of Gaussian PDF uncertainties
and quantile statistics, called $\xi_{1\sigma}$. Here we can define an
analogous expression in the space of data,

$\xi^{(\rm data)}_{n{{\sigma}^{\prime}}}=\frac{1}{N_{\rm
data}}\sum_{i}^{N_{\rm data}}\frac{1}{N_{\rm fits}}\sum_{l}^{N_{\rm
fits}}I_{[-n{{\sigma}^{\prime}}_{i}^{(l)},n{{\sigma}^{\prime}}_{i}^{(l)}]}\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}_{i}\right]^{(l)}-{{f}^{\prime}}_{i}\right),$
(125)

where ${{\sigma}^{\prime}}_{i}^{(l)}=\sqrt{{{C}^{\prime}}^{(\rm
replica)}_{ii}}$ is the standard deviation of the theory predictions estimated
from the replicas of fit $l$ and $I_{[a,b]}(x)$ is the indicator function,
which is one when $a\leq x\leq b$ and zero otherwise. In other words,
$\xi^{(\rm data)}_{n{{\sigma}^{\prime}}}$ is counting how often the difference
between the prediction from the MAP estimator and the true observable value is
within the $n\sigma$-confidence interval of the replicas, assuming they’re
Gaussian. Since ${{C}^{\prime}}^{(\rm replica)}$ is primarily driven by the
replica fluctuations, we assume that it is roughly constant across fits, or
independent upon the specific instance of observational noise. This allows us
to write $\xi^{(\rm data)}_{n{{\sigma}^{\prime}}}$ for a specific data point
in the limit of infinite fits, each to a different instance of the data as

$\xi^{(\rm
data)}_{n{{{\sigma}^{\prime}}_{i}}}=\int_{-\infty}^{\infty}I_{[-n{{\sigma}^{\prime}}_{i},n{{\sigma}^{\prime}}_{i}]}\,\left(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}_{i}\right]^{(l)}-{{f}^{\prime}}_{i}\right)\rho(\eta)\,{\rm
d}(\eta)\,,$ (126)

where $\mathbf{E}_{\\{u_{*}\\}}\left[\mathcal{G}_{i}\right]^{(l)}$ has
implicit conditional dependence on $\eta$. If the distribution of
$\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}_{i}\right]^{(l)}-{{f}^{\prime}}_{i}$
is Gaussian, centered on zero, we can define
${{{\hat{\sigma}}^{\prime}}}_{i}=\sqrt{{{C}^{\prime}}^{(\rm central)}_{ii}}$.
In this case

$\xi^{(\rm data)}_{n{{{\sigma}^{\prime}}_{i}}}={\rm
erf}\left(\frac{n{{\sigma}^{\prime}}_{i}}{{{\hat{\sigma}}^{\prime}}_{i}\sqrt{2}}\right),$
(127)

which is simply the standard result of integrating a gaussian over some finite
symmetric interval.

The analogy between $\mathcal{R}_{bv}$ and $\xi^{(\rm
data)}_{n{{\sigma}^{\prime}}}$ is clear, the ratios of uncertainties are both
attempts to quantify Eq. 124 whilst keeping effects due to using finite
statistics under control. Whilst with $\mathcal{R}_{bv}$ we take the average
over test data before taking the ratio, $\xi^{(\rm
data)}_{n{{\sigma}^{\prime}}}$ instead takes the ratio of the diagonal
elements - ignoring correlations. Since the predictions from the model will be
compared with experimental central values, taking into account experimental
error, we find it more natural to calculate $\xi^{(\rm
data)}_{n{{\sigma}^{\prime}}}$ in the basis which diagonalises the
experimental covariance of the test data as in Eq. 118. If we assume that in
this new basis, that both $\frac{{{C}^{\prime}}^{(\rm
replica)}_{ii}}{\left(\sigma^{\prime}_{i}\right)^{2}}$ and
$\frac{{{C}^{\prime}}^{(\rm
central)}_{ii}}{\left(\sigma^{\prime}_{i}\right)^{2}}$ are approximately
constant for all eigenvectors of the experimental covariance matrix, then we
recover the approximation

$\xi^{(\rm data)}_{n{{\sigma}^{\prime}}}\sim{\rm
erf}\left(\frac{n\mathcal{R}_{bv}}{\sqrt{2}}\right).$ (128)

Whilst it is clear that Eq. 128 is reliant on a fair few assumptions which may
not hold, we will use the comparison of $\xi^{(\rm
data)}_{n{{\sigma}^{\prime}}}$ with $\mathcal{R}_{bv}$ to consider how valid
these assumptions may be.

###  4.4 Closure estimators - Linear problems

Once again we return to the linear model framework set out in Sec. 3.2. We can
perform an analytical closure test in this framework, and check our
understanding of the closure estimators. Consider the true observable values
for the test data are given by

${{f}^{\prime}}={{\mathcal{G}}^{\prime}}w$ (129)

where $w\in X$, and we assume that the number of (non-zero) parameters in the
underlying law is less than or equal to the number of parameters in the model,
$N_{\rm law}\leq N_{\rm model}$. Using the previous results from Sec. 3.2, we
can write down the difference between the true observables and the predictions
from the MAP estimator (or the expectation of the model predictions across
model replicas - in the linear model these are the same)

$\begin{split}\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}&={{\mathcal{G}}^{\prime}}(\tilde{u}-w)\\\
&={{\mathcal{G}}^{\prime}}\tilde{C}_{X}\mathcal{G}^{T}C_{Y}^{-1}\,\eta\,,\end{split}$
(130)

where we recall that $\mathcal{G}$ is the forward map to the training
observables and $y_{0}$ are the corresponding training central values. Because
the training observables do not necessarily coincide with the data used to
compute the estimators, we have two different maps $\mathcal{G}$ and
$\mathcal{G}^{\prime}$. Calculating the covariance across training data of
$\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}}$
gives

${{C}^{\prime}}^{(\rm
central)}={{\mathcal{G}}^{\prime}}\tilde{C}_{X}{{\mathcal{G}}^{\prime}}^{T}\,,$
(131)

so the full expression for $\mathbf{E}_{y_{0}}[{\rm bias}]$ is given by

$\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{1}{N_{\rm data}}{\rm
Tr}\left[{{\mathcal{G}}^{\prime}}\tilde{C}_{X}{{\mathcal{G}}^{\prime}}^{T}{{C_{Y}}^{\prime}}^{-1}\right].$
(132)

We note that if the test data is identical the data the model was fitted on,
we recover an intuitive result $\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{N_{\rm
model}}{N_{\rm data}}$. In the case of a polynomial model, where the parameter
of the model are the coefficients of the polynomial function, the maximum
value which $N_{\rm model}$ can take whilst $\mathcal{G}$ still has linearly
independent rows is $N_{\rm data}$ and in this case the
$\mathbf{E}_{y_{0}}[{\rm bias}]$ takes its maximum value of 1. The central
predictions from the model exactly pass through each data point.

We can perform a similar exercise on the model replica predictions. The
difference between the predictions from model replica $(k)$ and the
expectation value of the model predictions is

$\begin{split}{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)-\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]&={{\mathcal{G}}^{\prime}}(u_{*}^{(k)}-\tilde{u})\\\
&={{\mathcal{G}}^{\prime}}\tilde{C}_{X}\mathcal{G}^{T}C_{Y}^{-1}\,\epsilon\,.\end{split}$
(133)

Since $\epsilon$ and $\eta$ follow the same distribution, it is clear that

${{C}^{\prime}}^{(\rm replica)}={{C}^{\prime}}^{(\rm central)},$ (134)

which simply means that

${\rm variance}=\mathbf{E}_{y_{0}}[{\rm bias}].$ (135)

We recall that when the map is linear, the NNPDF MC methodology generates
replicas which are sampled from the posterior distribution of the model given
the data. We have shown here that provided the underlying law belongs to the
model space, the posterior distribution of the model predictions satisfy the
requirement that $\mathcal{R}_{bv}=1$.

We note that due to the invariance of the trace under cyclic permutations, we
can rearrange Eq. 132 as

$\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{1}{N_{\rm data}}{\rm
Tr}\left[\tilde{C}_{X}{{\mathcal{G}}^{\prime}}^{T}{{C_{Y}}^{\prime}}^{-1}{{\mathcal{G}}^{\prime}}\right]\,,$
(136)

where the term
${{\mathcal{G}}^{\prime}}^{T}{{C_{Y}}^{\prime}}^{-1}{{\mathcal{G}}^{\prime}}$
can be understood as the covariance matrix of the posterior distribution in
model space given the test data, with zero prior knowledge of the model, which
we denote as $\tilde{C}_{X}^{\prime}$:

$\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{1}{N_{\rm data}}{\rm
Tr}\left[\tilde{C}_{X}{\tilde{C}_{X}}^{\prime-1}\right]\,,$ (137)

where we emphasise that the covariance matrices $\tilde{C}_{X}$ and
${{\tilde{C}_{X}}^{\prime}}$ are obtained from completely independent Bayesian
inferences with no prior information on the model parameters, unlike in Eq. 63
where a sequential marginalisation causes ${{\tilde{C}_{X}}^{\prime}}$ to
depend on $\tilde{C}_{X}$.

Alternatively, if we perform a sequential marginalisation of the data, and use
the result in Eq. 63, but then take ${{C_{Y}}^{\prime}}^{-1}\to 0$, i.e. there
is no information on the observables in the test set, then the covariance of
the posterior model distribution is

$\tilde{C}_{X}^{-1}=\mathcal{G}^{T}C_{Y}^{-1}\mathcal{G}\,,$ (138)

which is identical to the posterior model distribution given just the training
data - as one would expect. Now we can express bias (or variance) as

$\mathbf{E}_{y_{0}}[{\rm bias}]=\frac{1}{N_{\rm data}}{\rm
Tr}\left[{{\tilde{C}_{Y}}^{\prime}}{{C_{Y}}^{\prime}}^{-1}\right]\,,$ (139)

where ${{\tilde{C}_{Y}}^{\prime}}$ is the covariance of the posterior
distribution of ${{y}^{\prime}}$ with no prior information on that data. This
might seem peculiar because in determining ${{\tilde{C}_{Y}}^{\prime}}$ we
took the limit ${{C_{Y}}^{\prime}}^{-1}\to 0$, which encodes the fact that we
had no prior information on the unseen data, however in Eq. 139 we require
${{C_{Y}}^{\prime}}^{-1}$ to be finite. We rationalise Eq. 139 as a comparison
between the posterior distribution in the space of data of some unseen
observables to an independently determined prior from performing the relevant
experiment which measures the same observables. The two distributions can be
compared when the independently measured experimental data is published.

#### Underparameterised model

Note that if we were to choose the number of model parameters such that
$N_{\rm law}>N_{\rm model}$, then the variance would be unaffected, since the
underlying law parameters cancel. However, the bias would now contain an extra
term from the extra parameters in the underlying law, schematically:

$\begin{split}(\mathbf{E}_{\\{u_{*}\\}}\left[{{\mathcal{G}}^{\prime}}\left(u_{*}^{(k)}\right)\right]-{{f}^{\prime}})_{i}=\sum_{1\leq
j\leq N_{\rm
model}}{{\mathcal{G}}^{\prime}}_{ij}(\tilde{u}-w)_{j}-\sum_{N_{\rm
model}<j\leq N_{\rm law}}{{\mathcal{G}}^{\prime}}_{ij}w_{j},\end{split}$ (140)

which would mean that $\mathcal{R}_{bv}\neq 1$. This demonstrates that
requiring $\mathcal{R}_{bv}=1$ demands that the model space is suitably
flexible, if the underlying law is parameterised then this can be summarised
as requiring $w\in X$. Note that in the underparameterised regime the model
replicas are still drawn from the posterior distribution, however because
$w\notin X$ we have somehow invalidated the assumptions that go into the
relation between model predictions and the data-space prior.

Although $\mathcal{R}_{bv}$ was largely chosen on practical grounds, we see
that it is still a stringent test that our assumptions are correct and that
the distribution our model replicas are drawn from is meaningful, this is what
we mean when we say faithful uncertainties.

An unfortunate trade-off when using $\mathcal{R}_{bv}$ is that it cannot be
used as a diagnostic tool, and is instead used simply for validation. For
example, if $\mathcal{R}_{bv}>1$, then we cannot know whether there was a
problem with the fitting methodology used to generate the model replicas or a
deeper issue such as an inflexible model.

##  5 Numerical Setup and Results

In this section we first introduce the experimental setup used to run the
closure tests, and then discuss the actual results. Following the study
performed in Ref. [7], we first analyse the relative size of the different
components of PDFs uncertainty, comparing the changes between methodologies
used to produce the NNPDF3.1 [12] and NNPDF4.0 [1] sets of PDFs respectively.
We then move to the data space estimators $\mathcal{R}_{bv}$ and $\xi^{(\rm
data)}_{1{{\sigma}^{\prime}}}$, which have been computed only for NNPDF4.0
222As pointed out before, the computation of the expectation value over
training data, defined in Eq. 114, is made possible by the efficiency of the
new framework [13].. The results here act both as a proof of principle of new
estimators presented in this paper but also as part of a suite of
methodological validation tools, see also the ”future tests” [14], used to
understand the PDF uncertainties of the recent NNPDF4.0 set of PDFs. For the
purpose of understanding how the results here were produced, we will briefly
describe the key features of the NNPDF4.0 methodology, but refer the reader to
Ref. [1] for a full discussion on how these methodological choices were made.

###  5.1 Closure test setup

Using neural networks to fit PDFs has been discussed many times in previous
NNPDF publications, see for example [7, 12]. A new feature of NNPDF4.0 is
that, for the default fit performed in the evolution basis, a single neural
network parameterises all 8 PDF flavours
$\\{g,\Sigma,V,V_{3},V_{8},T_{3},T_{8},T_{15}\\}$ at the initial scale. The
PDF for a single flavour $j$, at the initial scale $Q_{0}=1.65~{}{\rm GeV}$ is
given by

$f_{j}(x,Q_{0})=NN(x,\ln x|u)_{j}*x^{1-\alpha_{j}}*(1-x)^{\beta_{j}},$ (141)

where $\alpha$ and $\beta$ are the preprocessing exponents, which control the
PDF behaviour at $x\to 0$ and $x\to 1$ respectively and $NN(x,\ln x|u)_{j}$ is
the $j^{\rm th}$ output of the neural network, which takes $x$ and $\ln x$ as
input. As discussed in Sec. 3.1, an ensemble of models is fitted, each one is
a MAP estimator of the corresponding pseudo-data it is fitted on. An
optimization algorithm is used to try and find the parameters which maximise
the likelihood. There are clearly many choices with respect to
hyperparameters, the discussion of how these choices have been made is beyond
the scope of this paper and left to the full NNPDF4.0 release [1]. A summary
of the hyperparameters used to produce results presented in this paper are
provided in Tab. 3.

As input to the closure test, a single replica was drawn randomly from a
previous NNPDF fit to experimental data. This has generally more structure
than the final central PDF and it is therefore a more general choice than any
final central fit. We refer to this as the underlying law and the
corresponding predictions as the true observable values. An example of the
gluon input is provided in Fig. 4. In principle any function could be used as
underlying law, however it makes sense to use a realistic input.

Figure 4: The green line is the input underlying law for the gluon PDF, which
is sampled from the ensemble of replicas from a fit to data. The central value
and the 68% confidence interval for those replicas are plotted as the orange
line and band respectively.

The observables used in the fits are a subset of the full NNPDF4.0 dataset.
For convenience, we chose to fit the PDFs on a variant of the NNPDF3.1 dataset
used in Ref. [15], which is described in detail in a study of the
determination of the strange PDF [16]. The datasets used in the calculation of
statistical estimators are the new datasets which have been included in
NNPDF4.0 and discussed in detail in the main release [1]. For a full summary
of the observables used in the test data and a visual representation of the
kinematic region of both the training and testing data, see App. B.

The partitioning of the available data into fitting and testing should not
affect the interpretation of the closure test results. However, one could
consider splitting the data into fitting and testing in a way which was
physically motivated e.g. the partitioning could have been stratified such
that the kinematic coverage of the data of each partition was approximately
equal. Alternatively, since the data is generated from the theory predictions
produced by the input underlying law, one could even produce completely
artificial data using a different set of FK tables.

In order to compute the expectation value over training data defined in Eq.
114, we generate 30 different sets of experimental central values (or L1
data), as discussed in Sec. 3.3, for the fitted 3.1-like dataset. Each set of
experimental central values was then fitted following NNPDF4.0 methodology
[1], producing 40 pseudo-data replicas.

###  5.2 Different components of the PDF uncertainty

As already discussed in Ref. [7], fitting to L0, L1 and L2 pseudo-data allows
us to validate different aspects of the fitting procedure. In an L0 fit, we
fit multiple time the exact same set of data, which corresponds to the theory
prediction from the chosen model. The fitted pseudo-data is the result
obtained by applying the forward map to the model. It is clear that in this
case the quality of the fit can be improved at will, provided that the
parametrization is flexible enough and the minimization algorithm is
efficient. There are indeed multiple solutions that reproduce exactly the data
set, while interpolating between the data points. In an L1 fit the data have
been fluctuated around the theoretical prediction – mimicking thereby the
central values of experimental measurements. The true model no longer
reproduces the data; instead it will yield a $\chi^{2}$ of order one. The
pseudo-data are held fixed, fluctuations from one replica to the next are due
to the existence of multiple solutions that hold a similar value for the
residual $\chi^{2}$ at the end of the minimization process. Finally, in the L2
fits, the fluctuations of the data are reproduced by the replicas, and
propagated to the model function when fitting the data for each replica. Since
the NNPDF fitting methodology has changed in the latest release, adopting the
procedure described in Ref. [17], it is important to compare the uncertainties
at L0, L1 and L2 that are obtained with the NNPDF4.0 methodology with the
corresponding ones obtained with the NNPDF3.1 framework.

An example of these relative errors on fitted PDFs is shown in the upper panel
of Fig. 5 where results from L0, L1 and L2 closure tests are displayed on the
same plot in the case of the gluon distribution. Each fit is normalized to the
corresponding central value. We note how the L0 and L1 uncertainties tend to
increase in the $x$ regions where less experimental data are available, namely
at small and large-$x$, where the model is left unconstrained and has more
freedom to fluctuate, while they are considerably reduced in the data region
where the contribution of the L2 error, induced by the actual experimental
data, becomes more important. The fact that the data uncertainty is not always
the dominant component of the PDF uncertainty was already stressed in Ref.
[7]. Improved methodologies should therefore aim to reduce the L0 and L1
error.

In this respect, it is interesting to compare these results with what we find
with the NNPDF3.1 methodology. The corresponding plots are shown in the lower
panel of Fig. 5: unlike the NNPDF4.0 case, here the PDF uncertainty is always
dominated by the L0 uncertainty, even in those kinematic regions where
experimental data are present. We can conclude that moving to the NNPDF4.0
methodology, thanks to the optimized hyperparameters listed in Tab. 3, we
observe a marked reduction of the L0 uncertainty, see Fig. 6, where the latter
is plotted on the same panel for both NNPDF4.0 and NNPDF3.1. The better
efficiency of the NNPDF4.0 methodology can be also appreciated by looking at
Fig. 7, where we plot the distribution across replicas of the L0 closure test
$\chi^{2}$ as a function of the training epoch and genetic algorithm
generation for the NNPDF4.0 and NNPDF3.1 methodology respectively. As
mentioned before, in a L0 closure test the quality of the fit can in principle
be improved at will, assuming the use of an efficient optimizer. For the final
$\chi^{2}$ of the central value of each fit, plotted as a black dashed line,
we find $0.002$ and $0.012$ in the NNPDF4.0 and NNPDF3.1 fitting code
respectively, showing the better efficiency of the new methodology.

It should be noted how Fig. 5 only provide a qualitative assessment of the
relative size of the different components of PDFs uncertainty. Despite being
useful to assess how the methodology has improved with respect to the previous
one, they do not provide any quantitative estimation of the faithfulness of
PDF uncertainty. This is best achieved in data space, using the new estimators
introduced in the previous sections.

Figure 5: Relative PDF error for the gluon distribution in the NNPDF4.0 (upper
panel) and NNPDF3.1 (lower panel) methodology, plotted in logarithmic (left)
and linear scale (right).

Figure 6: Level 0 uncertainty in the NNPDF4.0 and NNPDF3.1 methodology.

Figure 7: $\chi^{2}$ distribution across replicas as a function of the
training epoch and genetic algorithm generation for the NNPDF4.0 (left) and
NNPDF3.1 (right) methodology. The black dotted line in each plot represents
the final $\chi^{2}$ of the central value of the corresponding fit, equal to
$0.002$ and $0.012$ respectively.

###  5.3 Data space estimators

Results for $\mathcal{R}_{bv}$ and the corresponding uncertainty calculated on
the test data are shown in the first column of Tab. 1. The uncertainties,
which take into account the finite size of the replicas and fits samples, have
been computed by performing bootstrap, [18], where we randomly sample from
both fits and replicas and re-calculate $\mathcal{R}_{bv}$. The value and
error presented in the table is then the mean and standard deviation across
bootstrap samples. We checked that the distribution of the estimator across
bootstrap samples is indeed Gaussian. We also checked that increasing the
number of fits and replicas reduced the bootstrap error but the central values
were the same within the estimated bootstrap uncertainties. We see that
overall $\mathcal{R}_{bv}$ is consistent with 1, which gives a good indication
that, at least for the unseen data used in this study, the uncertainties are
faithful

$\mathcal{R}_{bv}$ | $\xi_{1\sigma}$ | ${\rm erf}(\mathcal{R}_{bv}/\sqrt{2})$
---|---|---
$1.03\pm 0.05$ | $0.69\pm 0.02$ | $0.67\pm 0.03$
Table 1: In the first column we show the bias-variance ratio,
$\mathcal{R}_{bv}$, for unseen data, summarised in Tab. 2. The uncertainty is
estimated by performing a bootstrap sample across fits and replicas and
calculating the standard deviation. We see that overall $\mathcal{R}_{bv}$ is
consistent with 1, within uncertainties. In the second and third columns we
compare the measured value of $\xi_{1}\sigma$ and the estimated value from
$\mathcal{R}_{bv}$. The two values are consistent, which suggests the
approximation that the ratio of uncertainties is approximately the same across
all data is not completely invalidated.

In Fig. 8 we compare qualitatively the distribution of bias across fits, to
the distribution of the difference between replica predictions and expectation
values of predictions (in units of the covariance) across different fits and
replicas. The square root ratio of the mean of these two distributions is
precisely $\mathcal{R}_{bv}$.

Figure 8: The green histogram is the distribution of the total bias across
fits, the orange histogram is the distribution of the difference between the
replica and central predictions squared, in units of the covariance across all
fits and replicas. This gives a qualitative picture of the full distribution,
in Tab. 1 we compare the square root of the mean of each distribution, getting
the quoted value for $\mathcal{R}_{bv}$.

As discussed in Sec. 4.3, one can define an analogous estimator in data space,
based upon $\xi_{n\sigma}$, which was defined on a grid of points in $x$ and
$Q^{2}$ in PDF space in [7]. There is not a one-to-one correspondence between
this and $\mathcal{R}_{bv}$, but a loose approximation using Eq. 127. In the
second and third columns of Tab. 1 we compare the estimated $\xi_{1\sigma}$
from substituting $\mathcal{R}_{bv}$ into Eq. 127 and to the measured value.
Despite the assumptions entering each of the two estimators differing, we see
good agreement between the $\xi_{1\sigma}$ estimated from $\mathcal{R}_{bv}$
and that measured directly. We find this result reassuring, since it indicates
not only that the total uncertainty averaged across all data is faithful, but
also that the uncertainty on each data point seems faithful. If the results
differed it would indicate some kind of imbalance, where some components of
the uncertainty are correctly represented by the replicas but other directions
are not. Finally we note how not only are the measured value and estimated
value from $\mathcal{R}_{bv}$ self consistent, but they are also consistent
with $0.68$, which further supports the argument that the model uncertainties
are faithful.

##  6 Summary

We have presented a theoretical framework for treating inverse problems from a
Bayesian perspective. In particular, the framework provides a more formal
description of what it means when we talk about propagating experimental
uncertainties into the space of models by focusing on the posterior
probability measure in model space given the available data. Strictly
speaking, there is no fitting required to obtain the expression for the
posterior distributions in the space of the data or the model, instead these
are obtained by marginalising the joint distribution. We note that sampling
from the posterior distribution of the model is, in general, highly non-
trivial; however we can show that at least for linear problems the NNPDF MC
methodology produces a sample of models which are distributed according to the
posterior model distribution predicted from the Bayesian formalism.
Furthermore, we provide evidence that even for non-linear models this result
at least holds as a good approximation close to the MAP estimator.

We then use this formal framework to think about some of the estimators which
we use as part of the NNPDF closure test. In particular we derive bias and
variance from decomposing an out of sample error function, which is understood
from a classic fitting point of view. The estimators are then related back to
the posterior distributions in the Bayesian framework. We note that the
estimators themselves are not perfect and suffer from only testing the model
uncertainties locally (in regions where the test data probes). Furthermore,
the estimators only give an approximate overall picture, and cannot be used to
diagnose where the problem arises if we find evidence that the model
uncertaintes are not faithful.

Given the framework set out here, future work should be undertaken to
generalise the closure estimators to model space. This would likely involve a
combination of the closure estimators presented here and the extension of the
Bayesian framework to infinite-dimensional spaces.

We give some closure results, as a proof of principle. The results presented
here are examined in more detail alongside the full NNPDF4.0 release [1] but
serve as an example of how the data space estimators can be practically
included even in a rather complex setting. The NNPDF4.0 methodology passes the
closure test according to the estimators that we analysed, providing evidence
that for unseen data the current NNPDF methodology appears to provide faithful
uncertainties. The estimators are not limited to this specific application,
and the results presented demonstrate how the data space estimators can be
incorporated into any inverse problem. As previously mentioned, a more general
set of estimators in model space would be the gold-standard and give us
confidence in our uncertainties for future observables which probe regions
which are not covered by either the training or testing data.

In the closure test framework described in this study, the prior on the data
is fully consistent with the generation of the observable central values from
the true observables and uncertainties by construction, which is likely not
the case in real world fits. Something which should be investigated is how we
can guarantee faithful uncertainties when this is no longer the case. For
example one could consider a closure test where the generated data was
inconsistent with the prior. Most observed inconsistentency between the prior
and the data central values is likely due to missing theoretical
uncertainties, but once all sources of theoretical uncertainties have been
accounted for there could still be tension in the data. The advantage of
viewing inverse problems from a Bayesian perspective is access to
methodologies which deal with inconsistent data (or unknown systematics) in a
Bayesian framework, for example in these cosmological studies [19, 20]. These
methods potentially offer a more formal approach to dealing with inconsistent
data, rather than ad hoc procedures.

The inclusion of MHOUs in the likelihood was justified from a Bayesian
perspective [4], but here we have drawn a connection between the model
replicas and the posterior distribution in model space, which explains why
theoretical uncertainties must be included in both the sampling of the
pseudodata replicas and the likelihood. Therefore we emphasise that the
framework set out here is not only useful for understanding model
uncertainties with the current methodology, but also for motivating future
methodological development from a Bayesian perspective. We defer more detailed
investigations of the potential of this approach to future work.

## Acknowledgements

We thank the members of the NNPDF collaboration for discussions and comments
on the present work. L. D. D. is supported by the U.K. Science and Technology
Facility Council (STFC) grant ST/P000630/1, T. G. is supported by NWO via a
ENW-KLEIN-2 project and M. W. by the STFC grant ST/R504737/1.

##  Appendix A Gaussian integrals

Theory errors can be included in this framework by allowing the distribution
of observables around the theory prediction to have a finite width, e.g. by
replacing the Dirac delta

$\delta(y-\mathcal{G}u)$ (142)

in Eq. 12 with a Gaussian

$\theta(y,u|\mathcal{G})\propto\exp\left[-\frac{1}{2}\left(y-\mathcal{G}u\right)^{T}C_{T}^{-1}\left(y-\mathcal{G}u\right)\right]\,.$
(143)

For the purposes of this study, we do not want to provide a realistic estimate
of theory errors. Instead we will be assuming that the errors are uncorrelated
and identical for all data points

$C_{T}=\sigma^{2}\mathds{1}\,,$ (144)

and we will be interested in the limit where $\sigma^{2}\to 0$.

###  A.1 Integrating out the data

Marginalizing with respect to y in this case yields

$\displaystyle\pi_{M}(u|y_{0},u_{0},\mathcal{G})$
$\displaystyle\propto\pi_{M}^{0}(u|u_{0})\,\int
dy\,\pi_{D}^{0}(y|y_{0})\theta(y,u|\mathcal{G})\,.$ (145)

The argument of the exponential in the integrand is a quadratic form in y,

$A=\left(y-y_{0}\right)^{T}C_{D}^{-1}\left(y-y_{0}\right)+\left(y-\mathcal{G}u\right)^{T}C_{T}^{-1}\left(y-\mathcal{G}u\right)\,.$
(146)

The integral can be easily evaluated by completing the square,

$A=\left(y-\tilde{y}\right)^{T}\tilde{C}_{D}^{-1}\left(y-\tilde{y}\right)+R_{D}\,.$
(147)

Comparing Eqs. 146 and 147 at order $y^{2}$ and y, yields

$\displaystyle\tilde{C}_{D}^{-1}$
$\displaystyle=\frac{1}{\sigma^{2}}\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)\,,$
(148) $\displaystyle\tilde{y}$
$\displaystyle=\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)^{-1}\left(\mathcal{G}u+\sigma^{2}C_{D}^{-1}y_{0}\right)\,,$
(149)

and therefore

$\displaystyle\tilde{y}^{T}\tilde{C}_{D}^{-1}\tilde{y}$
$\displaystyle=\frac{1}{\sigma^{2}}\left(\mathcal{G}u\right)^{T}\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)^{-1}\left(\mathcal{G}u\right)+y_{0}^{T}C_{D}^{-1}\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)^{-1}\left(\mathcal{G}u\right)+$
$\displaystyle\quad+\left(\mathcal{G}u\right)^{T}C_{D}^{-1}\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)^{-1}y_{0}+\sigma^{2}y_{0}^{T}C_{D}^{-1}\left(\mathds{1}+\sigma^{2}C_{D}^{-1}\right)^{-1}C_{D}^{-1}y_{0}\,.$
(150)

Note that the four terms in the equation above are ordered in increasing
powers of $\sigma^{2}$ and ultimately we will be interested in the limit
$\sigma^{2}\to 0$, which reproduces the Dirac delta in $\theta(y,u)$. Plugging
Eq. 150 in Eq. 147 and again comparing to Eq. 146, we find

$\begin{split}R_{D}&=\frac{1}{\sigma^{2}}\left(\mathcal{G}u\right)^{T}\left[\mathds{1}-\frac{1}{\mathds{1}+\sigma^{2}C_{D}^{-1}}\right]\left(\mathcal{G}u\right)-y_{0}^{T}C_{D}^{-1}\left(\mathcal{G}u\right)-\left(\mathcal{G}u\right)^{T}C_{D}^{-1}y_{0}+\\\
&\quad+y_{0}^{T}C_{D}^{-1}y_{0}+\mathcal{O}(\sigma^{2})\,,\end{split}$ (151)

Expanding for small values of $\sigma^{2}$ the terms of order $1/\sigma^{2}$
cancel; keeping only finite terms in the limit $\sigma^{2}\to 0$ we finally
obtain

$R_{D}=\left(\mathcal{G}u-y_{0}\right)^{T}C_{D}^{-1}\left(\mathcal{G}u-y_{0}\right)\,.$
(152)

This is exactly the result that we obtained earlier when

$\theta(y,u|\mathcal{G})=\delta(y-\mathcal{G}u)\,.$ (153)

It should not come as a surprise since in the limit where $\sigma^{2}\to 0$
the Gaussian distribution that we chose to describe the fluctutations of the
data around the theory predictions reduces indeed to a Dirac delta. The
posterior for the model is exactly the one we computed in Sect. 2. We do not
learn anything new from this exercise, but it is a useful warm-up for the next
example. The integral over y can now be performed easily, since it is yet
again a Gaussian integral.

###  A.2 Integrating out the model

Using the same approach as above, we now want to marginalise with respect to
the model in order to obtain the posterior distribution of the data:

$\displaystyle\pi_{D}(y|y_{0},u_{0},\mathcal{G})$
$\displaystyle\propto\pi_{D}^{0}(y|y_{0})\,\int
du\,\pi_{M}^{0}(u|u_{0})\theta(y,u|\mathcal{G})\,.$ (154)

We follow exactly the same procedure outlined above, starting from the
argument of the exponential

$A=\left(u-u_{0}\right)^{T}C_{M}^{-1}\left(u-u_{0}\right)+\left(y-\mathcal{G}u\right)^{T}C_{T}^{-1}\left(y-\mathcal{G}u\right)\,,$
(155)

we complete the square and rewrite it in the form

$A=\left(u-\tilde{u}\right)^{T}\tilde{C}_{M}^{-1}\left(u-\tilde{u}\right)+R_{M}\,.$
(156)

It can be readily checked that in this case

$\displaystyle\tilde{C}_{M}^{-1}$
$\displaystyle=\frac{1}{\sigma^{2}}\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{D}^{-1}\right)\,,$
(157) $\displaystyle\tilde{u}$
$\displaystyle=\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)\left(\mathcal{G}^{T}y+\sigma^{2}C_{M}^{-1}u_{0}\right)\,.$
(158)

In order to evaluate $R_{M}$, we need

$\begin{split}\tilde{u}^{T}\tilde{C}_{M}^{-1}\tilde{u}&=\frac{1}{\sigma^{2}}y^{T}\mathcal{G}\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)^{-1}\mathcal{G}^{T}y+u_{0}^{T}C_{M}^{-1}\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)^{-1}\mathcal{G}^{T}y+\\\
&\quad+y^{T}\mathcal{G}\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)^{-1}C_{M}^{-1}u_{0}+\mathcal{O}(\sigma^{2})\,.\end{split}$
(159)

Noting that

$\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)^{-1}=\frac{1}{\sigma^{2}}C_{M}-\frac{1}{\sigma^{2}}C_{M}\mathcal{G}^{T}\left(\mathcal{G}\frac{1}{\sigma^{2}}C_{M}\mathcal{G}^{T}+\mathds{1}\right)^{-1}\mathcal{G}\frac{1}{\sigma^{2}}C_{M}\,,$
(160)

we have, in the limit where $\sigma^{2}\to 0$

$\mathcal{G}\left(\mathcal{G}^{T}\mathcal{G}+\sigma^{2}C_{M}^{-1}\right)^{-1}\mathcal{G}^{T}=\mathds{1}-\sigma^{2}\left(\mathcal{G}C_{M}\mathcal{G}^{T}\right)^{-1}+\mathcal{O}(\sigma^{4})\,.$
(161)

Collecting all terms we find

$R_{M}=\left(y-\mathcal{G}u_{0}\right)^{T}\left(\mathcal{G}C_{M}\mathcal{G}^{T}\right)^{-1}\left(y-\mathcal{G}u_{0}\right)\,.$
(162)

Performing the Gaussian integral over $u$ in Eq. 154, we obtain the posterior
distribution of the data

$\pi_{D}^{y}(y)\propto\exp\left[-\frac{1}{2}\left(y-y_{0}\right)^{T}C_{D}^{-1}\left(y-y_{0}\right)-\frac{1}{2}\left(y-\mathcal{G}u_{0}\right)^{T}\left(\mathcal{G}C_{M}\mathcal{G}^{T}\right)^{-1}\left(y-\mathcal{G}u_{0}\right)\right]\,.$
(163)

As in the case above, we note that this is a Gaussian distribution,

$\pi_{D}^{y}(y)\propto\exp\left[-\frac{1}{2}\left(y-\tilde{y}\right)^{T}\tilde{C}_{D}^{-1}\left(y-\tilde{y}\right)\right]\,,$
(164)

where the mean and the covariance are given in Eqs. 54 and 55. We can rewrite
those expressions as

$\displaystyle\tilde{y}$ $\displaystyle=\mathcal{G}\tilde{u}\,,$ (165)
$\displaystyle\tilde{C}_{D}^{-1}$
$\displaystyle=C_{D}^{-1}+\left(\mathcal{G}C_{M}\mathcal{G}^{T}\right)^{-1}\,.$
(166)

In order to simplify the notation we introduce

$\hat{C}_{M}=\left(\mathcal{G}C_{M}\mathcal{G}^{T}\right)\,,$ (167)

and then

$\tilde{C}_{D}=\hat{C}_{M}-\hat{C}_{M}\left(\hat{C}_{M}+C_{D}\right)^{-1}\hat{C}_{M}\,.$
(168)

##  Appendix B Closure test setup details

###  B.1 Data

The full list of datasets included in the test set are shown in Tab. 2. The
central values are not actually used in the closure test, however we use the
experimental uncertainities in the calculation of both $\mathcal{R}_{bv}$ and
$\xi^{(\rm data)}_{1{{\sigma}^{\prime}}}$. The corresponding predictions
generated from the underlying law are used as the true observable values.
Neither of the data-space closure estimators rely on the central values of the
test datasets.

Data set | Ref.
---|---
DY E906 $\sigma^{d}_{\rm DY}/\sigma^{p}_{\rm DY}$ (SeaQuest) | [21]
ATLAS $W,Z$ 7 TeV ($\mathcal{L}=4.6$ fb-1) | [22]
ATLAS DY 2D 8 TeV | [23]
ATLAS high-mass DY 2D 8 TeV | [24]
ATLAS $\sigma_{W,Z}$ 13 TeV | [25]
ATLAS $W^{+}$+jet 8 TeV | [26]
ATLAS $\sigma_{tt}^{\rm tot}$ 13 TeV ($\mathcal{L}=$139\text{\,}{\mathrm{fb}}^{-1}$$) | [27]
ATLAS $t\bar{t}$ lepton+jets 8 TeV | [28]
ATLAS $t\bar{t}$ dilepton 8 TeV | [29]
ATLAS single-inclusive jets 8 TeV, R=0.6 | [30]
ATLAS dijets 7 TeV, R=0.6 | [31]
ATLAS direct photon production 13 TeV | [32]
ATLAS single top $R_{t}$ 7, 8, 13 TeV | [33, 34, 35]
CMS dijets 7 TeV | [36]
CMS 3D dijets 8 TeV | [37]
CMS $\sigma_{tt}^{\rm tot}$ 5 TeV | [38]
CMS $t\bar{t}$ 2D dilepton 8 TeV | [39]
CMS $t\bar{t}$ lepton+jet 13 TeV | [40]
CMS $t\bar{t}$ dilepton 13 TeV | [41]
CMS single top $\sigma_{t}+\sigma_{\bar{t}}$ 7 TeV | [42]
CMS single top $R_{t}$ 8, 13 TeV | [43, 44]
LHCb $Z\to\mu\mu,ee$ 13 TeV | [45]
Table 2:  Observables included in the test data. We wish to stress that the
observable central values themselves are not used, however the experimental
uncertainties are used in the definition of the closure estimators, and the
corresponding predictions from either the underlying law or the closure fits.

For completeness, in Fig. 9 the kinematic coverage of the training datasets,
which as mentioned is the NNPDF3.1-like dataset used in [16], and the test
datasets shown in Tab. 2 is plotted.

Figure 9: The kinematic coverage of the training and test data used to train
the models and produce results presented in this paper. The kinematics of the
two sets of data with this particular split overlaps but there are also
kinematic regions which the test dataset probes, for which there are no
training data.

###  B.2 Models

A summary of the hyperparameters for the neural networks used during the
closure fits is given in table Tab. 3. These hyperparameters were chosen as
part of an extensive hyperparameters scan, which will be explained in detail
in NNPDF4.0, for this study we simply provide the values of the
hyperparameters as a point of reference.

hyperparameter | value
---|---
architecture | 2-25-20-8
activation | tanh-tanh-linear
minimiser | NAdam
max training length | 17000 epochs
pre-processing exponents | not fitted
Table 3: Hyperparameters for neural networks used in this study. The parameter
choices, and how these choices were made are discussed in Ref. [1].

##  Appendix C Understanding NNPDF3.0 data estimators

In the closure test presented in NNPDF3.0 [7] there was a data-space estimator
which aimed to measure the level of over or under fitting,
$\Delta_{\chi^{2}}$. Here we discuss how $\Delta_{\chi^{2}}$ can emerge from
the bias-variance decomposition and then use the linear model to try and
understand it in the context of viewing the ensemble of model replicas as a
sample from the posterior distribution of the model given the data.

Despite the link between the estimators emerging from the decomposition of
$\mathcal{E}^{\rm out}$ and the posterior distribution for data which is not
used to inform the model parameters, if we perform the same decomposition as
in Sec. 4.1 but set ${{y_{0}}^{\prime}}=y_{0}$ then we find that the cross
term in the final line of Eq. LABEL:eq:EoutDecomposition does not go to zero
when the expectation across data is taken because there is a dependence on
$y_{0}$ in both the model predictions and the noisey data. As a result we have
to modify Eq. 113 to be

$\mathbf{E}_{y_{0}}[\mathcal{E}^{\rm in}]=\mathbf{E}_{y_{0}}[{\rm
bias}]+\mathbf{E}_{y_{0}}[{\rm variance}]+\mathbf{E}_{y_{0}}[{\rm
noise}]+\mathbf{E}_{y_{0}}[{\rm noise\,cross\,term}]\,,$ (169)

where we refer now to the right hand side of Eq. 169 as $\mathcal{E}^{\rm in}$
because it’s evaluated on the data used to inform the model replicas.

Now we examine the definition of $\Delta_{\chi^{2}}$ introduced in [7],
defined as the difference between the $\chi^{2}$ between the expectation value
of the model predictions and the level one data, and the $\chi^{2}$ between
the underlying observable values and the level one data. In [7] the
denominator was also set to be the second term in the numerator, however here
we slightly re-define $\Delta_{\chi^{2}}$ to instead simply be normalised by
the number of data points:

$\begin{split}\Delta_{\chi^{2}}&=\\\ \frac{1}{N_{\rm
data}}&\Big{[}\left(\mathbf{E}_{\\{u_{*}\\}}\left[\mathcal{G}\left(u_{*}^{(k)}\right)\right]-y_{0}\right)^{T}C_{Y}^{-1}\left(\mathbf{E}_{\\{u_{*}\\}}\left[\mathcal{G}\left(u_{*}^{(k)}\right)\right]-y_{0}\right)\\\
&\,-\left(f-y_{0}\right)^{T}C_{Y}^{-1}\left(f-y_{0}\right)\Big{]}\\\ &={\rm
bias}+{\rm noise\,cross\,term}\,,\end{split}$ (170)

where in the second line we show how $\Delta_{\chi^{2}}$ itself can be
decomposed to be equal to two of the terms in Eq. 169.

Constant values of $\Delta_{\chi^{2}}$ define elliptical contours in data
space centered on the level one data. $\Delta_{\chi^{2}}=0$, in particular,
defines a contour which is centered on the level one data and passes through
the underlying law. When viewing $\Delta_{\chi^{2}}$ from a classical fitting
perspective, if $\Delta_{\chi^{2}}<0$ then the expectation value of the model
predictions fit the level one data better than the underlying observables -
which indicates an overfitting of the shift, $\boldsymbol{\eta}$. Similarly,
$\Delta_{\chi^{2}}>0$ indicates some underfitting of the level one data.

If we return to the linear model we can write the analytic value of
$\Delta_{\chi^{2}}$. Firstly, since ${{y_{0}}^{\prime}}=y_{0}$ we can simplify
Eq. 132

$\begin{split}\mathbf{E}_{y_{0}}[{\rm bias}]&=\frac{1}{N_{\rm data}}{\rm
Tr}\left[\mathcal{G}\tilde{C}_{X}\mathcal{G}^{T}C_{Y}^{-1}\right]\\\
&=\frac{1}{N_{\rm data}}{\rm
Tr}\left[\tilde{C}_{X}\tilde{C}_{X}^{-1}\right]\\\ &=\frac{N_{\rm
model}}{N_{\rm data}}\,,\end{split}$ (171)

because $\tilde{C}_{X}\tilde{C}_{X}^{-1}$ is an $N_{\rm model}\times N_{\rm
model}$ identity matrix. Similarly we can write down the cross term

$\begin{split}\mathbf{E}_{y_{0}}[{\rm noise\,cross\,term}]&=\frac{-2}{N_{\rm
data}}\mathbf{E}_{y_{0}}\left[(\mathcal{G}\tilde{C}_{X}\mathcal{G}^{T}C_{Y}^{-1}\eta)^{T}C_{Y}^{-1}\eta\right]\\\
&=\frac{-2}{N_{\rm data}}{\rm
Tr}[\mathcal{G}\tilde{C}_{X}\mathcal{G}^{T}C_{Y}^{-1}]\\\ &=-2\frac{N_{\rm
model}}{N_{\rm data}}\end{split}$ (172)

which leaves us with

$\mathbf{E}_{y_{0}}[\Delta_{\chi^{2}}]=-\frac{N_{\rm model}}{N_{\rm data}}\,.$
(173)

The point is that the linear model has already been shown to be a sample from
posterior distribution of the model given the data. But from the classical
fitting point of view we would say this model has overfitted.

As such, we do not report any results with $\Delta_{\chi^{2}}$ here, because
when $\mathcal{R}_{bv}=1$, it doesn’t add much to the discussion. It may still
be useful as a diagnostic tool when $\mathcal{R}_{bv}\neq 1$, which as
discussed could be for a variety of reasons - including fitting inefficiency.
It also may be used as a performance indicator for deciding between two
fitting methodologies: if both fits are shown to have $\mathcal{R}_{bv}=1$,
the methodology with smaller magnitude of $\Delta_{\chi^{2}}$ could be
preferential. The same could be said for bias and variance however, bias in
particular is clearly closely related to $\Delta_{\chi^{2}}$.

## References

  * [1] R. D. Ball et al., The Path to Proton Structure at One-Percent Accuracy, arXiv:2109.02653.
  * [2] A. Stuart, Inverse problems: A bayesian perspective, Acta Numerica 19 451–559.
  * [3] NNPDF Collaboration, R. Abdul Khalek et al., A first determination of parton distributions with theoretical uncertainties, Eur. Phys. J. C (2019) 79:838, [arXiv:1905.04311].
  * [4] NNPDF Collaboration, R. Abdul Khalek et al., Parton Distributions with Theory Uncertainties: General Formalism and First Phenomenological Studies, Eur. Phys. J. C 79 (2019), no. 11 931, [arXiv:1906.10698].
  * [5] G. Backus and F. Gilbert, The resolving power of gross earth data, Geophysical Journal of the Royal Astronomical Society 16 169–205.
  * [6] M. S. Advani, A. M. Saxe, and H. Sompolinsky, High-dimensional dynamics of generalization error in neural networks, Neural Networks 132 (2020) 428–446, [arXiv:1710.03667].
  * [7] NNPDF Collaboration, R. D. Ball et al., Parton distributions for the LHC Run II, JHEP 04 (2015) 040, [arXiv:1410.8849].
  * [8] NNPDF Collaboration, R. D. Ball, E. R. Nocera, and R. L. Pearson, Nuclear Uncertainties in the Determination of Proton PDFs, Eur. Phys. J. C 79 (2019), no. 3 282, [arXiv:1812.09074].
  * [9] R. D. Ball, E. R. Nocera, and R. L. Pearson, Deuteron Uncertainties in the Determination of Proton PDFs, Eur. Phys. J. C 81 (2021), no. 1 37, [arXiv:2011.00009].
  * [10] R. Abdul Khalek, J. J. Ethier, J. Rojo, and G. van Weelden, nNNPDF2.0: quark flavor separation in nuclei from LHC data, JHEP 09 (2020) 183, [arXiv:2006.14629].
  * [11] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to machine learning for physicists, Physics Reports 810 (2019) 1–124, [arXiv:1803.08823]. A high-bias, low-variance introduction to Machine Learning for physicists.
  * [12] R. D. Ball, V. Bertone, S. Carrazza, L. D. Debbio, S. Forte, P. Groth-Merrild, A. Guffanti, N. P. Hartland, Z. Kassabov, J. I. Latorre, and et al., Parton distributions from high-precision collider data, The European Physical Journal C 77 (Oct, 2017) [arXiv:1706.00428].
  * [13] NNPDF Collaboration, R. D. Ball et al., An open-source machine learning framework for global analyses of parton distributions, arXiv:2109.02671.
  * [14] J. Cruz-Martinez, S. Forte, and E. Nocera, Future tests of parton distributions, Acta Physica Polonica B 52 (2021), no. 3 243, [arXiv:2103.08606].
  * [15] R. D. Ball, S. Carrazza, L. D. Debbio, S. Forte, Z. Kassabov, J. Rojo, E. Slade, and M. Ubiali, Precision determination of the strong coupling constant within a global pdf analysis, The European Physical Journal C 78 (May, 2018) [arXiv:1802.03398].
  * [16] F. Faura, S. Iranipour, E. R. Nocera, J. Rojo, and M. Ubiali, The strangest proton?, The European Physical Journal C 80 (Dec, 2020) [arXiv:2009.00014].
  * [17] S. Carrazza and J. Cruz-Martinez, Towards a new generation of parton densities with deep learning models, Eur. Phys. J. C 79 (2019), no. 8 676, [arXiv:1907.05075].
  * [18] B. Efron and R. J. Tibshirani, An introduction to the bootstrap.  CRC press, 1994.
  * [19] J. Luis Bernal and J. A. Peacock, Conservative cosmology: combining data with allowance for unknown systematics, Journal of Cosmology and Astroparticle Physics 2018 (Jul, 2018) 002–002, [arXiv:1803.04470].
  * [20] M. P. Hobson, S. L. Bridle, and O. Lahav, Combining cosmological data sets: hyperparameters and bayesian evidence, Monthly Notices of the Royal Astronomical Society 335 (Sep, 2002) 377–388, [astro-ph/0203259].
  * [21] J. Dove, B. Kerns, R. E. McClellan, S. Miyasaka, D. H. Morton, K. Nagai, S. Prasad, F. Sanftl, M. B. C. Scott, A. S. Tadepalli, and et al., The asymmetry of antimatter in the proton, Nature 590 (Feb, 2021) 561–565, [arXiv:2103.04024].
  * [22] 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. C77 (2017), no. 6 367, [arXiv:1612.03016].
  * [23] ATLAS Collaboration, M. Aaboud et al., Measurement of the Drell-Yan triple-differential cross section in $pp$ collisions at $\sqrt{s}=8$ TeV, JHEP 12 (2017) 059, [arXiv:1710.05167].
  * [24] ATLAS Collaboration, G. Aad et al., Measurement of the double-differential high-mass Drell-Yan cross section in pp collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 08 (2016) 009, [arXiv:1606.01736].
  * [25] ATLAS Collaboration, G. Aad et al., Measurement of $W^{\pm}$ and $Z$-boson production cross sections in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B759 (2016) 601–621, [arXiv:1603.09222].
  * [26] ATLAS Collaboration, M. Aaboud et al., Measurement of differential cross sections and $W^{+}/W^{-}$ cross-section ratios for $W$ boson production in association with jets at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 05 (2018) 077, [arXiv:1711.03296]. [Erratum: JHEP 10, 048 (2020)].
  * [27] ATLAS Collaboration, G. Aad et al., Measurement of the $t\bar{t}$ production cross-section in the lepton+jets channel at $\sqrt{s}=13$ TeV with the ATLAS experiment, Phys. Lett. B 810 (2020) 135797, [arXiv:2006.13076].
  * [28] ATLAS Collaboration, G. Aad et al., Measurements of top-quark pair differential cross-sections in the lepton+jets channel in $pp$ collisions at $\sqrt{s}=8$ TeV using the ATLAS detector, Eur. Phys. J. C76 (2016), no. 10 538, [arXiv:1511.04716].
  * [29] ATLAS Collaboration, M. Aaboud et al., Measurement of top quark pair differential cross-sections in the dilepton channel in $pp$ collisions at $\sqrt{s}$ = 7 and 8 TeV with ATLAS, Phys. Rev. D 94 (2016), no. 9 092003, [arXiv:1607.07281]. [Addendum: Phys.Rev.D 101, 119901 (2020)].
  * [30] ATLAS Collaboration, M. Aaboud et al., Measurement of the inclusive jet cross-sections in proton-proton collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 09 (2017) 020, [arXiv:1706.03192].
  * [31] ATLAS Collaboration Collaboration, G. Aad et al., Measurement of dijet cross sections in $pp$ collisions at 7 TeV centre-of-mass energy using the ATLAS detector, JHEP 1405 (2014) 059, [arXiv:1312.3524].
  * [32] M. Aaboud, G. Aad, B. Abbott, J. Abdallah, O. Abdinov, B. Abeloos, S. Abidi, O. AbouZeid, N. Abraham, H. Abramowicz, and et al., Measurement of the cross section for inclusive isolated-photon production in pp collisions at s=13 tev using the atlas detector, Physics Letters B 770 (Jul, 2017) 473–493, [arXiv:1701.06882].
  * [33] ATLAS Collaboration, G. Aad et al., Comprehensive measurements of $t$-channel single top-quark production cross sections at $\sqrt{s}=7$ TeV with the ATLAS detector, Phys. Rev. D 90 (2014), no. 11 112006, [arXiv:1406.7844].
  * [34] ATLAS Collaboration, M. Aaboud et al., Measurement of the inclusive cross-sections of single top-quark and top-antiquark $t$-channel production in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 04 (2017) 086, [arXiv:1609.03920].
  * [35] ATLAS Collaboration, M. Aaboud et al., Fiducial, total and differential cross-section measurements of $t$-channel single top-quark production in $pp$ collisions at 8 TeV using data collected by the ATLAS detector, Eur. Phys. J. C 77 (2017), no. 8 531, [arXiv:1702.02859].
  * [36] CMS Collaboration, S. Chatrchyan et al., Measurements of differential jet cross sections in proton-proton collisions at $\sqrt{s}=7$ TeV with the CMS detector, Phys.Rev. D87 (2013) 112002, [arXiv:1212.6660].
  * [37] CMS Collaboration, A. M. Sirunyan et al., Measurement of the triple-differential dijet cross section in proton-proton collisions at $\sqrt{s}=8\,\text{TeV}$ and constraints on parton distribution functions, Eur. Phys. J. C 77 (2017), no. 11 746, [arXiv:1705.02628].
  * [38] CMS Collaboration, A. M. Sirunyan et al., Measurement of the inclusive $\mathrm{t}\overline{\mathrm{t}}$ cross section in pp collisions at $\sqrt{s}=5.02$ TeV using final states with at least one charged lepton, JHEP 03 (2018) 115, [arXiv:1711.03143].
  * [39] CMS Collaboration, A. M. Sirunyan et al., Measurement of double-differential cross sections for top quark pair production in pp collisions at $\sqrt{s}=8$ TeV and impact on parton distribution functions, Eur. Phys. J. C77 (2017), no. 7 459, [arXiv:1703.01630].
  * [40] CMS Collaboration, A. M. Sirunyan et al., Measurement of differential cross sections for the production of top quark pairs and of additional jets in lepton+jets events from pp collisions at $\sqrt{s}=$ 13 TeV, Phys. Rev. D97 (2018), no. 11 112003, [arXiv:1803.08856].
  * [41] CMS Collaboration, A. M. Sirunyan et al., Measurements of $\mathrm{t\overline{t}}$ differential cross sections in proton-proton collisions at $\sqrt{s}=$ 13 TeV using events containing two leptons, JHEP 02 (2019) 149, [arXiv:1811.06625].
  * [42] CMS Collaboration, S. Chatrchyan et al., Measurement of the Single-Top-Quark $t$-Channel Cross Section in $pp$ Collisions at $\sqrt{s}=7$ TeV, JHEP 12 (2012) 035, [arXiv:1209.4533].
  * [43] CMS Collaboration, V. Khachatryan et al., Measurement of the t-channel single-top-quark production cross section and of the $\mid V_{tb}\mid$ CKM matrix element in pp collisions at $\sqrt{s}$= 8 TeV, JHEP 06 (2014) 090, [arXiv:1403.7366].
  * [44] CMS Collaboration, A. M. Sirunyan et al., Cross section measurement of $t$-channel single top quark production in pp collisions at $\sqrt{s}=$ 13 TeV, Phys. Lett. B 772 (2017) 752–776, [arXiv:1610.00678].
  * [45] LHCb Collaboration, R. Aaij et al., Measurement of the forward Z boson production cross-section in pp collisions at $\sqrt{s}=13$ TeV, JHEP 09 (2016) 136, [arXiv:1607.06495].