validation

Monte Carlo Validation

Stochastic evaluation of the P(t) master formula: sampling flux and parameter distributions, estimating the integral, and propagating uncertainty.

Last updated: 2026-07-20

Deterministic quadrature, described in numerical simulation, evaluates the Schubart Master Formula for a single fixed set of inputs. But the inputs to a real evaluation — the flux spectrum, the cross-section, and the material parameters — are themselves uncertain, and a distribution over inputs produces a distribution over outputs. Monte Carlo methods are the natural tool for handling this: they estimate the integral by random sampling and, at the same time, propagate input uncertainty into an uncertainty on $P(t)$ . This chapter describes the methodology generically; the numbers used are illustrative placeholders, never claimed as measured results.

Monte Carlo estimation of the integral

The double integral can be estimated statistically rather than by a fixed grid. One draws $N$ sample points $(E_k,\theta_k)$ from a chosen sampling distribution over the integration domain, evaluates the integrand $f(E,\theta,t)$ at each, and forms a weighted average. For uniform sampling over a domain of measure $V$ , the estimator is

\hat{P}_N(t) = \frac{V}{N}\sum_{k=1}^{N} f(E_k,\theta_k,t),

which is unbiased: its expectation equals the exact integral for any $N$ . What improves with $N$ is not the bias but the variance. If $\text{Var}[f]$ denotes the variance of the integrand over the sampling distribution, the standard error of the estimator is

\text{SE}\big[\hat{P}_N\big] = V\,\sqrt{\frac{\text{Var}[f]}{N}} \;\propto\; \frac{1}{\sqrt{N}} ,

the characteristic $1/\sqrt{N}$ scaling that is the signature — and the cost — of Monte Carlo integration. This convergence is slow: reducing the error by a factor of ten demands a hundredfold increase in samples. Crucially, though, the exponent is independent of the dimension of the domain, which is why Monte Carlo overtakes tensor-grid quadrature when the integrand is high-dimensional or awkwardly shaped, and why it handles sharply peaked domains gracefully.

Variance reduction

Because the error is governed by $\text{Var}[f]/N$ , the same accuracy can be bought either by raising $N$ or by lowering the effective variance of the estimator — and the latter is usually far cheaper. Importance sampling replaces the uniform draw with a distribution $g$ concentrated where the integrand is large, for instance near the peak overlap of $\Phi$ and $\sigma(E)$ , and reweights each sample by $f/g$ . When $g$ is chosen well, the reweighted integrand is nearly constant and its variance collapses, so that a modest $N$ achieves what would otherwise require orders of magnitude more uniform samples. Stratified sampling divides the domain into sub-regions and samples each in proportion to its contribution, guaranteeing coverage of the important regions and eliminating the clustering that plagues purely random draws. Control variates subtract a correlated quantity with known integral, cancelling much of the variance. These are the classical techniques catalogued in the Monte Carlo literature, and at least one of them is usually essential to make an evaluation of the master formula practical rather than merely possible.

Sampling the input distributions

The distinctive use of Monte Carlo in this framework is not merely to evaluate the integral but to propagate uncertainty. Each uncertain input is described by a probability distribution: the flux $\Phi$ might be sampled around a measured mean spectrum with an assumed spread; the cross-section $\sigma(E)$ might be sampled within its measurement error bars; and $\Delta p$ , $v_{\text{ph}}$ , and $\eta_{\text{conv}}$ might each be drawn from distributions reflecting how well they are known. The procedure has an outer and an inner loop. In the outer loop, one draws a complete set of input values from their joint distribution — respecting any correlations between them, since treating correlated inputs as independent can badly misstate the output spread. In the inner loop, the integral is evaluated for that fixed input set, by quadrature or by a nested Monte Carlo integration. Each outer draw thus yields one realisation of $P(t)$ , and repeating the outer loop many times produces an ensemble of outputs whose spread reflects the combined input uncertainty. This nested structure is the essence of forward uncertainty propagation, and it is far more informative than a single deterministic number.

Reading the output distribution

The ensemble of $P(t)$ values is then summarised. Its mean estimates the expected output; its standard deviation quantifies the uncertainty; and its percentiles give confidence intervals. Percentiles are worth reporting alongside the standard deviation because the output distribution is generally not Gaussian — the master formula is a product of factors, and products of positive random variables tend toward skewed, long-tailed distributions for which a symmetric mean-plus-or-minus- $\sigma$ interval understates the upper tail. One might report, purely as an illustration of the method, that a placeholder run yields a mean output with a spread of, say, order ten to twenty percent of the mean, tightening as the input uncertainties are reduced. Such figures are notional and stand only to show the form of the result; they must be re-derived for any real system and must never be quoted as if they were experimental measurements. The width of this output distribution is the quantity most relevant to honest reporting, and its interpretation as a standard deviation connects directly to the statistical discussion in sigma analysis.

Convergence and reliability of the estimate

A Monte Carlo result carries its own error bar, and that error must be shown to be small enough. The standard practice is to monitor the running estimate as $N$ increases and to confirm that both the mean and its statistical error have stabilised; a running-mean plot that has flattened, together with an error bar that has shrunk below the reporting tolerance, is the practical stopping criterion. Because the statistical error falls only as $1/\sqrt{N}$ , halving it requires quadrupling the sample count, so the variance-reduction techniques above are usually essential for efficiency rather than optional. Cross-checking the Monte Carlo mean against a converged deterministic quadrature on the same fixed inputs is a valuable consistency test: the two methods should agree within their combined errors, and a disagreement of more than a couple of standard errors signals a bug in one of them — a mismatched domain, a dropped Jacobian factor, or a sampling distribution that never visits part of the domain. The quality of the pseudo-random stream matters too; a poor generator introduces correlations that violate the independence assumption behind the $1/\sqrt{N}$ law.

What validation does and does not establish

It is important to be precise about what this procedure validates. Monte Carlo validation establishes that the integral has been evaluated correctly and that input uncertainties have been propagated faithfully. It does not, by itself, establish that the input models are physically correct. If the assumed flux spectrum, cross-section, or momentum-transfer mechanism is wrong, the computation will faithfully propagate a wrong assumption to a confidently stated but incorrect output — the familiar “garbage in, garbage out” failure, dressed in a precise-looking error bar. Validation of the numerics is therefore necessary but not sufficient; the physical inputs must be validated separately by measurement. Where a proposed input — for example a speculative $\Delta p$ mechanism — has not been confirmed, the resulting output distribution should be presented as conditional on that assumption, and labelled accordingly. This is also why the statistical spread produced here must never be read as an experimental discovery significance: a tight output distribution certifies that the model is internally consistent given its assumptions, not that the assumptions match the world. The relationship between the statistical spread computed here and the cross-section notation is disentangled in sigma analysis.