Stochastic sampling methodology

SANDY can work as a nuclear data sampling tool applicable to sensitivity analysis (SA) and uncertainty quantification (UQ) problems. The code exploits the basic theory of stochastic sampling to propagate nuclear data covariances through the nuclear models under study.

Definition of best-estimates and covariances

Nuclear parameters $\mathbf{x}=\left[x_1, x_2, \dots, x_m\right]^T$ — e.g., cross sections, resonances, decay constants, particle energy-angular distributions, ... — are stored in the ENDF-6 nuclear data libraries as best-estimate or evaluated data

$$x_{i,BE} = E \left[ x_i \right] = \int x_i p(\mathbf{x}) d\mathbf{x}$$

where $E$ is the expectation operator which performs the integral of the selected $i$-th parameter $x_i$ over the joint probability distribution of $x$, $p(\mathbf{x})$. Best-estimate parameters are often used as inputs of computational models.

UQ studies aim at propagating the model input uncertainties and/or covariances to the model responses. Along this working line, the ENDF-6 format allowed storing the covariance matrices $\mathbf{\Sigma}$ of the best-estimate nuclear data provided in evaluated files.

A covariance matrix $\mathbf{\Sigma}$ of the best-estimate data $\mathbf{x}_{BE}$ contains the so-called variance terms on the diagonal.

$$\Sigma_{i,i} \equiv cov(x_i, x_i) \equiv V(x_i) = E \left[ (x_i - x_{i,BE})^2 \right]$$

From the variance term of a parameter $x_i$ one can extract its standard deviation (or uncertainty)

$$\Delta x_i \equiv \sqrt{V(x_i)} = cov(x_i, x_i)^{1/2}$$

The off-diagonal terms contain the covariances, a generalization of the variance integral to two different variables

$$\Sigma_{i,j} \equiv cov(x_i, x_j) = E \left[ (x_i - x_{i,BE})(x_j - x_{j,BE}) \right]$$

The covariance is a measure of the joint variability of two variables to the same source of error. A covariance $\Sigma_{i,j}$ can be expressed in terms of the standard deviations of its parameters introducing the Pearson correlation coefficient $\rho_{i,j}$:

$$cov(x_i, x_j) = \rho_{i,j} \Delta x_i \Delta x_j$$

The Pearson correlation coefficient is a parameter that varies between -1 and 1 and that can be interpreted as it follows:

• $\rho_{i,j}=1$ implies that $x_i$ linearly increases as $x_j$ increases;
• $\rho_{i,j}=-1$ implies that $x_i$ linearly increases as $x_j$ decreases;
• $\rho_{i,j}=0$ implies that there is no linear correlation between the variables.

A last important concept, which is often used in the rest of the manual, is the relative covariance.

$$rcov(x_i, x_j) = \frac{cov(x_i, x_j)}{x_{i,BE}x_{j,BE}}$$

Stochastic Sampling

It is common to refer with the term stochastic sampling to computational approaches that rely on a repeated random sampling of chosen parameters to obtain statistical outcomes of selected responses.

Monte Carlo sampling methods are based on drawing sets of samples from the model input joint probability density function (PDF), as

$$\mathbf{X}=\begin{bmatrix} x_1^{(1)} & x_1^{(2)} & ... & x_1^{(n)} \\ x_2^{(1)} & x_2^{(2)} & ... & x_2^{(n)} \\ \vdots & \vdots & \ddots & \vdots \\ x_m^{(1)} & x_m^{(2)} & ... & x_m^{(n)} \\ \end{bmatrix}$$

