Nuclear data evaluation with Bayesian networks

Bayesian inference is a framework for inference under uncertainty. The central formula in Bayesian statistics is given by Bayes theorem, e.g., Sivia (1996),

$$P(H\,|\,E)=\frac{P(E\,|\,H)P(H)}{P(E)}\,,$$

(1)

where $H$ is an hypothesis and $E$ represents an observation. For instance, $H$ could be an hypothesis from the set $\{H_{1}:\textrm{`It rained'},H_{2}:\textrm{`It did not rain'}\}$ and the evidence $E$ an observation from a set of possible observations $\{\textrm{`Floor is wet'},\textrm{`Floor is dry'}\}$. The notation $P(.)$ denotes a probability, a number between zero and one, which is interpreted as degree of belief in Bayesian statistics. The probability $P(H)$ is referred to as prior probability and represents the degree of belief that hypothesis $H$ is true without taking into account the evidence $E$. The probability $P(E\,|\,H)$ denotes the probability for making the observation $E$ given that $H$ is true. For instance, $P(\textrm{Floor is wet}\,|\,\textrm{It rained})=1$ if we know that rain always causes the floor to be wet. The probability $P(E)$ is referred to as marginal likelihood or model evidence and is the probability of making a specific observation $P(E)$ based on the modeling assumptions. The probability $P(H\,|\,E)$ is referred to as posterior probability and represents the degree of belief that $H$ is true after the observation $E$ has been made.

The sets we have introduced above are discrete and we can exhaustively enumerate their elements. However, for continuous quantities, such as temperatures or cross sections in nuclear data, the notion of probability needs to be replaced by the concept of probability density function. Given a probability density function $\rho(x)$, the probability that a value is within the range enclosed by $a$ and $b$ is given by

$$P(a\leq x\leq b)=\int_{a}^{b}\rho(x)dx\,.$$

(2)

Due to this relationship, probability density functions must be non-negative everywhere and their integral over the full range one. We do not emphasize the distinction between probabilities and probability density functions in the following because this aspect is not crucial for the concept of Bayesian networks.

Even though the Bayesian theorem in its basic form is at the core of Bayesian inference as it formalizes the mechanics of learning from experience, it only deals with a very simple case. In the following we generalize the perspective on the inference problem to prepare the discussion of Bayesian networks. The distinction between hypotheses and evidence is to some extent arbitrary. For instance, the possible observation ‘Floor is wet’ we mentioned as an example is also a hypothesis, but we named it evidence because it is an hypothesis that we observed to be true. In contrast, which hypothesis $H$ among the possible ones is true was not directly observed and Bayes theorem allows us to indirectly improve our knowledge about their likelihood by using the evidence.

Therefore, at a higher level of abstraction, we can equally say that we have a system of hypotheses and the Bayesian inference framework provides a mechanism to update our knowledge about their likelihood by observing some of the hypotheses to be true or false. Assume that we have a number of variables $H_{1},H_{2},H_{3},\dots$ with each one being an hypothesis from a certain class of hypotheses. For instance, $H_{1}$ could be the value of the elastic cross section, $H_{2}$ the value of the non-elastic cross section and $H_{3}$ the value of the total cross section. Their joint prior probability density function is denoted by $P(H_{1},H_{2},H_{3})$. To give an example of the form of the Bayesian update formula in the case of three variables, assume that we know the elastic cross section $H_{1}$. The Bayesian update formula takes then the form

$$P(H_{2},H_{3}\,|\,H_{1})=\frac{P(H_{1}\,|\,H_{2},H_{3})P(H_{2},H_{3})}{P(H_{1})}\,.$$

(3)

Even though the Bayesian update formula can be written down for an arbitrary number of variables in this way, typical inference tasks with many variables pose computational challenges. A relevant inference query is to find an assignment of values to $H_{2}$ and $H_{3}$ that maximizes the value of the posterior probability density function $P(H_{2},H_{3}\,|\,H_{1})$. This assignment is called a maximum a posteriori probability (MAP) estimate.

The idea of Bayesian networks is to exploit the structure of the relationships between variables. Every joint probability distribution can be decomposed into a product of conditional probability distributions. For instance, for three variables, we can write the joint probability density function as

$$P(H_{1},H_{2},H_{3})=P(H_{1})P(H_{2}\,|\,H_{1})P(H_{3}\,|\,H_{1},H_{2})\,.$$

(4)

In this product, going from left to right, each conditional probability distribution is conditioned on all variables that appeared before. These conditional dependencies can be visualized as a Bayesian network:

Variables (or a collection of variables of the same type) are associated with nodes and directed arrows go from the variables being conditioned on to the dependent variables. Bayesian networks are not allowed to have cycles, hence it is impossible to arrive at the same node twice by following the arrows along their pointing direction.

In the product in eq. 4 conditional dependencies can be dropped if variables are conditionally independent. Two random variables $X$ and $Y$ are conditional independent given a third random variable $Z$ if

$$P(X\,|\,Y,Z)=P(X\,|\,Z)\;\;\textrm{and}\;\;P(Y\,|\,X,Z)=P(Y\,|\,Z)\,.$$

(5)

In words, knowledge about the value of $Y$ does not improve the knowledge about the value of $X$ if the value of $Z$ is already known.

In the example of the three hypotheses: If $H_{2}$ and $H_{3}$ are conditionally independent given $H_{1}$, eq. 4 simplifies to

$$P(H_{1},H_{2},H_{3})=P(H_{1})P(H_{2}\,|\,H_{1})P(H_{3}\,|\,H_{1})\,.$$

(6)

and the corresponding Bayesian network loses one arrow:

Please note that there is not a unique way to factorize a joint probability distribution. For instance, for two variables we can choose $P(H_{1},H_{2})=P(H_{1})P(H_{2}\,|\,H_{1})$ or $P(H_{1},H_{2})=P(H_{2})P(H_{1}\,|\,H_{2})$. Consequently, there are several Bayesian network topologies to describe the same joint probability distribution which usually differ in both the number of arrows and their orientations. An objective in Bayesian network modeling is therefore to find a topology which is simple, e.g., possesses a smaller number of connections or a topology with arrows reflecting the causal relationships between variables. The specific topology of Bayesian network can be exploited in Bayesian inference.

Bayesian networks may be classified according to whether variables are continuous or discrete in the network or the type of prior distribution imposed on variables, e.g., Gaussian distributions. Three general inference tasks can be identified, e.g., Koller and Friedman (2009):

1. 1.

   Determinig the links between variables and their orientation,

2. 2.

   Estimating the parameters of the conditional prior distributions, and,

3. 3.

   Inferring the values of unobserved variables.

In this paper, we restrict ourselves to Bayesian networks with all prior conditional distributions being multivariate normal and possibly non-linear relationships between variables. We deal only with the inference of unobserved variables (3), which is the essential inference task in nuclear data evaluation. To develop the concepts and notation, we first review the Generalized Least Squares method, already linking it to Bayesian networks, and afterwards extend the discussion to nested and non-linear relationships.

The Bayesian version of the Generalized Least Squares (GLS) method, e.g., Muir (1989), is usually regarded as a method to update the parameters of a model based on data from experiments affected by statistical and systematic errors. In the nuclear data context, parameters may be the cross sections associated with the points of an energy mesh or the parameters of a nuclear physics model. This view suggests that model parameters and experimental errors are quantities of different nature both from the point of physics and Bayesian statistics. However, from the point of Bayesian statistics, all these quantities are variables with uncertain values and can be treated in the same way. In order to introduce the Bayesian network interpretation, we present the GLS formulas without regard to the particular meaning of the variables to emphasize this symmetry.

The GLS method is based on the multivariate normal distribution, which is of the form

$$\mathcal{N}(\vec{x}\,|\,\vec{\mu},\mathbf{\Sigma})=\frac{1}{\sqrt{(2\pi)^{N}\det\mathbf{\Sigma}}}\\ \exp\left(-\frac{1}{2}(\vec{x}-\vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}-\vec{\mu})\right)\,,$$

(7)

and characterized by a mean vector $\vec{\mu}$ and a covariance matrix $\mathbf{\Sigma}$. We use the notation $\vec{x}\sim\mathcal{N}(\vec{\mu},\mathbf{\Sigma})$ to indicate that the random vector $\vec{x}$ is distributed according a multivariate normal distribution with mean vector $\vec{\mu}$ and covariance matrix $\mathbf{\Sigma}$, hence its probability density function given by eq. 7.

Assume we have a set of variables $\{y_{i}\}_{i=1..N}$ whose values are uncertain and we assemble these variables to a vector $\vec{y}$ . Some variables are either completely or to a certain degree determined by other variables. For instance, the values of the parameters of a deterministic nuclear model determine completely the resulting predictions. In contrast, for a stochastic nuclear model, the model parameters certainly influence the predictions but do not completely determine them due to the stochastic nature of the simulation. In a nuclear experiment, the measurement is influenced by the underlying true values of nuclear properties but not completely determined due to various measurement errors.

To account for these dependencies, we partition the variables in $\vec{y}$ into two subvectors, $\vec{y}_{I}$ and $\vec{y}_{J}$. We refer to variables whose indices are in the set $I$ as independent variables and assume them to be governed by a multivariate normal distribution, i.e., $\vec{y}_{I}\sim\mathcal{N}(\vec{u}_{I},\mathbf{U}_{I,I})$. The variables with their indices being in $J$ are referred to as dependent variables and they are assumed to be functions of the independent variables,

$$\vec{y}_{J}=f(\vec{y}_{I},\vec{\tau})=g(\vec{y}_{I})+\vec{\tau}\,.$$

(8)

We introduced the random variable $\vec{\tau}$ to allow for stochastic perturbation of the otherwise deterministic link. The vectors $\vec{\tau}$ and $\vec{y}_{J}$ are of the same size. Please note that the variables in $\vec{\tau}$ have also to be regarded as independent variables, even though their indices are not included in $I$ for notational reasons. In nuclear data evaluation, $\vec{\tau}$ typically contains the statistical errors. We assume $\vec{\tau}$ to be governed by a multivariate normal distribution, i.e., $\vec{\tau}\sim\mathcal{N}(\vec{u}_{J},\mathbf{U}_{J,J})$, and the random variable $\vec{y}_{I}$ to be independent of $\vec{\tau}$. With these specifications, the joint probability distribution $P(\vec{y}_{J},\vec{y}_{I},\vec{\tau})$ factorizes therefore into

$$P(\vec{y}_{J},\vec{y}_{I},\vec{\tau})=P(\vec{y}_{I})P(\vec{\tau})P(\vec{y}_{J}\,|\,\vec{y}_{I},\vec{\tau})$$

(9)

with the following corresponding Bayesian network:

A non-linear relationships $g(\vec{y}_{I},\vec{\tau})$ in eq. 8 can be cast into the GLS framework by constructing a linear Taylor approximation. A linear Taylor approximation of eq. 8 is of the form

$$\vec{y}_{J}=\vec{y}_{\textrm{ref},J}+\mathbf{T}\,\left(\vec{y}_{I}-\vec{y}_{\textrm{ref},{I}}\right)+(\vec{\tau}-\vec{\tau}_{\textrm{ref}})$$

(10)

where $\vec{y}_{\textrm{ref},J}=f(\vec{y}_{\textrm{ref},I},\vec{\tau}_{\textrm{ref}})$ and $\mathbf{T}$ is the Jacobian matrix of $f$ evaluated at $\vec{y}_{\textrm{ref},I}$. The Jacobian matrix is somestimes also referred to as sensitivity matrix. At this point it may be regarded as unnecessary to introduce a potentially non-zero vector $\vec{\tau}_{\textrm{ref}}$ because $\vec{\tau}$ is an additive (i.e., also linear) contribution to the result. However, this term will become important later in the treatment of nested relationships.

Depending on the specific construction of $\mathbf{T}$, this equation may implement interpolation using splines, a Legendre polynomial, a Fourier polynomial or any other type of polynomial. For instance, for spline interpolation $\mathbf{T}$ may map the values at the knot points given in $\vec{y}_{I}$ to the locations of the data points $\vec{y}_{J}$. In the case of a Legendre polynomial, the Jacobian $\mathbf{T}$ maps the values of the Legendre coefficients given in $\vec{y}_{I}$ to the function values at the positions of the data points whose values are given in $\vec{y}_{J}$. Equation 10 is therefore general and can accomodate many different types of interpolation schemes. Importantly, this equation as a building block does not only cover interpolation but, for instance, is also used to distribute normalization errors given in $\vec{y}_{I}$ to the data points in $\vec{y}_{J}$.

It will be helpful to have a compact formula to express $\vec{y}$, i.e., both $\vec{y}_{I}$ and $\vec{y}_{J}$, as a function of the independent variables $\vec{y}_{I}$ and $\vec{\tau}$. In other words, we want to extend eq. 10 so that its result is the full vector $\vec{y}$. To this end, we introduce a vector $\vec{z}$ of the same size as $\vec{y}$ whose components are given by

$$\vec{z}_{I}=\vec{y}_{I}\;\;\;\;\textrm{and}\;\;\;\;\vec{z}_{J}=\vec{\tau}\,.$$

(11)

This vector $\vec{z}$ contains all independent variables, which completely determine the values in $\vec{y}_{J}$.

Analogously, we also introduce a vector $\vec{z}_{\textrm{ref}}$ with elements given by

$$\vec{z}_{\textrm{ref},I}=\vec{y}_{\textrm{ref},I}\;\;\;\;\textrm{and}\;\;\;\;\vec{z}_{\textrm{ref},J}=\vec{\tau}_{\textrm{ref}}\,.$$

(12)

Finally, we regard $\mathbf{T}$ as the respective block in a larger matrix $\mathbf{S}$ whose blocks are defined as

$$\mathbf{S}_{I,I}=\mathbb{1},\;\;\;\;\mathbf{S}_{I,J}=\mathbf{0}\;\;\;\;\mathbf{S}_{J,I}=\mathbf{T},\;\;\;\;\mathbf{S}_{J,J}=\mathbb{1}\,.$$

(13)

Please note that even though we use the term block, the ordering of the rows and columns does not matter, as long as the index sets $I$ and $J$ do not contain common indices and all indices from $1$ to $N$ are present in $A=I\cup J$. The same is true for the vectors $\vec{z}$ and $\vec{z}_{\textrm{ref}}$.

The introduction of these quantities allows us to rewrite eq. 10 to cover both dependent and independent variables,

$$\vec{y}=\vec{y}_{\textrm{ref}}+\mathbf{S}\,\left(\vec{z}-\vec{z}_{\textrm{ref}}\right)\,.$$

(14)

Regarding the Bayesian network interpretation, this equation propagates the values of the variables without parent nodes, i.e., $\vec{y}_{I}$ and $\vec{\tau}$ to all the other variables. Now with all variables and their mutual relationships defined, we can deal with the question of inference in the Bayesian framework.

Bayesian statistics can be seen as a framework to update the knowledge about the possible values of variables based on observation of some variables. Using our notation, the Bayesian inference formula to achieve this update is given by

$$P(\vec{y}_{I}\,|\,\vec{y}_{J})=\frac{1}{P(\vec{y}_{I})}P(\vec{y}_{J}\,|\,\vec{y}_{I})P(\vec{y}_{I})$$

(15)

where we marginalized over $\vec{\tau}$ to treat it as a nuisance variable, i.e., $P(\vec{y}_{J}\,|\,\vec{y}_{I})=\int P(\vec{\tau})P(\vec{y}_{J}\,|\,\vec{y}_{I},\vec{\tau})d\vec{\tau}$.

The GLS method is a special case of the general Bayesian update formula to update the knowledge about $\vec{y}_{I}$ based on the observed values of $\vec{y}_{J}$ in eq. 8 if the following assumptions are met:

• •

  The function $f(.)$ in eq. 8 is linear, i.e., it can be exactly represented by eq. 10

• •

  The prior distribution of $\vec{y}_{I}$ and $\vec{\tau}$ are multivariate normal and no prior correlations exist between elements in $\vec{y}_{I}$ and elements in $\vec{\tau}$

The resulting posterior distribution $P(\vec{y}_{I}\,|\,\vec{y}_{J})$ is then also multivariate normal and the posterior mean vector and covariance matrix can be obtained by analytic matrix formulas.

To state the GLS formulas, we introduce the prior mean vector $\vec{u}$ and prior covariance matrix $\mathbf{U}$ which characterize the prior on $\vec{y}_{I}$ and $\vec{\tau}$ simultaneously. The priors on these vectors have been already introduced along with eq. 8, which we repeat here for the convenience of the reader: $\vec{y}_{I}\sim\mathcal{N}(\vec{u}_{I},\mathbf{U}_{I,I})$ and $\vec{\tau}\sim\mathcal{N}(\vec{u}_{J},\mathbf{U}_{J,J})$. The matrices $\mathbf{U}_{I,I}$ and $\mathbf{U}_{J,J}$ are now considered to be the blocks of the compound covariance matrix $\mathbf{U}$ at the positions defined by index sets $I$ and $J$, respectively. Due to the independence assumption between $\vec{y}_{I}$ and $\vec{\tau}$, we have $\mathbf{U}_{I,J}=\mathbf{U}_{J,I}^{T}=\mathbf{0}$. Let us assume that we have observed the values of the vector $\vec{y}_{J}$, i.e., $\vec{y}_{J}=\vec{r}$. Observed nodes in the Bayesian network are colored grey:

The update formula to obtain the posterior covariance matrix $\mathbf{U}_{I,I}^{\prime}$ is given by

$$\mathbf{U}_{I,I}^{\prime}=\left(\mathbf{S}_{J,I}^{T}\mathbf{U}_{J,J}^{-1}\mathbf{S}_{J,I}+\mathbf{U}_{I,I}^{-1}\right)^{-1}\,.$$

(16)

and the update formula for the posterior mean vector by

$$\vec{u}_{I}^{\prime}=\vec{y}_{\textrm{ref},I}+\mathbf{U}_{I,I}^{\prime}\Big{(}\mathbf{S}_{J,I}^{T}\mathbf{U}_{J,J}^{-1}(\vec{r}-\vec{y}_{\textrm{ref},J})\\ +\mathbf{U}_{I,I}^{-1}(\vec{u}_{I}-\vec{y}_{\textrm{ref},I})\Big{)}\,.$$

(17)

A derivation of eqs. 16 and 17 can be found in appendix A.

These equations can be written in more compact form. The assumption of zero cross-covariance elements between $\vec{y}_{I}$ and $\vec{\tau}$, i.e., $\mathbf{U}_{I,J}=\mathbf{U}_{J,I}^{T}=\mathbf{0}$, allows us to write

$$\mathbf{U}_{I,I}^{\prime}=\left(\mathbf{S}_{A,I}^{T}\mathbf{U}_{A,A}^{-1}\mathbf{S}_{A,I}\right)^{-1}\,.$$

(18)

If in addition taking into account the structure of $\mathbf{S}$ in eq. 13 and defining the vector $\vec{v}$ with elements

$$\vec{v}_{I}=\vec{u}_{I},\,\vec{v}_{J}=\vec{r}\,,$$

(19)

we can rewrite eq. 17 as

$$\vec{u}_{I}^{\prime}=\vec{y}_{\textrm{ref},I}+\mathbf{U}_{I,I}^{\prime}\mathbf{S}_{A,I}^{T}\mathbf{U}_{A,A}^{-1}(\vec{v}-\vec{y}_{\textrm{ref}})\,.$$

(20)

The vector $\vec{u}_{I}^{\prime}$ and the covariance matrix $\mathbf{U}_{I,I}^{\prime}$ are the updated distribution parameters of the multivariate normal distribution associated with $\vec{y}_{I}$ due to observing the values in $\vec{y}_{J}$, i.e., $\vec{y}_{I}\sim\mathcal{N}(\vec{u}_{I}^{\prime},\mathbf{U}_{I,I}^{\prime})$

To discuss the application of the GLS method for a Bayesian network with more nodes, we consider the following scenario for the purpose of illustration. So far we did not impose further constraints on the prior of $\vec{\tau}$ except being a multivariate normal distribution. If, however, $\vec{\tau}$ can be decomposed into two a priori mutually independent blocks $\vec{\tau}_{K},\vec{\tau}_{L}$, i.e., $\mathbf{U}_{K,L}=(\mathbf{U}_{L,K})^{T}=\mathbf{0}$, we may split up eq. 10 into two equations:

	$\displaystyle\vec{y}_{K}$	$\displaystyle=\vec{y}_{\textrm{ref},K}+\mathbf{T}_{K,I}\,\left(\vec{y}_{I}-\vec{y}_{\textrm{ref},{I}}\right)+(\vec{\tau}_{K}-\vec{\tau}_{\textrm{ref},K})\,,$		(21)
	$\displaystyle\vec{y}_{L}$	$\displaystyle=\vec{y}_{\textrm{ref},L}+\mathbf{T}_{L,I}\,\left(\vec{y}_{I}-\vec{y}_{\textrm{ref},{I}}\right)+(\vec{\tau}_{L}-\vec{\tau}_{\textrm{ref},L})\,.$		(22)

If we regard the values in $\vec{y}_{I}$ as fixed, the remaining contribution to $\vec{y}_{K}$ and $\vec{y}_{L}$ given by $\vec{\tau}_{K}$ and $\vec{\tau}_{L}$, respectively, are independent of each other due to our independence assumption for the blocks in $\vec{\tau}$ above. Due to this conditional independence of $\vec{y}_{K}$ and $\vec{y}_{L}$ there is no arrow between $\vec{y}_{K}$ and $\vec{y}_{L}$. Considering this conditional independence and the values in $\vec{y}_{K}$ and $\vec{y}_{L}$ being given as a function of $\vec{y}_{I}$, we have the following Bayesian network:

In this scenario, one may assume that only the variables in $\vec{y}_{K}$ have been observed, i.e., the variables in the vector $\vec{y}$ referred to by index set K.

Keeping this example in mind, we want to discuss the general form of the GLS equations for a Bayesian network without intermediate nodes in which some of the terminal nodes have been observed, which amounts to knowing the values of some variables in the vector $\vec{y}_{J}$ in eq. 10. The treatment of nested relationship will be discussed later in section II.5.

We denote by $D$ the index set referring to the observed values in $\vec{y}$, i.e., $\vec{y}_{D}$ is the subvector with observed values. For instance, in the example above we can identify $D=K$. Let us denote by $\neg{D}$ the index set with all indices that are not in $D$, i.e., $\neg{D}=A\setminus D$ with $A=\{1,\dots,N\}$. In this case, we can consider eq. 14 to realize that the variables in the vector $\vec{z}_{\neg D}$, which includes all variables in $\vec{y}_{I}$ and a subset of the variables in $\vec{\tau}$, can be regarded as the independent variables. The dependent and measured variables are now those in $\vec{y}_{D}$. For this more general case, the structure of the update formulas remains the same.

For the posterior covariance matrix we have

$$\mathbf{U}_{\neg D,\neg D}^{\prime}=\left(\mathbf{S}_{D,\neg D}^{T}\mathbf{U}_{D,D}^{-1}\mathbf{S}_{\neg D,D}+\mathbf{U}_{\neg D,\neg D}^{-1}\right)^{-1}\,,$$

(23)

and for the posterior mean vector

$$\vec{u}_{\neg D}^{\prime}=\vec{z}_{\textrm{ref},\neg D}+\mathbf{U}_{\neg D,\neg D}^{\prime}\Big{(}\mathbf{S}_{D,\neg D}^{T}\mathbf{U}_{D,D}^{-1}(\vec{r}-\vec{y}_{\textrm{ref},D})\\ +\mathbf{U}_{\neg D,\neg D}^{-1}(\vec{u}_{\neg D}-\vec{z}_{\textrm{ref},\neg D})\Big{)}\,.$$

(24)

If we assume that the random variables in $\vec{z}_{D}$ are independent of those in $\vec{z}_{\neg D}$, we can again simplify the equations. Later when dealing with Bayesian networks with more nodes, this independence assumption will be fulfilled because variables of the form $\vec{\tau}$ associated with different nodes are a priori independent and all variables associated with each single node are either observed or unobserved at the same time. We do not consider situations in which only a part of the variables associated with a node is observed.

If the observed values are $\vec{y}_{D}=\vec{r}$, the vector $\vec{v}$ is now defined by

$$\vec{v}_{\neg D}=\vec{u}_{\neg D},\;\vec{v}_{D}=\vec{r}\,.$$

(25)

In addition we define a vector $\vec{v}_{\textrm{ref}}$ as

$$\vec{v}_{\textrm{ref},\neg D}=\vec{z}_{\textrm{ref},\neg D},\;\vec{v}_{\textrm{ref},D}=\vec{y}_{\textrm{ref},D}\,.$$

(26)

With these specifications, the posterior covariance matrix of the independent variables $\vec{z}_{\neg D}$ can be written as

$$\mathbf{U}_{\neg D,\neg D}^{\prime}=\left(\mathbf{S}_{A,\neg D}^{T}\mathbf{U}_{A,A}^{-1}\mathbf{S}_{A,\neg D}\right)^{-1}\,,$$

(27)

and the posterior mean vector as

$$\vec{u}_{\neg D}^{\prime}=\vec{z}_{\textrm{ref},\neg D}+\mathbf{U}_{\neg D,\neg D}^{\prime}\mathbf{S}_{D,A}^{T}\mathbf{U}_{A,A}^{-1}(\vec{v}-\vec{v}_{\textrm{ref}})\,.$$

(28)

The outer inversion in eq. 27 has never to be evaluated explicitly. In section II.3, we discuss how blocks of the posterior covariance matrix of interest can be computed and how to draw samples from the posterior distribution. The multiplication with $\mathbf{U}_{\neg D,\neg D}^{\prime}$ in eq. 28 can be reformulated as the task to solve a system of linear equations,

	$\displaystyle\mathbf{A}\vec{x}=\vec{b}\;\;\;\textrm{with}\;\;\;\mathbf{A}$	$\displaystyle=\mathbf{S}_{A,\neg D}^{T}\mathbf{U}_{A,A}^{-1}\mathbf{S}_{A,\neg D}$
	$\displaystyle\vec{b}$	$\displaystyle=\mathbf{S}_{D,A}^{T}\mathbf{U}_{A,A}^{-1}(\vec{v}-\vec{v}_{\textrm{ref}})\,,$

with the solution vector $\vec{x}$ representing the result of the matrix product in eq. 28. As all matrices $\mathbf{U}^{-1}$, $\mathbf{S}$, and $\mathbf{S}_{A,\neg D}^{T}\mathbf{U}_{A,A}^{-1}\mathbf{S}_{A,\neg D}$ are typically sparse, we employ sparse matrix algorithms to greatly speed up the computation of the posterior expectation vector $\vec{u}_{\neg D}^{\prime}$.

Please note that so far we have only recapitulated one version of the formulas of the Generalized Least Squares method and introduced specific notation. These formulas will be the basis for inference in Bayesian networks with linear relationships and with a slight modification explained in section II.6 also in the case of non-linear relationships. Using these formulas implies that observed nodes are always attached to a vector $\vec{\tau}$ with non-zero prior uncertainties. We feel that this assumption is not a serious limitation in the context of nuclear data evaluation because measurements are always associated with a statistical uncertainty component. If this is not the case, the prior uncertainty associated with the elements of $\vec{\tau}$ can be made very small which effectively eliminates the impact of these noise node on the results although it may be detrimental to numerical stability.

Deterministic relationships, i.e., those associated with a zero diagonal element in the covariance matrix $\mathbf{U}$, must be treated as a special case because $\mathbf{U}^{-1}$ is not defined anymore. Assume that $Z$ is the subset of $\neg D$ associated with zero prior uncertainties and $\neg Z$ is the subset of $\neg D$ with non-zero prior uncertainties. The approach to deal with deterministic relationships is to use $\neg Z$ instead of $\neg D$ in eq. 28 and eq. 27. Afterward the full posterior mean vector $\vec{u}_{\neg D}$ and covariance matrix $\mathbf{U}_{\neg D,\neg D}$ can be obtained by imputing the prior mean vector and covariance matrix to the missing blocks,

$$\vec{u}_{Z}^{\prime}=\vec{u}_{Z}\;\;\;\textrm{and}\;\;\;\mathbf{U}_{Z,Z}^{\prime}=\mathbf{0},\mathbf{U}_{Z,\neg{Z}}=(\mathbf{U}_{\neg Z,Z})^{T}=\mathbf{0}\,.$$

(29)

In order to obtain the posterior expectation $\vec{u}_{D}$ of $\vec{z}_{D}$, we note that the known vector $\vec{r}$ of the observed node is given by the sum of the posterior expectation $\vec{u}_{D}$ and the posterior expectations in $\vec{u}_{\neg D}$ propagated to the observed node by means of eq. 14. With the result of the propagation given by

$$g_{D}(\vec{u}_{\neg D}):=\vec{y}_{\textrm{ref},\neg D}+\mathbf{S}_{\neg D,D}(\vec{u}_{\neg D}^{\prime}-\vec{z}_{\textrm{ref},\neg D})$$

(30)

we can therefore compute the posterior expectation $\vec{u}_{D}^{\prime}$ by

$$\vec{u}_{D}^{\prime}=\vec{r}-g_{D}(\vec{u}_{\neg D})\,.$$

(31)

To compute the posterior covariance matrix blocks $\mathbf{U}_{D,D}$ and $\mathbf{U}_{\neg D,D}$, we remark that the following relationships for the covariance matrix of random vectors $\vec{x},\vec{y},\vec{z}$ hold:

	$\displaystyle\operatorname{Cov}\left[\vec{x},\vec{y}\right]$	$\displaystyle=\left(\operatorname{Cov}\left[\vec{y},\vec{x}\right]\right)^{T},$
	$\displaystyle\operatorname{Cov}\left[\vec{x}+\vec{y},\vec{z}\right]$	$\displaystyle=\operatorname{Cov}\left[\vec{x},\vec{z}\right]+\operatorname{Cov}\left[\vec{y},\vec{z}\right],$
	$\displaystyle\operatorname{Cov}\left[\mathbf{S}\vec{x},\vec{y}\right]$	$\displaystyle=\mathbf{S}\operatorname{Cov}\left[\vec{x},\vec{y}\right]\,,$

where the notation $\operatorname{Cov}\left[\vec{x},\vec{y}\right]$ denotes the covariance matrix between random vectors $\vec{x}$ and $\vec{y}$.

Introducing the abbreviation $\operatorname{Var}\left[\vec{x}\right]:=\operatorname{Cov}\left[\vec{x},\vec{x}\right]$ and calculating the variance of both sides of the equation $\vec{z}_{D}=\vec{r}-g_{D}(\vec{z}_{\neg D})$, we find that $\operatorname{Var}\left[\vec{z}_{D}\right]=\operatorname{Var}\left[g_{D}(\vec{z}_{\neg D})\right]$ and in terms of the corresponding posterior covariance matrix blocks:

$$\mathbf{U}_{D,D}^{\prime}=\mathbf{S}_{D,\neg D}\mathbf{U}_{\neg D,\neg D}^{\prime}\left(\mathbf{S}_{D,\neg D}\right)^{T}$$

(32)

The application of the variance operator on both sides of the rearranged equation $\vec{z}_{D}+g_{D}(\vec{z}_{\neg D})=\vec{r}$ yields

$$\mathbf{U}_{D,\neg D}^{\prime}=\left(\mathbf{U}_{\neg D,D}^{\prime}\right)^{T}=\mathbf{S}_{D,\neg D}\mathbf{U}_{\neg D,\neg D}^{\prime}\,.$$

(33)

In summary, with eqs. 27, 28, 29, 30, 31, 32 and 33 we obtain the posterior expectation $\vec{u}^{\prime}$ and covariance matrix $\mathbf{U}^{\prime}$ for the full vector of independent variables $\vec{z}$ by conditioning on the known vector $\vec{r}$ of one (or several) observed nodes for both stochastic and deterministic linear functional relationship between variables. The posterior expectation $\vec{y}_{\textrm{post}}$ associated with the dependent nodes in $\vec{y}$ can be obtained by using eq. 14 with the posterior expectation of independent variables:

$$\vec{y}_{\textrm{post}}=\vec{y}_{\textrm{ref}}+\mathbf{S}\left(\vec{u}^{\prime}-\vec{z}_{\textrm{ref}}\right)$$

(34)

The posterior covariance matrix $\mathbf{Y}_{\textrm{post}}$ of dependent nodes is given by

$$\mathbf{Y}_{\textrm{post}}=\mathbf{S}\mathbf{U}^{\prime}\mathbf{S}^{T}\,.$$

(35)

The important feature of these formulas is that all matrices involved are typically sparse in nuclear data evaluation if we do not absorb functional relationships, such as those for normalization errors into the covariance matrix but explicitly keep them in $\mathbf{S}$. The matrix $\mathbf{S}$ contains the information about the connection structure of the Bayesian network. In this paper, we do not need specialized algorithms for Bayesian networks to exploit the connection structure but can rely on general libraries for sparse matrix operations instead. In the next section, we elaborate on ways to deal with the large posterior covariance matrix by exploiting sparsity. Afterwards, we focus on the Bayesian network interpretation and discuss the proper handling of nested and non-linear relationships.

When working with ten or even hundred thousands of variables, it is often not practical to compute or store the full posterior covariance matrix explicitly. Even if we just need to evaluate the posterior covariance block $\mathbf{U}_{I,I}^{\prime}$ of independent variables, it may be still too large in practice.

In nuclear data evaluation, the ability to compute at least the posterior covariance matrix block associated with the evaluated cross sections is very important because this information is needed for linear error propagation through application codes. However, sometimes linear error propagation based on the covariance matirx is not feasible with a reasonable amount of computational resources. For instance, if we have to use a very fine energy mesh, e.g., several tens of thousands of mesh points, in energy regions with quickly oscillating cross sections, even the handling of the posterior covariance matrix block associated with cross sections becomes impractical. In such cases, we may prefer to draw samples from the posterior distribution and to subsequently propagate these samples through application codes according to the Total Monte Carlo (TMC) method Koning and Rochman (2008); Rochman et al. (2009, 2010).

Therefore, feasible approaches to evaluate blocks of the posterior covariance matrix and to sample from the posterior distribution using the full posterior covariance matrix are both important in the field of nuclear data.

Following we explain the computation of blocks from the posterior covariance matrix and the sampling using the full posterior covariance matrix. For both requirements, we exploit the key insight that the inverse of the posterior covariance matrix in eq. 27,

$$\left(\mathbf{U}_{\neg D,\neg D}^{\prime}\right)^{-1}=\mathbf{S}_{A,\neg D}^{T}\mathbf{U}_{A,A}^{-1}\mathbf{S}_{A,\neg D}\,,$$

(36)

can be very sparse even though the posterior covariance matrix is not. The reason is that off-diagonal elements in the posterior covariance matrix are non-zero if the associated variables are correlated, irrespective of whether this correlation is mediated by a third variable. In contrast to that, elements in the inverse posterior covariance matrix are zero if the associated variables are conditionally independent given all other variables, e.g., (Koller and Friedman, 2009, Chap. 7.1.3). In other words, if two variables are only correlated because they are both influenced by one or a set of other variables, they are conditionally independent. The approaches described in the following to obtain random samples from the posterior distribution and to calculate blocks of the full posterior covariance matrix rely on the property that the inverse posterior covariance matrix is sparse. This sparseness property is present in all the Bayesian network examples we present in sections III.1, III.2 and III.3 thanks to the sparse Gaussian process construction explained in section II.7 and the specific handling of normalization errors.

First we discuss one way to obtain samples from the posterior distribution that is amenable to optimization by exploiting sparsity. To sample from a multivariate normal distribution $\vec{x}\sim\mathcal{N}(\vec{\mu},\mathbf{\Sigma})$, we can first compute a Cholesky decomposition of the inverse covariance matrix $\mathbf{\Sigma}^{-1}=\mathbf{L}^{T}\mathbf{L}$ with $\mathbf{L}$ being an upper triagonal matrix. Making use of the fact that the elements in $\vec{z}=\mathbf{L}\vec{x}$ are distributed according to a standard normal distribution, i.e., $z_{i}\sim\mathcal{N}(0,1)$, we can generate a vector $\vec{z}$ with elements from a standard normal distribution and solve the set of linear equations $\mathbf{L}\vec{y}=\vec{z}$. Because $\mathbf{L}$ is an upper triagonal matrix, this system of linear equations can be efficiently solved using backward substitution, e.g., (Press, 2007, Chap. 2.2). A sample from the multivariate normal distribution is then given by $\vec{s}=\vec{\mu}+\vec{y}$.

Therefore, we can make a sparse Cholesky decomposition of $\left(\mathbf{U}_{\neg D,\neg D}^{\prime}\right)^{-1}$ in order to apply the approach described in the previous paragraph to obtain samples from the posterior distribution. The essential idea of a sparse Cholesky decomposition of a sparse matrix $\mathbf{A}$ is to find a permutation matrix $\mathbf{P}$ so that the Cholesky decomposition of the permuted matrix $\mathbf{P}\mathbf{A}\mathbf{P}^{T}$ is sparse. Several algorithms exist to find such a so-called fill-reducing permutation, e.g., Davis (2006), and a popular one is the minimum-degree algorithm, e.g., George and Liu (1989). Without a permutation applied before the Cholesky decomposition, the resulting Cholesky factor $\mathbf{L}$ is in general not sparse. We employed the Matrix package Bates and Maechler (2021) available in the programming language R, which itself relies on the CHOLMOD routines Chen et al. (2008) written in C to perform the sparse Cholesky decomposition. The sparsity of the Cholesky factor $\mathbf{L}$ does not only help to keep the storage requirement low but also speeds up the solution of the system of linear equations $\mathbf{L}\vec{x}=\vec{z}$. Finally, because we only have obtained a sample of the variables associated with indices $\neg D$, we need to compute the missing variables associated with indices $D$ using eq. 31. Please note that the values in the resulting vector are realizations of the independent variables (those in $\vec{z}$ in the previous section). To obtain a sample of the values at dependent nodes, we have to propagate these sampled values associated with independent nodes by means of eq. 14.

Although the computation of the full posterior covariance matrix may be intractable, it is still valuable to be able to compute posterior covariances between quantities of interest. The sparse Cholesky decomposition of the inverse posterior covariance matrix also helps to solve this task. We know that $\left(\mathbf{U}_{\neg D,\neg D}^{\prime}\right)^{-1}=\mathbf{P}^{T}\mathbf{L}^{T}\mathbf{L}\mathbf{P}$ with $\mathbf{P}$ the permutation matrix and $\mathbf{L}$ an upper triagonal matrix. We denote by $\vec{e}_{i}$ a vector with all elements being zero except the element at the i${}^{\textrm{th}}$ position being one. The element $c_{ij}$ in the $i^{\textrm{th}}$ row and $j^{\textrm{th}}$ column of the posterior covariance matrix is given by

$$\vec{e}_{i}^{T}\left(\mathbf{U}_{\neg D,\neg D}^{\prime}\right)\vec{e}_{j}=\vec{e}_{i}^{T}\left((\mathbf{U}_{\neg D,\neg D}^{\prime})^{-1}\right)^{-1}\vec{e}_{j}=\\ \vec{e}_{i}^{T}\left(\mathbf{P}^{T}\mathbf{L}^{T}\mathbf{L}\mathbf{P}\right)^{-1}\vec{e}_{j}=\vec{e}_{i}^{T}\mathbf{P}^{-1}\mathbf{L}^{-1}\left(\mathbf{L}^{T}\right)^{-1}\left(\mathbf{P}^{T}\right)^{-1}\vec{e}_{j}=\\ \left((\mathbf{L}^{T})^{-1}\left(\mathbf{P}^{T}\right)^{-1}\vec{e}_{i}\right)^{T}\left(\mathbf{L}^{T}\right)^{-1}\left(\mathbf{P}^{T}\right)^{-1}\vec{e}_{j}=\\ \left((\mathbf{L}^{T})^{-1}\mathbf{P}\vec{e}_{i}\right)^{T}\left(\mathbf{L}^{T}\right)^{-1}\mathbf{P}\vec{e}_{j}\,,$$

(37)

where we made use of the fact that for permutation matrices $\mathbf{P}^{T}=\mathbf{P}^{-1}$. The multiplications $\left(\mathbf{L}^{T}\right)^{-1}\mathbf{P}\vec{e}_{i}$ and $\left(\mathbf{L}^{T}\right)^{-1}\mathbf{P}\vec{e}_{j}$ can be restated as the solution of the set of linear equations

$$\mathbf{L}^{T}\vec{r}_{i}=\mathbf{P}\vec{e}_{i}\;\;\;\textrm{and}\;\;\;\mathbf{L}^{T}\vec{r}_{j}=\mathbf{P}\vec{e}_{j}$$

(38)

which allows the computation of the covariance element as $c_{ij}=\vec{r}_{i}^{T}\vec{r}_{j}$. Blocks of the posterior covariance matrix can be evaluated by assemblig the (column) vectors $\vec{e}_{i}$ with indices of interest to matrices $\mathbf{E}_{I}=(\vec{e}_{i_{1}},\vec{e}_{i_{2}},\dots)$ and $\mathbf{E}_{J}=(\vec{e}_{j_{1}},\vec{e}_{j_{2}},\dots)$, and then solve the set of linear equations

$$\mathbf{L}^{T}\mathbf{R}_{I}=\mathbf{P}\mathbf{E}_{I}\;\;\;\textrm{and}\;\;\;\mathbf{L}^{T}\mathbf{R}_{J}=\mathbf{P}\mathbf{E}_{J}$$

(39)

to finally compute $\mathbf{C}_{I,J}=(\mathbf{R}_{I})^{T}\mathbf{R}_{J}$. This covariance block is associated with a subset of the independent variables in $\vec{z}_{\neg D}$. To calculate the covariance matrix block associated with the independent variables in $\vec{z}_{D}$, we can use eq. 32. The product in this equation can be computed by solving the set of linear equations

$$\left(\mathbf{U}_{\neg D,\neg D}^{\prime}\right)^{-1}\mathbf{M}=\mathbf{S}_{D,\neg D}$$

(40)

and then computing $\mathbf{S}_{D,\neg D}\mathbf{M}$. Instead of using the index set $D$, it is perfectly possible to use only a subset of D comprising only the indices associated with variables of interest.

Besides computational advantages thanks to exploiting sparsity, the strength of Bayesian networks is their interpretability. Instead of formalizing an estimation problem in nuclear data evaluation in terms of an experimental covariance matrix and parameter prior covariance matrix, hence converting all functional relationships between variables into elements of a covariance matrix, we maintain and work explicity with these relationships. Consequently, we are also able to obtain updated posterior distributions of quantities, which are often marginalized out, such as systematic errors of experiments. This offers a convenient way to spot inconsistencies in the posterior distribution. For instance, a posterior expectation of a normalization error not consistent with its prior specification is a clear warning that some modeling assumptions are wrong.

To emphasize the advantages of the Bayesian network interpretation and also to prepare the discussion of nested relationships in the next section, we give a simple example of how the fitting of model parameters to an experimental dataset with $M$ measurement points affected by statistical and normalization errors may be modeled as a Bayesian network. We then show how to compute the quantities that are required to evaluate the GLS formuals in eqs. 27 and 28.

The relationship of measured values, associated systematic and statistical errors, and the model parameters can be written as

$$\vec{\sigma}_{\textrm{exp}}=g(\vec{p})+\mathbf{J}_{\textrm{norm}}\eta+\mathbf{J}_{\textrm{stat}}\vec{\tau}\,,$$

(41)

with $g(\vec{p})$ assumed to be a linear model, e.g., the linearized version of a nuclear physics model,

$$g(\vec{p})=\vec{y}_{\textrm{ref}}^{\textrm{mod}}+\mathbf{J}_{\textrm{mod}}(\vec{p}-\vec{p}_{\textrm{ref}})\,.$$

(42)

The matrix $\mathbf{J}_{\textrm{mod}}$ denotes the Jacobian matrix of the model. The mapping matrix $\mathbf{J}_{\textrm{norm}}$ of dimension $M\times 1$ maps the normalization error $\eta$ systematicly to each data point. We did not use the same symbol $\mathbf{T}$ as in eq. 10 to denote the Jacobian matrices because the Jacobian matrices here have to be regarded as submatrices of $\mathbf{T}$ and therefore also as submatrices of $\mathbf{S}$ defined in eq. 13 as we will discuss in a moment.

For an absolute normalization error, all elements of $\mathbf{J}_{\textrm{norm}}$ are one. Finally, $\mathbf{J}_{\textrm{stat}}$ is given by the identity matrix to map the statistical errors to all the data points. To be fully consistent with the notation introduced in section II.2, we need to introduce the reference vectors associated with the Taylor approximation,

$$\vec{\sigma}_{\textrm{exp}}=\vec{\sigma}_{\textrm{ref}}^{\textrm{exp}}+\mathbf{J}_{\textrm{mod}}(\vec{p}-\vec{p}_{\textrm{ref}})\\ +\mathbf{J}_{\textrm{norm}}(\eta-\eta_{\textrm{ref}})+\mathbf{J}_{\textrm{stat}}(\vec{\tau}-\vec{\tau}_{\textrm{ref}})\,,$$

(43)

with $\vec{\sigma}_{\textrm{ref}}^{\textrm{exp}}=\vec{y}_{\textrm{ref}}^{\textrm{mod}}+\mathbf{J}_{\textrm{norm}}\eta_{\textrm{ref}}+\mathbf{J}_{\textrm{stat}}\vec{\tau}_{\textrm{ref}}$. The corresponding Bayesian network is given by:

The vectors $\vec{y}$ and $\vec{z}$ are

$$\vec{y}=\begin{pmatrix}\vec{\sigma}_{\textrm{exp}}\\ \vec{p}\\ \eta\end{pmatrix}\;\;\;\textrm{and}\;\;\;\;\vec{z}=\begin{pmatrix}\vec{\tau}\\ \vec{p}\\ \eta\end{pmatrix}\,.$$

(44)

The respective reference vectors for a Taylor expansion are given by

$$\vec{y}_{\textrm{ref}}=\begin{pmatrix}\vec{\sigma}_{\textrm{ref}}^{\textrm{exp}}\\ \vec{p}_{\textrm{ref}}\\ \eta_{\textrm{ref}}\\ \end{pmatrix}\;\;\;\;\textrm{and}\;\;\;\;\vec{z}_{\textrm{ref}}=\begin{pmatrix}\vec{\tau}_{\textrm{ref}}\\ \vec{p}_{\textrm{ref}}\\ \eta_{\textrm{ref}}\\ \end{pmatrix}$$

(45)

and the matrix $\mathbf{S}$ as required for eq. 14 is of the form

$$\mathbf{S}=\begin{pmatrix}\mathbf{J}_{\textrm{stat}}&\mathbf{J}_{\textrm{mod}}&\mathbf{J}_{\textrm{norm}}\\ \mathbf{0}&\mathbb{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbb{1}\end{pmatrix}\,.$$

(46)

Please note that $\mathbf{J}_{\textrm{stat}}=\mathbb{1}$ and the vector $\vec{\tau}$ has the same purpose as in the last section, i.e., to introduce stochasticity into the otherwise deterministic link which connects the model parameters and the normalization error with the experimental measurement vector.

As the values of the experimental measurement vector $\vec{\sigma}_{\textrm{exp}}$ are observed, we can use the GLS update formula in eqs. 27 and 28 to obtain the joint posterior distributions of $\vec{p}$ and $\eta$. To that end, the index set $D$ refers to the indices of the elements associated with the block $\vec{\sigma}_{\textrm{exp}}$ in $\vec{y}$.

If we assume the prior covariance matrix for the model parameter vector $\vec{p}$ to be diagonal, the full prior covariance matrix $\mathbf{U}$ summarizing at the same time our knowledge on model parameters, the normalization error, and the statistical errors is diagonal as well. The mapping matrix $\mathbf{S}$ is very sparse, except possibly for the block $\mathbf{J}_{\textrm{mod}}$ associated with the mapping of model parameters to the corresponding predictions. However, the matrix $\mathbf{J}_{\textrm{mod}}$ is very rectangular for a nuclear model because there are usually much fewer model parameters than predicted values. If using a mathematical fitting function, we can opt for piecewise linear interpolation or similar local interpolation schemes that yield a sparse Jacobian. Prior knowledge on possible shapes of the function, such as its smoothness, can be introduced via the prior specification. Section II.7 elaborates on a sparse Gaussian process construction for this purpose.

In the previous section, we discussed a very simple Bayesian network that directly connected the independent variables to the dependent variables that are measured. However, the full potential of Bayesian networks can only be exploited if we permit nested relationships. Before we discuss inference in the presence of nested relationships in general, we give a motivating example with a nested relationship to develop intuition.

To this end, we extend the discussion of the previous section to incorporate a positivity constraint. We keep the linear model specification in eq. 42, which may represent a piecewise linear model, spline or any other sort of polynomial depending on the form of $\mathbf{J}_{\textrm{mod}}$. Because a linear model may produce negative predictions, we apply a transformation to the model output $\vec{y}_{\textrm{pred}}=f(\vec{y}_{\textrm{mod}})$ with the effect of $f$ being element-wise exponentiation to enforce positivity. This transformation represents a non-linearity in the model description. In the current discussion, we linearize the non-linear relationship and focus on the nested property. In the next section, we show how to deal exactly with the non-linearity.

In vicinity of a vector $\vec{y}_{\textrm{pred}}^{\textrm{ref}}$, the linear Taylor approximation of $f(\vec{y}_{\textrm{mod}})$ can be written as

$$\vec{y}_{\textrm{pred}}=\vec{y}_{\textrm{pred}}^{\textrm{ref}}+\mathbf{J}_{\textrm{pos}}(\vec{y}_{\textrm{mod}}-\vec{y}_{\textrm{mod}}^{\textrm{ref}})\,,$$

(47)

with $\mathbf{J}_{\textrm{pos}}$ being a diagonal matrix with the diagonal elements $J_{\textrm{mod},ii}=\exp(z_{i})$ with $z_{i}$ being the i${}^{\textrm{th}}$ element of $\vec{y}_{\textrm{mod}}^{\textrm{ref}}$. Please note that we did not add a random variable $\vec{\tau}$ to eq. 47 because we deal with a deterministic relationship. For the application of eq. 11 and later equations, such a node must be formally present. Formally, we therefore set all elements in the respective blocks in the reference vector $\vec{z}_{\textrm{ref}}$, the prior estimate $\vec{u}$ and covariance matrix $\mathbf{U}$ to zero, which is equivalent to removing the node altogether.

The Taylor approximation for the compound function can be obtained by plugging the functional form of $g(p)$ into the Taylor approximation of $f$:

$$\vec{y}_{\textrm{pred}}=\vec{y}_{\textrm{pred}}^{\textrm{ref}}+\mathbf{J}_{\textrm{pos}}\big{(}\vec{y}_{\textrm{mod}}^{\textrm{ref}}+\mathbf{J}_{\textrm{mod}}(\vec{p}-\vec{p}_{\textrm{ref}})-\vec{y}_{\textrm{mod}}^{\textrm{ref}}\big{)}=\\ =\vec{y}_{\textrm{pred}}^{\textrm{ref}}+\mathbf{J}_{\textrm{pos}}\mathbf{J}_{\textrm{mod}}(\vec{p}-\vec{p}_{\textrm{ref}})\,.$$

(48)

Finally, the Taylor approximation of this compound function can be combined with the normalization and statistical error, as has been already done in eq. 43,

$$\vec{\sigma}_{\textrm{exp}}=\vec{\sigma}_{\textrm{exp}}^{\textrm{ref}}+\mathbf{J}_{\textrm{pos}}\mathbf{J}_{\textrm{mod}}(\vec{p}-\vec{p}_{\textrm{ref}})\\ +\mathbf{J}_{\textrm{norm}}(\eta-\eta_{\textrm{ref}})+\mathbf{J}_{\textrm{stat}}(\vec{\tau}-\vec{\tau}_{\textrm{ref}})\,,$$

(49)

with $\vec{\sigma}_{\textrm{exp}}^{\textrm{ref}}=\vec{y}_{\textrm{pred}}^{\textrm{ref}}+\mathbf{J}_{\textrm{norm}}\eta_{\textrm{ref}}+\mathbf{J}_{\textrm{stat}}\vec{\tau}_{\textrm{ref}}$. The Bayesian network representing the relationships between the variables in this case is illustrated here:

The quantities to evaluate eqs. 28 and 27 are therefore given by

$$\vec{y}=\begin{pmatrix}\vec{\sigma}_{\textrm{exp}}\\ \vec{y}_{\textrm{mod}}\\ \vec{p}\\ \eta\end{pmatrix}\;\;\;\textrm{and}\;\;\;\vec{z}=\begin{pmatrix}\vec{\tau}\\ \vec{0}\\ \vec{p}\\ \eta\end{pmatrix}\,.$$

(50)

The second element in $\vec{z}$ is zero because we regard the link from the node $\vec{p}$ to $\vec{y}_{\textrm{mod}}$ as deterministic, which is formally solved as stating that the associated parentless node is always zero and also its reference point, prior expectation and associated prior covariance matrix element in $\mathbf{U}$ are zero as well.

The reference vectors of the Taylor expansion are given by

$$\vec{y}_{\textrm{ref}}=\begin{pmatrix}\vec{\sigma}_{\textrm{exp}}^{\textrm{ref}}\\ \vec{y}_{\textrm{mod}}^{\textrm{ref}}\\ \vec{p}_{\textrm{ref}}\\ \eta_{\textrm{ref}}\\ \end{pmatrix}\;\;\;\textrm{and}\;\;\;\vec{z}_{\textrm{ref}}=\begin{pmatrix}\vec{\tau}_{\textrm{ref}}\\ \vec{0}\\ \vec{p}_{\textrm{ref}}\\ \eta_{\textrm{ref}}\end{pmatrix}$$

(51)

and the sensitivity matrix is of the form

$$\mathbf{S}=\begin{pmatrix}\mathbb{1}&\mathbf{0}&\mathbf{J}_{\textrm{pos}}\mathbf{J}_{\textrm{mod}}&\mathbf{J}_{\textrm{norm}}\\ \mathbf{0}&\mathbb{1}&\mathbf{J}_{\textrm{mod}}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbb{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbb{1}\\ \end{pmatrix}\,.$$

(52)

The matrix $\mathbf{S}$ is again very sparse, with the potential exception of the blocks containing $\mathbf{J}_{\textrm{mod}}$.

Because the values of the experimental measurement vector $\vec{\sigma}_{\textrm{exp}}$ are known, the indices in the index set $D$ introduced in the paragraph above eq. 23 refer to the positions of elements of $\vec{\sigma}_{\textrm{exp}}$ in the vector $\vec{y}$. The respective overall Jacobian matrix $\mathbf{S}$ was constructed by hand as an example of how it is done in principle. For more complex Bayesian networks with more deeply nested functions, the manual approach is impractical. In such scenarios, a programmatic recipe to construct $\mathbf{S}$ based on the Jacobian matrices of the individual functional relationships is pertinent.

To this end, we remind ourselves of the notation established in section II.2. There we introduced a vector $\vec{y}$ that comprised the subvector of independent variables $\vec{y}_{I}$ and the dependent variables $\vec{y}_{J}$ whose values are partially determined as a function of the independent variables, $\vec{y}_{J}=f(\vec{y}_{I})+\vec{\tau}$. In the Bayesian network picture, the vector $\vec{\tau}$ contains all the parentless nodes directly attached to the dependent nodes. We also combined the vectors $\vec{y}_{I}$ and $\vec{\tau}$ to the vector $\vec{z}$ that contained the variables associated with all independent nodes. The mental distinction between the parentless nodes in $\vec{y}_{I}$ and those associated with variables in $\vec{\tau}$ was to structurally ensure that we can apply the GLS method by using eq. 27 and eq. 28.

The GLS formula for the posterior expectation eq. 28 relies on the availability of the vector $\vec{y}_{\textrm{ref}}$ containing the propagated values of the reference values $\vec{z}_{\textrm{ref}}$ through the network and the overall Jacobian matrix $\mathbf{S}$. An example of the overall Jacobian matrix for a nested relationship was given in eq. 52.

To describe the necessary propagation and computation of the overall Jacobian, we first introduce the notion of what we call a mapping. Each individual functional relationship between two blocks of variables can be characterized by a set of source indices $\mathcal{S}$ and a set of target indices $\mathcal{T}$ and a function $f$. The two sets must not have indices in common, i.e., $\mathcal{S}\cap\mathcal{T}=\emptyset$.

Each individual mapping takes the elements in $\vec{y}$ indexed by $\mathcal{S}$, applies a function $f(.)$ to it in order to obtain the contribution to the values of the elements in $\vec{y}$ indexed by $\mathcal{T}$. In other words, the application of a mapping on a vector $\vec{y}$ produces a new vector $\vec{y}^{\prime}$ of the same size given by

$$\vec{y}_{\mathcal{T}}^{\prime}=\vec{y}_{\mathcal{T}}+f(\vec{y}_{\mathcal{S}})\;\;\;\;\textrm{and}\;\;\;\;\vec{y}_{\mathcal{\neg T}}^{\prime}=\vec{y}_{\neg T}$$

(53)