Given input parameters $\mathbf{x}=\left[x_1, x_2, \dots, x_m\right]^T$ with a joint PDF $p(\mathbf{x})$ and covariance matrix $\mathbf{\Sigma}$`, $\mathbf{X}$ is the sample matrix containing on its columns $n$ independent sets of samples of the type $\mathbf{x}^{(k)}=\left[x_1^{(k)}, x_2^{(k)}, \dots, x_m^{(k)}\right]^T$ with $k \in 1,\dots,n$. The samples are taken with replacement, that is, each point in the model input domain is not excluded after being sampled.

One can prove that for any given parameter $i$, when $n\rightarrow\infty$, then

$$\overline{x}i \equiv \frac{1}{n}\sum{k=1}^n x_i^{(k)} \rightarrow E\left[x_i\right]$$

$$s^2_{x_i} \equiv \frac{1}{n}\sum_{k=1}^n \left( x_i^{(k)} - E\left[x_i\right]\right)^2 \rightarrow V(x_i)$$

Analogously, as $n\rightarrow\infty$, the samples' distribution and covariance matrix converge to $p(\mathbf{x})$ and $\mathbf{\Sigma}$, respectively.

Now, let us assume that the user is interested in a model of the type $y =f(\mathbf{x})$, where $y$ is the model response value, or simply response. More specifically, the user wants to quantify the uncertainty of $y$ produced by $\mathbf{x}$ and assuming that the model $f$ does not carry any error.

The average - or best estimate - and uncertainty of the response are simply

$$ y_{BE} = E\left[f(\mathbf{x})\right] = \int f(\mathbf{x}) p(\mathbf{x}) d\mathbf{x}$$

$$ \Delta y = E \left[ (f(\mathbf{x}) - y_{BE})^2 \right] ^{1/2}$$

The Monte Carlo sampling approach for uncertainty quantification uses a straightforward approach to reproduce the integrals above. For any given set $k$ of sampled input parameters $\mathbf{x}^{(k)}$ the model response is quantified as $y^{(k)} = f(\mathbf{x}^{(k)})$.

We often refer to $y^{(k)}$ as to a perturbed response, since it reflects the variation of $y$ generated by the perturbations on $\mathbf{x}$ introduced by the stochastic sampling.

Then, the $n$ perturbed responses can be grouped into a matrix

$$\mathbf{Y}=\begin{bmatrix} y^{(1)} & y^{(2)} & ... & y^{(n)} \\ \end{bmatrix}$$

Notice that this matrix contains only one row, since so far we only considered models with a single response value.

The observable best estimate and uncertainty can be approximated with

$$ \overline{y} \equiv \frac{1}{n}\sum_{k=1}^n y^{(k)}$$

$$ s^2_{y} \equiv \frac{1}{n}\sum_{k=1}^n \left( y^{(k)} - E\left[y\right]\right)^2$$

which converge to the real solutions for $n\rightarrow\infty$.

Should the model response be a vector $\mathbf{y}=\mathbf{F}(\mathbf{x})$, where $\mathbf{y}=\left[ y_1, y_2, \dots, y_l \right]$, then the matrix of perturbed responses becomes

$$ \mathbf{Y}=\begin{bmatrix} y_1^{(1)} & y_1^{(2)} & ... & y_1^{(n)} \\ y_2^{(1)} & y_2^{(2)} & ... & y_2^{(n)} \\ \vdots & \vdots & \ddots & \vdots \\ y_l^{(1)} & y_l^{(2)} & ... & y_l^{(n)} \\ \end{bmatrix}$$

Any two values $y_i$ and $y_j$ could be two completely different model responses, or the same response evaluated in two different points of the output phase-space, e.g. different time, energy or space. In addition, the covariance coefficient of $y_i$ with $y_j$ can be evaluated as

$$ lim_{n\rightarrow\infty} \frac{1}{n}\sum_{k=1}^n \left( y_i^{(k)} - E\left[y_i\right]\right)\left( y_j^{(k)} - E\left[y_j\right]\right) = cov(y_i,y_j)$$

In short, the stochastic sampling for uncertainty propagation simply consists in running the same model multiple times, every time replacing the input parameters with a set of samples. This approach brings numerous advantages compared to other uncertainty propagation techniques, such as perturbation theory:

• its application to any model is straightforward, as it does not interact with the model itself but only with the model inputs;

• there is no need for developing complicated algorithms to calculate adjoint functions for the model and responses under study;

• it does not apply only locally around the input best estimates, but it can cover the whole input and output phase space;

• it can represent any-order effect of the model, as well as interactions between model inputs.

The major drawback of this method lies in the large amount of sets of samples that must be drawn in order to grant the convergence of the response statistics, that is, to reduce the statistical error :math:\epsilon on the response best estimate. The central limit theorem [ref] tells us that this error is proportional to $1/\sqrt{n}$ and $lim_{n\rightarrow\infty} \epsilon = 0$

The huge improvements of computer performances in the last decades, combined with model simplifications and dimensionality reductions help reduce the computational time of the solvers, thus making Monte Carlo sampling a practical option for nuclear data uncertainty propagation.