Theoretical grounds for the propagation of uncertainties in Monte Carlo particle transport

0.9 MB·

by  Paolo Saracco

1
Theoretical grounds for the propagation of
uncertainties in Monte Carlo particle transport
Paolo Saracco, Maria Grazia Pia, Matej Batic
Abstract—Weintroduceatheoreticalframeworkforthecalcu- at the LHC (Large Hadron Collider). Hardly any effort has
lation of uncertainties affecting observables produced by Monte beeninvestedsofarinestimatingthepredictivecapabilitiesof
4 Carlo particle transport, which derive from uncertainties in
simulation codes, such as EGS5 [22], EGSnrc [23], FLUKA
1 physical parameters input into simulation. The theoretical de-
0 velopments are complemented by a heuristic application, which [24], [25], Geant4 [26], [27], ITS5 [28], MCNP [29], [30]
2 illustrates the method of calculation in a streamlined simulation or PENELOPE [31], commonly used in these experiments:
n environment. this ability would be useful in experimental scenarios where
a Index Terms—Monte Carlo, simulation, uncertainty quantifi- direct validation of simulation use cases would be difficult
J cation or not practically feasible, for instance in some space sci-

6 ence projects, astroparticle physics experiments and medical
1 physics investigations, as well as in the process of detector
I. INTRODUCTION
design,wherethehardwarethatissimulatedmaynotyetexist.
] UNCERTAINTY quantification (UQ) is a fast growing
n This paper defines a theoretical foundation for the calcula-
a sector in interdisciplinary research: its applications span tionoftheuncertaintiesaffectingsimulatedobservables,which
- political science [1], computational biology [2], climate sci-
a areaconsequenceoftheuncertaintiesaffectingtheinputtothe
t ence [3], economic and financial processes [4], industrial and simulation itself. This capability is the basis for establishing
a
civil engineering [5], as well as many other disciplines. In a
d the predictive reliability of Monte Carlo transport codes in
. broad sense uncertainty quantification is a domain of applied experimental practice.
s
mathematics; the variety of its applications has promoted
c The quantification of the uncertainties that affect the re-
i a large number of approaches and methods to address the
s sults of Monte Carlo simulation as a consequence of the
y problem. uncertainties associated with its physical input is a vast and
h Uncertainty quantification is an issue in scientific comput-
complexproblem,whichrequiresextensivescientificresearch.
p ing. Reviews of ongoing research in this field can be found in
This paper is not intended to present an exhaustive solu-
[
[6]–[10]; interested readers can find further information in the tion to the problem, nor to document applications to real-

1 bibliography of the cited references. Software systems, such
life experimental scenarios simulated with general purpose
v as DAKOTA [11], PSUADE [12] and similar codes reviewed
Monte Carlo transport codes; its scope is limited to setting
5 in [9], have been developed to facilitate this task, mainly
4 a theoretical ground, which to the best of our knowledge has
focusing on methods and algorithms for sensitivity analysis
0 never been previously documented in the literature, to enable
and statistical evaluation of uncertainties.
4 further conceptual and mathematical progress in this field in
. Uncertainty quantification is especially relevant to physics
1 view of future experimental applications.
simulation, where the ability to estimate the reliability of
0 A preliminary report of this study is documented in [32].
4 simulated results is critical to establish it as a predictive
1 instrument for experimental research. Nevertheless, relatively
: limited attention has been invested so far into the problem of II. OVERVIEWOFTHEPROBLEMDOMAIN
v
i quantifying the uncertainties of the outcome of Monte Carlo Uncertaintyquantificationinthecontextofacomputational
X
particle transport in general terms.
system is the process of identifying, characterizing, and quan-
r Investigations of uncertainty quantification in the domain
a tifying those factors that could affect the accuracy of the
of Monte Carlo particle transport mainly concern applications
computational results [10].
tonuclearpowersystems,suchas[13]–[16].Commonexperi-
Uncertainties can arise from many sources; the computa-
mentalpracticeinotherapplicationareas,suchashighenergy
tional model propagates them into uncertainties in the results.
physics experiments, focuses on the validation of specific
This problem is usually referred to as forward uncertainty
use cases by direct comparison of simulation results and
quantification. In experimental practice one encounters also
experimental measurements: representative examples of this
the problem of backward uncertainty quantification, i.e. the
practicecanbefoundin[17]–[21],whichconcernexperiments
assessmentoftheuncertaintiesthatmaybepresentinamodel:
this issue is of raising interest for its applicability in robust
ManuscriptreceivedJune15,2013.
P. Saracco and M. G. Pia are with INFN Sezione di Genova, Genova, design, as it concerns the ability of making a rational choice
16146 Italy (telephone: [phone], e-mail: [email protected], among different conceptual designs that can be drafted.
[email protected]).
Uncertainties in the results can derive from the conceptual
M.BaticwaswithINFN,Genova,16146ItalyandnowiswithSinergise,
Ljubljana,Slovenia. model upon which a computational system is constructed,
2
the formulation of the model in the software and the actual models: theoretical calculations that describe particle interac-
computation process. Possible sources of uncertainties are tions with matter are usually tabulated in data libraries, or
characterized in [33]: transformed into look-up tables in the initialization phase of
1) parameter uncertainty identifies situations where some the Monte Carlo simulation to reduce the computational bur-
of the computer code inputs are unknown; den of the simulation. These data implicitly contain modeling
2) model inadequacy may derive from structural uncer- assumptions and approximations; the associated uncertainties,
taintyrelated,forinstance,toapproximationsintheused although treated as parameter uncertainties, may encompass a
model, and from algorithmic uncertainty related to the component of model inadequacy.
numerical methods employed to solve the model; The problem domain encompasses correlated and uncorre-
3) residual variability occurs when the process itself may lated uncertainties. In the context of Monte Carlo simulation
be inherently unpredictable and stochastic, or the model for particle transport, correlated uncertainties may originate
itself is not fully specified; from systematic errors in the physical data required for the
4) parametric variability concerns use cases where some simulation, while uncorrelated ones may be associated with
of the inputs are intentionally uncontrolled or unspec- statistical errors of input experimental parameters. Mixed and
ified, therefore they contribute a further component of more complex situations may also occur, for instance where
uncertainty to the predicted process. experimental data with their associated uncertainties are used
Parameter uncertainty plays a major role in Monte Carlo in connection to physical models, whose validity is in turn
simulations of particle transport, since in practice all the variable.
physicalinputtothesimulationisaffectedbyuncertainties:the
cross sections determining particle interactions either derive
III. STRATEGYOFTHISSTUDY
from interpolations of experimental measurements, which are
affected by uncertainties, or from analytical models, which This paper deals with forward propagation of uncertainties
in turn involve experimental or theoretical uncertainties. In in Monte Carlo simulation.
addition, parameter uncertainties are involved in the model The conceptual abstractions involved in the process of
of the experimental set-up: they concern the geometrical uncertainty propagation and their functional relationships are
dimensions and material composition of the apparatus, and illustrated in Fig. 1. Input parameters of the simulation are
the conditions of its operation, such as electromagnetic fields, affected by uncertainties: conceptually, this situation corre-
pressure and temperature. sponds to executing simulations with many possible input
Algorithmic uncertainty, as defined by [33], derives from parameters,whichvarycompatiblewiththeirassociatedprob-
the numerical methods employed to solve the model: in this ability distributions, and produce many possible outputs. The
acception it is implicit in the Monte Carlo transport process, statistical properties of the outcome of the simulations repre-
whichisastatisticaldeterminationofsomephysicalquantities sent the effect of input uncertainties.
(densities, fluxes, energy deposition etc.) of experimental Thefollowingsectionelaboratesamathematicalfoundation
interest through sampling methods. forthecalculationofuncertaintiesinobservablesproducedby
Structural uncertainties may be present in the physics aMonteCarlosimulation,whichderivefromtheuncertainties
models on which the transport is based: for instance, some affecting quantities fed into the particle transport process.
effects may be neglected a priori in the formulation of the The adopted approach aims at a complete determination of
models(e.g.thosederivingfromthemolecularcompositionof the probability distribution function (PDF) of the output ob-
materials, or their solid structure), or the models may involve servable; in this respect, it differs from other approaches to
assumptions that are questionable in some energy range of uncertaintyquantification,suchasthosepursuedbyDAKOTA,
the transported particles. Condensed history methods fall into PSUADE and similar software tools, which focus on the
this category: they are approximate Monte Carlo methods to management of large ensembles of calculations used in a
deal with physics processes affected by infrared divergence statistical evaluation of uncertainty.
and involving a large number of collisions with very small Although the theoretical foundations and calculation meth-
changes in direction and in energy, which would require a ods developed in this paper are generally applicable to any
prohibitive investment of computing resources to be treated inputparameters,thefollowingdiscussionsaremainlyfocused
exactly. Moreover, the concept itself of particle transport is on the quantification of the effects of physical uncertainties,
an approximation to the real world: Monte Carlo codes for and their interplay with the algorithmic uncertainties associ-
particle transport are based on the assumption of classical ated with the Monte Carlo sampling process.
particles moving between localized interaction points, while Theoretical investigations are carried out in a simplified
quantum mechanics is implicit in the determination of the calculationenvironment,whichretainstheessentialconceptual
cross sections involved in the process. features characterizing the problem domain. The scheme for
Whereas these concepts provide guidance in the identifi- uncertainty propagation is initially developed for quantifying
cation of possible sources of uncertainties, a rigid scheme the effects on an observable associated with the uncertainties
of classifications can hardly capture the complexity of the affectingasingleparameter.Theconceptualfoundationestab-
problem domain. For instance, parameter and structural un- lished in the present paper is the basis for further calculations
certainties are often intertwined in Monte Carlo codes for in more complex environments, representing more realistic
particle transport, as input data often actually embed physics experimental scenarios, which are intended to follow.
3
The theoretical elaboration described in section IV estab-
Input physical data
lishes that it is in principle possible to disentangle the effect (with uncertainties)
ofparameteruncertaintiesfromalgorithmicuncertainties:itis
possibletodetermineaprobabilitydistributionoftheoutcome
Physical data sample 1 Physical data sample 2 … Physical data sample N
of the original physical problem deriving from parameter MC
uncertainties alone. By algorithmic uncertainties we mean
…
those related only to the numerical methods employed to Simulation 1 Simulation 2 Simulation NMC
solve the underlying transport equations, i.e. from Monte
Carlo sampling, as defined in section II. The effect of Monte
Observable distribution 1 Observable distribution 2 … Observable distribution N
Carlo sampling consists of blurring the distribution resulting MC
from the propagation of parameter uncertainties with some
statistical noise.
Statistical analysis
The next stage of investigation consists of devising a
Observable ± error
procedure, based on these findings, to calculate in practice
the probability distribution of the possible outcomes of an Fig.1. Conceptualschemeoftheuncertaintyquantificationprocess.
observable within a desired margin of statistical error. The
possibility of performing this calculation by means of a small
number of Monte Carlo simulations is examined in section V, A. Propagation of parameter uncertainties
while other practical aspects of uncertainty quantification are
If one runs many Monte Carlo simulations, each one en-
discussed in section VII.
compassing N events, varying one of the physical parameters
The discussion of the theoretical aspects of the problem
involved in particle transport (identified here with Σ) in some
is supported by a heuristic investigation, documented in sec-
interval with probability distribution f(Σ), then the final
tionVI,whichillustratesempiricallytheprocessofuncertainty
distribution G(x) for the value x of any desired observable
propagation and the conclusions deriving from its mathemat-
X, which depends on the value of Σ, is expected to be
ical foundation.
In the context of this paper it is assumed that uncertainties
G(x)(cid:39) (1)
affecting the physical input to Monte Carlo particle transport (cid:34) (cid:35)(cid:115)
are known. Their quantitative knowledge is established by (cid:90) +∞ dΣf(Σ)exp −(x−x0(Σ))2 N
evaluationproceduresappliedtodatalibrariesofexperimental 2σ2 /N 2πσ2
−∞ x0 x0
origin, such as some components of [34], and dedicated
validation tests, such as those documented in [35]–[41]. This where we finally integrated over the possible values of the
assumption does not reflect the real status of general purpose physical parameter Σ with their probabilities. Here x is a
Monte Carlo codes for particle transport: although some of stochasticvariablebecauseitrepresentsthesampledstatistical
them have been widely used for decades, the validation of meanofthecontributionstotherequestedobservablefromthe
their physics modeling is still in progress, and quantitative N events and σ2 is its sampled variance. This result derives
x0
assessments of the uncertainties associated with their basic from the Central Limit Theorem, if N is sufficiently large to
physics parameters are scarcely documented in the literature. make σ2 independent from N itself.
x0
This process is further hindered by the presence of epistemic When writing this equation, we made a further assumption:
uncertainties, i.e. intrinsic knowledge gaps [42]. The current namely that a function x (Σ) exists, which relates any input
0
incomplete quantification of the physical input uncertainties value of Σ to a corresponding value for the peak position of
in Monte Carlo particle transport codes does not hinder the the observable means distribution.
conceptual foundation of the theoretical scheme discussed This assumption relies on the interpretation of the process
in this paper. Nevertheless, such knowledge is required for of simulation as a surrogate for the solution of the Boltz-
practical applicability of the findings of this paper, i.e. for the mann transport equation. The function x (Σ) represents the
0
calculation of uncertainties affecting output observables. parametric dependence of the physical solution we produce
through the simulation to a variation of the input parameter.
Fullknowledgeofx (Σ)isequivalenttotheabilityofsolving
0
IV. THEORY the transport equation.
If this function is invertible - a fact that, as we shall see
later in Sect.VI-F, can be established when necessary - we
The subject under study is how the sources of uncertainty
may change the variable of integration from Σ to x ; using
in the input to the simulation, wherever they originate, are 0
Σ=Σ(x ), we obtain in the limit N →∞
mathematically mapped to uncertainties in the simulation re- 0
sults. We show that under some general conditions, which are
(cid:90) +∞
elucidated in the following, Monte Carlo simulation transfers G(x) = dΣf(Σ)δ(x−x (Σ)) (2)
0
thedistributionofinputdatauncertaintiesintothedistribution −∞
(cid:12) (cid:12)
of the results, preserving the same functional form or with a (cid:12)dΣ(x)(cid:12)
= (cid:12) (cid:12)f(Σ(x))
form derived from it, and adding some statistical fluctuation. (cid:12) dx (cid:12)
4
0.5
where we used the identity
(cid:72)a(cid:76)
(cid:114) 1 (cid:90) +∞ (cid:34) (x−x )2(cid:35) 0.4 Σ(cid:60)A(cid:62)
f(x)= lim dx f(x )exp − 0
σ→0 2πσ2 −∞ 0 0 2σ2
0.3
stemmingfromtheexponentialrepresentationofDiracδ func- (cid:76)A
tion. If the function x0(Σ) is invertible only over subintervals (cid:72)G 0.2
of the variability interval for Σ - that is the same physical
solutionx0 resultsfromdifferentvaluesoftheinputparameter 0.1
Σ - we must slightly modify (2) by summing over all the
possible solutions of the equation Σ=Σ(x0). 0.0
Equation(2)statesthatweareabletoexactlyknowhowthe 2.0 2.5 3.0 3.5 4.0 4.5 5.0
A(cid:72)arbitraryunits(cid:76)
inputprobabilitydistributionofphysicaldatatransfersintothe
final result, provided we are able to determine the unknown
0.5
function Σ(x). This conclusion is the theoretical foundation (cid:72)b(cid:76)
for uncertainty quantification.
0.4
We derived equation (2) in the context of Monte Carlo
simulation; nevertheless, through the limit N → ∞ the 0.3
formula has dropped any direct reference to the Monte Carlo (cid:76)A
(cid:72)G
simulation environment and expresses a general property of 0.2
forward propagation of uncertainty for an arbitrary determin-
istic problem. In other environments this expression is known 0.1
as the Markov formula [43], [44].
The theory here presented is valid independently from the 0.0
numberofunknowns:thefirstlineof(2)holdsevenifΣ −→ 2.0 2.5 3.0 3.5 4.0 4.5 5.0
A(cid:72)arbitraryunits(cid:76)
{Σ ,...,Σ } and x −→ {x ,...,x }.
1 M 1 K
Fig.2. TheexpectedPDFG(A)fortheareaofacircleasafunctionofA
B. Verification withrandomlychosenradii0.8≤R≤1.2(thickblackcurve),togetherwith
theempiricaldistributions(grey)comingfrom100000MonteCarloestimates
The correctness of the method outlined in section IV-A can ofthesamearea.In(a)theMonteCarlosimulationencompasses500samples,
be verified in a scenario where the outcome of a Monte Carlo in (b) 50000. In (a) the standard deviation on < A > at the fixed value
calculation is an exactly solvable problem. A=Amax isalsoshownwithtwodashedlines
As an example, we consider the calculation of the area of a
circle; we assume that the radius of the circle is known with
some uncertainty. The solution A(R) = πR2 is invertible as visible,andappearsmoreevidentintheplotresultingfromthe
R(A) = (cid:112)A/π for R > 0 (the only physically acceptable larger Monte Carlo sample; it is quantitatively confirmed with
solution), where A is the area and R is the radius. 0.05 significance by a statistical comparison of the results:
Onthebasisofequation(2),theexpectedprobabilitydistri- the p-values are respectively 0.083 and 0.910 for a Cramer-
bution function for the area of a circle with input uncertainty Von-Mises test with a bin width of 0.01. The blurring effect
on the measure of its radius is: deriving from the process of Monte Carlo calculation is also
manifest.
dR(A)
G(A)= f(R(A)) (3)
This verification test is also instrumental to highlight an
dA
issue,whichwillbediscussedmoreextensivelyinsectionVI:
where f(R) is the PDF for the unknown input radius. If, for
the interplay of different statistical errors in the process of
example, f(R) is a flat distribution with R ≤R ≤R ,
min max Monte Carlo sampling.
then
(cid:114) The upper plot of Fig. 2 also shows the standard deviation
1 π
G(A)= 2π Aθ(A−Amin)θ(Amax−A) (4) on < A > at the fixed value A = Amax with two dashed
lines.Thisshouldhelpthereadertounderstandthatwithinthe
We consider a test case where the radius is uniformly process of Monte Carlo sampling the exact result expressed
distributed in the interval 0.8 ≤ R ≤ 1.2 (in arbitrary units), in equation (2) implies the mixing of two different statistical
so that A = 2.011 and A = 4.524 (in arbitrary units). errors: one deriving from the finiteness of the number of
min max
The results of the Monte Carlo calculation of the area are events N in each MC simulation (see equation (1), and one
compared with the theoretical prediction of equation (4) of its stemming from the finiteness of the number of Monte Carlo
PDFinFig.2:thetwohistogramscorrespondtotheempirical really run. The effect of the former is clearly visible on the
distribution of the area resulting from 500 and 50000 Monte sides of the empirical distribution, which are exponentially
Carlosamplesrespectively,eachsampleencompassing100000 smearedaroundtheexactvaluesA andA byanamount
√ min max
estimates,whilethethickcurverepresentsthetheoreticalPDF σ = σ / N. The comparison of Fig. 2(a) and 2(b)
<A> A
asin(4).TheconsistencyofthePDFofthesimulatedobserv- makesalsovisuallyevidentthelimitingprocessingoingfrom
able with the theoretical prediction (A−1/2) is qualitatively (1)to(2):in(b)σ isreducedbyafactor10with(therefore
<A>
5
it cannot be visually displayed in the plot), because we used (see Fig. 9). In this case equation (2) reduces to
100 times more Monte Carlo samples. The oscillations of the
G(x)∝f(Σ(x)) , (5)
empirical distribution around its theoretical value within the
interval of allowed values Amin < A < Amax are, instead, whichrepresentsthePDFfortheobservablemeansasalinear
a consequence of the finiteness of the number NMC of the mapofthePDFfortheunknowninputphysicalparameter.The
numberofMonteCarlosimulationsreallyrun:theiramplitude twoparameterscharacterizingthelinearmapcanbepractically
is ruled by 1/N1/2. determined with any predefined margin of statistical accuracy
MC
The theoretical elaboration in section IV-A has reduced by running two Monte Carlo simulations for two distinct
the problem of uncertainty quantification from the computa- values of the input variable Σ.
tionally demanding production of the empirical distribution, Whileitisunlikelythatx (Σ)islinearingeneral,anyfunc-
0
requiring a large number of Monte Carlo simulation runs (as tion can be approximated by a linear relation in a sufficiently
depicted in Fig. 1), to the determination of the parameters small interval [45]. For the practical purpose of reducing
fixing the form of the expected theoretical distribution: in this equation (2) to the easily manageable form of equation (5),
verification test they are Amin, Amax and the exponent of its it is sufficient that x0(Σ) can be considered linear over an
powerform.Inthistestcasetheyareallknown,sincethegiven interval relevant to the simulation scenario: over the whole
problem is analytically solvable; in a realistic situation these variability interval, if the input distribution is supported on
parameters must be determined through the simulation itself. a bound interval, otherwise over some predefined confidence
In section VI-D we will show that it is possible to determ√ine interval.
theseparameterswithanaccuracyruledbyσ<x0> =σx0/ N The hypothesis of a linear approximation of the functional
through a small number of Monte Carlo runs. in this example form of Σ(x) in the experimental scenario under study can
the observable x0 coincides with the area A. be verified by running Monte Carlo simulations for a few
suitably chosen points in the variability interval for Σ. As
a minimal requirement for estimating the suitability of a
V. STATISTICALESTIMATION
linearapproximationonecouldenvisagethreesimulations,for
Equation (2) is fundamental to the purpose of uncertainty
instance at the two interval extrema and at the mode, while
quantification, as it relates the functional form of the output
a more powerful test of the hypothesis of linearity could be
PDF G(x) to those of f(Σ), which is a necessary input, and
performed with a larger number of simulation runs.
of Σ(x ) or x (Σ). Nevertheless, equation (2) is of limited
0 0 Thesameconceptualprocedureappliesif,insteadofalinear
practicalusefulnesstoquantifytheuncertaintiesofobservables
relationship, one considers other computationally manageable
resulting from Monte Carlo simulation, since it is equivalent
approximations, for instance a power representation of the
to solving the transport equation: if one were able to do this,
form
one would not need any Monte Carlo simulation at all. (cid:88)m
x (Σ)= a (Σ−<Σ>)k.
A remark is here in order: equation (2) is an exact relation, 0 k
even if we derived it starting from (1) which is valid only for k=0
a MC simulation; knowledge of x (Σ) enables to know the Theproceduredescribedabovetostatisticallydeterminethe
0
required form of G(x) from which any statistical information parameters of the theoretical probability distribution function
about the observable x can be extracted. Then in principle of an observable is computationally less expensive than deter-
whicheversolverforthetransportequationcanbeusedtothis mining the probability distribution function of the observable
purpose, not only MC simulation. Most part of the discussion from a large number of Monte Carlo simulations, each one
in this Section and in Sect.VII can be translated accordingly, performed with a different value of the input parameter con-
by substituting ”statistical accuracy” with ”numerical accu- sistent with its uncertainty.
racy” and ”MC simulation” with ”numerical solver for the
transport equation”. VI. APPLICATIONOFTHEMETHOD
A practical approach to address the problem consists of
We illustrate here an application of the method for un-
using the Monte Carlo simulation process itself to statistically
certainty quantification elaborated in the previous sections.
estimatetheunknownfunctionΣ(x ):thisispossible,because
0 The calculation is performed in a streamlined simulation
anysingleMonteCarlorun,executedwithafixedvalueofΣ,
environment, which is identified in the following as “toy
gives a statistical estimate of the value of x (Σ). Uncertainty
0 MonteCarlo”:althoughsimplified,thiscomputationalscenario
quantification requires the ability to perform this estimate
implementsalltheessentialelementsoftheconceptualframe-
withinsomepredefinedmarginoferror.Itisworthnotingthat
work of uncertainty quantification illustrated in Fig. 3.
the task discussed here is not the prediction of the simulated
This example of application serves also as an environment
PDF for x resulting from a large number of Monte Carlo
forheuristicanalysis,withthepurposeofidentifyingmethods
simulations, which would be computationally expensive, but
of more general applicability.
only the calculation of the parameters determining the final
exact PDF G(x) defined in (2).
A. Toy Monte Carlo
This task becomes particularly simple when the relation
x (Σ)islinear,orcanbetreatedasapproximatelylinear,asit The toy Monte Carlo is a lightweight simulation system,
0
happens for instance in the example we give in the following interfaced to the computational environment of Mathematica
6
[46]forfurtherelaborationofitsresults.Itisaninstrumentfor
verification and heuristic investigation of theoretical features
0.10
of Monte Carlo particle transport relevant to the problem
of uncertainty quantification; it is not intended to implement val 0.08
physics and experimental modeling functionality for real-life r
e
simulation application scenarios. nt
i 0.06
n
It consists of a random path generator, which is ruled by i
s
two constant parameters, describing the relative probability of ore 0.04
absorption (Σ ) and scattering processes (Σ ) a transported c
A S s
particlemayundergo.Theuserisfreetomodifytheseparam- (cid:37)
0.02
eters. This simplified situation physically corresponds to the
transport of neutral particles through a uniform medium, with
0.00
energy independent absorption and scattering cross sections, 3.69 3.73 3.77 3.81
the scattering occurring always in S-wave (i.e. scattering an- O x102(cid:72)arbitraryunits(cid:76)
tl
glesinthecenter-of-massframearesampledfromanuniform
distribution). Fig.4. Resultsfor1000ofMonteCarlosimulationsfortheobservableOtl,
eachencompassing106events,foranobservablescoredclosetotheprimary
Within this streamlined simulation context the units in
particle source (see text), produced with different values of the ΣS input
which lengths and macroscopic cross sections are expressed physicalparameter:ΣS =1(whitehistogram),ΣS =1.1(greyhistogram)
arenotrelevant;thereforeforsimplicityweusearbitraryunits andΣS =1.2(blackhistogram).
for these two quantities, with the only constraint of being one
the reciprocal of the other.
The experimental configuration consists of an isotropic
primary particle source located at the origin of a sphere of
uniform material. The radius of the sphere is 5 in arbitrary
units.
In the application example considered here we focus on
the track length as an observable (O in the following) of
tl
interest in the simulation; nevertheless the issues discussed in
the following are not specific to this observable, rather they
highlight general features of the process of estimating the un-
Fig. 5. Results for 1000 of Monte Carlo simulations for the observable
Otl,eachencompassing106 events,foranobservablescoredfarawayfrom
the primary particle source (see text), produced with different values of the
ΣS input physical parameter: ΣS = 1 (white histogram), ΣS = 1.1 (grey
histogram)andΣS =1.2(blackhistogram).
certainty associated with a result produced by the simulation.
The observable is scored in two sensitive volumes, consisting
ofacubeofside0.2inarbitraryunits:onewithcenterlocated
in (0.1,0.1,0.1), and one centered in (0.7,0.7,0.7).
The execution of the Monte Carlo simulation contributes to
theresultthroughthestatisticalsamplingprocess.Ifthephys-
ical parameters Σ and Σ are affected by uncertainties, we
A S
expect the simulation outcome to be sensitive to the effects of
parameter uncertainties mixed with algorithmic uncertainties.
We stress that this toy Monte Carlo is only used to verify
andtoexemplifyresultsthathavebeenderivedonatheoretical
groundinsectionIV:itsusedoesnotconstitutealimitationto
Fig.3. Thegeometryofthe”toyMonteCarlo”.Thesourceoftransported
particlesisatthecentre.Theenclosingspherehasbeenscaledtomakevisible the validity of the theoretical assessments established in this
theactivevolumesforthetracklengthobservableused. paper,ratheritrepresentsanauxiliarytooltoourinvestigation.
7
0.08 0.08
al al
v v
nter 0.06 (cid:72)a(cid:76) nter 0.06 (cid:72)b(cid:76)
i i
n n
i i
s 0.04 s 0.04
e e
r r
o o
c c
s 0.02 s 0.02
(cid:37) (cid:37)
0.00 0.00
3.54 3.68 3.82 3.96 3.54 3.68 3.82 3.96
O x102(cid:72)arbitraryunits(cid:76) O x102(cid:72)arbitraryunits(cid:76)
tl tl
0.08 0.08
al al
v v
nter 0.06 (cid:72)c(cid:76) nter 0.06 (cid:72)d(cid:76)
i i
n n
i i
s 0.04 s 0.04
e e
r r
o o
c c
s 0.02 s 0.02
(cid:37) (cid:37)
0.00 0.00
3.54 3.68 3.82 3.96 3.54 3.68 3.82 3.96
O x102(cid:72)arbitraryunits(cid:76) O x102(cid:72)arbitraryunits(cid:76)
tl tl
Fig.6. EmpiricaldistributionsofthetracklengthobservableOtl scoredclosetotheprimaryparticlesourceobtainedby105 MCsimulations.Eachsingle
run has been produced with a random value for ΣS in the interval 1.≤ΣS ≤1.2; each simulation consists of (a) 105 events, (b) 5105 events, (c) 106
events,(d)108 events.
B. Algorithmic uncertainties source mentioned in section VI-A, respectively. For graphical
reasons here and in the following we make use of the unusual
The events generated in the course of a Monte Carlo notationO ×10m,eveniftheobservablescoredisinarbitrary
tl
simulationcontributetodetermineobservablesofexperimental units: this notation has the purpose to indicate that the same
interest. Each particle trajectory encompassed in an event units (even if arbitrary) have been used in all the plots; the
is the result of a sequence of random interactions sampled indicatedmultipliermeansthatinFig.4thex-scalestartsfrom
from some statistics; the observable score from single events 0.0369 while in Fig. 5 starts from 0.001.
itself is a stochastic variable with some unknown distribution. If the number of generated events is sufficiently large, the
The outcome of the Monte Carlo simulation is the mean Central Limit Theorem ensures that, whatever is the unknown
of the contributions from single events. Nevertheless, we observabledistribution,thedistributionoftheobservablemean
cannotdirectlydeterminethemean,thevariance,thestandard estimatesproducedbyrunningmanyMonteCarlosimulations
deviation or any other properties of the distribution of the is approximately gaussian, with standard deviation
observable exactly, because we do not know the distribution √
σ =σ / N (6)
itself:wecanonlyestimatethesequantitiesfromthestatistical <Otl> Otl
sample.Asaconsequence,theobservablemeanestimateitself The two experimental configurations, characterized by dif-
is a stochastic variable. ferent distances between the source and the detector, give rise
As a first example, we illustrate the effect of algorithmic todifferentstatisticalerrors,becausealargernumberofevents
uncertainties, i.e. of uncertainties purely due to the Monte contribute to the observable score, if the sensitive volume is
Carlosamplingprocess.Weconsiderthreephysicalscenarios, closer to the source.
where Σ assumes the values 1.0, 1.1 and 1.2 respectively
S
(in arbitrary units), while Σ is fixed and has value 0.1 (in C. Effect of cross section uncertainties
A
arbitrary units). In this use case the physical parameters Σ Within the context of the toy Monte Carlo, we consider
S
andΣ areassumedtobeexemptfromuncertainties.Foreach a scenario where one of the input parameters is affected by
A
physical scenario we execute 1000 Monte Carlo runs, each uncertainties:itsvaluemayvaryoversomeintervalwithsome
one encompassing 106 events corresponding to the transport probability distribution. For simplicity we assume it to be
of a single primary particle. The resulting distributions of the a uniform distribution, which is corresponds to the typical
observable, corresponding to the three different values of Σ scenario of an epistemic uncertainty.
S
assumed in the simulation, are plotted in Fig. 4 and Fig. 5 In this application example Σ is sampled from a flat
S
for the sensitive volume at closer and larger distance from the distribution in the interval 1. ≤ Σ ≤ 1.2. We perform 104
S
8
0.06 0.06
al 0.05 (cid:72)a(cid:76) al 0.05 (cid:72)b(cid:76)
v v
r r
nte 0.04 nte 0.04
i i
n n
i 0.03 i 0.03
s s
e e
or 0.02 or 0.02
c c
s s
(cid:37) 0.01 (cid:37) 0.01
0.00 0.00
0.976 1.02 1.06 1.10 0.976 1.02 1.06 1.10
O x103(cid:72)arbitraryunits(cid:76) O x103(cid:72)arbitraryunits(cid:76)
tl tl
0.06 0.06
al 0.05 (cid:72)c(cid:76) al 0.05 (cid:72)d(cid:76)
v v
r r
nte 0.04 nte 0.04
i i
n n
i 0.03 i 0.03
s s
e e
or 0.02 or 0.02
c c
s s
(cid:37) 0.01 (cid:37) 0.01
0.00 0.00
0.976 1.02 1.06 1.10 0.976 1.02 1.06 1.10
O x103(cid:72)arbitraryunits(cid:76) O x103(cid:72)arbitraryunits(cid:76)
tl tl
Fig.7. EmpiricaldistributionsofthetracklengthobservableOtl scoredfromthesource(seetext)obtainedby105 MCsimulations.Eachsinglerunhas
beenproducedwitharandomvalueforΣS intheinterval1.≤ΣS ≤1.2;eachsimulationconsistsof(a)105 events,(b)5105 events,(c)106 events,(d)
108 events.
Monte Carlo simulation runs, each one with a different value uncertainties of the parameters on which they depend. These
of Σ ; each simulation encompasses a predefined number empiricalobservationsareconsistentwiththetheoreticalfoun-
S
N of events. As described in section VI-A, the track length dations discussed in section IV.
observableisscoredintwovolumesatdifferentdistancesfrom
the primary particle source.
D. Probability distribution function of an observable
Fig. 6 shows the distribution of the estimated observable
Equation (2) indicates that the task of uncertainty quan-
means close to the source resulting from the simulation; the
tification is reduced to the determination of the output PDF
four histograms correspond to different number N of events
G(x). From this distribution one can extract all the necessary
in each simulation run. For a relatively small number of
information to quantify the knowledge of the output physical
generated events, illustrated in the top row of histograms
variable x produced by the simulation, such as the desired
in Fig. 6, the distribution of the observable mean estimates
confidence intervals. In the following we show how one can
resembles a gaussian distribution. By increasing the number
determine the characteristics of the PDF of the simulated
ofgeneratedeventsfromFig.6(a)toFig.6(b),thewidthofthe
√ observable with predefined accuracy in the simple application
distributiondecreases(roughlywith N).Forlargervaluesof
scenariointroducedinsectionVI-C.Todothisweexaminethe
N thewidthstabilizestoaconstantvalueandthehistogramno
propertiesofthedistributionfromwhichsimulationresultsare
longer resembles a gaussian distribution: eventually it appears
actually sampled - eq. (1) - that we indicate here as G (x),
consistent with a flat distribution. emp
to distinguish it from the exact result G(x).
Fig. 7, corresponding to the observable scored at larger In the case of a flat distribution of the input parameter
distance from the primary particle source, exhibits the same f(Σ ) = 1/(Σ − Σ ) over an interval Σ ≤
S Smax Smin Smin
qualitative behavior, although the process of approaching a Σ ≤ Σ , the integral in equation (1) can be evaluated
S Smax
final stable distribution appears slower in this case. This is a analytically even for finite N, if we assume that σ = σ
consequence, as noted in section VI-B, of the larger statistical is independent from Σ and that Σ (x) is linear. Txh0e first
S S
error associated with sensitive volumes far from the source. assumption is justified, if one considers results at the lowest
√
This example shows qualitatively that, apart from some order in 1/ N [47]; we recall that σ is the standard
x0
additional statistical noise, the distribution of observables deviationoftherequiredobservable(nottheoneofitsmean),
resulting from the Monte Carlo simulation process retains a as it results from the simulation. The second assumption is
functional form that is related to the input distribution of the discussed in detail in section VI-F. Under these two condi-
9
tions the empirical distribution coming from the Monte Carlo
simulation resulting from (1) is:
(cid:34) (cid:32)√ (cid:33)
dΣ /dx N(x−a)
G (x) (cid:39) S erf √
emp 2(ΣSmax−ΣSmin) 2σ
(cid:32)√ (cid:33)(cid:35)
N(x−b)
−erf √
2σ
being
2 (cid:90) x
erf(x)= √ exp[−t2]dt;
π
−∞
or
(cid:34) (cid:32)√ (cid:33)
1 N(x−a)
G (x) (cid:39) erf √
emp 2(b−a) 2σ Fig.8. TheformofGemp(x)fordifferentvaluesofσ2/N:10−3 dotted
(cid:32)√ (cid:33)(cid:35) line,10−4 dashedlineand10−5 solidline.
N(x−b)
−erf √ (7)
2σ √ √
Inparticular,δ = 2σerf−1(ε)/ N.Inthisintervaltheerror
In equation (7) a,b are the extrema of the range of vari- function differs from ±1 more than 1−ε. The value of N,
ability of the observable mean, or a = x (Σ ) and
0 S,min i.e.thenumberofgeneratedevents,fixesthestatisticalindeter-
b=x (Σ ).
0 S,max mination on the evaluation of the true epistemic uncertainty
ThebehaviorofG (x)isilllustratedinFig.8fordifferent
emp [a,b] of the simulated observable. One can choose N such
values of σ2/N: one can observe that the larger the value of
that the intervals in x over which both error functions in (7)
simulated events N, the more closely the curve resembles a
differ from ±1 by a predetermined amount are as small as
flatdistribution,whichreflectstheflatdistributionoftheinput
requested:toreduceδonemustincreaseN.Inpracticalterms,
parameter.Fig.8showsthatthebehavioroftheempiricalPDF
thismeansthattheintervalofvariabilityoftheobservablecan
is directly ruled by the sole parameter σ2/N.
be determined with any predefined statistical accuracy by two
ThebehavioroftheanalyticallycalculatedG (x)appears
emp Monte Carlo runs, corresponding to the values Σ =Σ
S S,min
consistent with the outcome observed in the toy Monte Carlo.
and Σ = Σ of the input parameter, encompassing an
S S,max
The key point to note is that G (x) from equation (7)
emp adequate number of events N.
is the distribution from which observable means are actually
From the previous discussion one can evince that for
sampled: by running many Monte Carlo simulations one can
G (x) to represent an adequate approximation of a flat
emp
obtain a statistical sample from it, as is illustrated in Figs. 6
distribution the condition
and 7. The histograms in these figures, which derive from √
the execution of a finite number of Monte Carlo simulations, σ<O> =σ/ N (cid:28)b−a (8)
are approximations of analytical curves as shown in Fig. 8.
shouldbesatisfied.Thisconditionrelatesthescaleofσ =
Theformtheyassumederivesfromtheinterplaybetweentwo √ <O>
σ/ N (that is the accuracy of a single Monte Carlo run) and
different statistical errors: the one coming from each Monte
the range of variability of the observable mean.
Carlo run (σ2/N) and the one coming from the finiteness
of the number of Monte Carlo runs, as noted earlier in
SectionIV-B. F. Approximationoftheparametricdependenceofanobserv-
able from an input unknown
E. Epistemic uncertainty of an observable In the previous discussion we have assumed a linear rela-
In some experimental scenarios knowledge of the range of tionship x0(ΣS) between the peak position of the observable
variability of the simulation outcome, rather than complete means distribution and the input parameter. This assumption
quantificationoftheuncertaintyoftheobservablemean,would can be verified by running a few Monte Carlo simulations
be sufficient. This range of variability can be considered as for some values of the input parameters within its range of
the epistemic uncertainty of the observable produced by the variability.
simulation; for such a concept to be meaningful in a context The results of performing this procedure in the context of
of uncertainty quantification, one should be able to estimate a our toy Monte Carlo are illustrated in Fig. 9. Each plot shows
priorithestatisticalindeterminationbywhichitisaffected.In the outcome of the simulation for 21 equidistant fixed values
the application scenario considered here this requirement cor- ofΣS intherangeofvariability1.0≤ΣS ≤1.2togetherwith
respondstoknowingthestatisticalindeterminationofextrema their 3σ confidence intervals. We display 3σ intervals instead
of the interval [a,b]. of conventional σ intervals for better visibility in the plots.
For this purpose we observe that for any fixed value of z Thefourplotscorrespondtodifferentpositionsofthedetector
and 0<ε<1, a value δ >0 exists, such that withrespecttotheprimaryparticlesource,withcenterlocated
(cid:12) √ √ (cid:12) at 0.1, 0.5, 0.9 and 1.3 arbitrary units from the source. The
(cid:12)(cid:12)erf( N(x−z)/( 2σ)(cid:12)(cid:12)<ε ⇐⇒ |x−z|<δ. number of events generated in each simulation is N = 108.
10
(cid:45) (cid:45)
(cid:72)(cid:76)100arbitraryunits 333...777368 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)(cid:76)1000arbitraryunits 222...000357 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
(cid:180) (cid:45) (cid:180) (cid:45)
O(cid:60)(cid:62)tl3.71 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)a(cid:76) O(cid:60)(cid:62)tl2.01 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)b(cid:76)
(cid:45) (cid:45)
1.04 1.08 1.12 1.16 1.04 1.08 1.12 1.16
(cid:83)S(cid:72)arbitraryunits(cid:76) (cid:83)S(cid:72)arbitraryunits(cid:76)
(cid:45)
(cid:76)units 6.16 (cid:45) (cid:45) (cid:45) (cid:76)units 2.53 (cid:45) (cid:45) (cid:45)
ary (cid:45) (cid:45) (cid:45) ary (cid:45) (cid:45) (cid:45) (cid:45)
arbitr 6.11 (cid:45) (cid:45) (cid:45) arbitr 2.51 (cid:45) (cid:45) (cid:45)
(cid:72)10000 6.05 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)10000 2.49 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
O(cid:60)(cid:62)(cid:180)tl5.99 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)c(cid:76) O(cid:60)(cid:62)(cid:180)tl2.48 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72)d(cid:76)
(cid:45)
1.04 1.08 1.12 1.16 1.04 1.08 1.12 1.16
(cid:83)S(cid:72)arbitraryunits(cid:76) (cid:83)S(cid:72)arbitraryunits(cid:76)
Fig.9. ThetestontheformofthefunctionΣS(x):weplotfor21equidistantandfixedvaluesofΣS theoutcomesofaMCrun,eachencompassing108
events,togetherwiththecorrespondingstatisticalerrors.Thefourplotscorrespondtodifferentdistanceofthedetectorfromtheprimaryparticlesource:the
centerofthedetectorisplacedat(a):0.1,(b):0.5,(c):0.9and(d)1.3(inarbitraryunits).Intheplotweshowalsothebestlinearfittothedata(seetext).
Theresultofalinearfittothesimulateddataissuperimposed tainty quantification is reduced to determining the parameters
to the plots. A linear relation x (Σ ) qualitatively appears a required for its calculation: this problem in turn consists of
0 S
justified approximation, although the statistical degradation of determining the unknown function x (Σ), which represents
0
thequalityofalinearfitisclearlyvisible,whentheobservable the response of the observable means to the variation of the
isscoredatlargerdistancefromthesource.Itisquantitatively input parameter, or its inverse Σ(x ).
0
supportedbythep-valuesofthelinearfit:0.995,0.993,0.997 From a Monte Carlo user perspective, the process involves
and0.985respectively,forincreasingdistancefromthesource. the investigation of two properties of the output.
Atincreasingdistancesofthedetectorfromthesourcealarger
The first concerns the extrema of the output variability
numberofeventscanbegenerated,ifnecessary,toreducethe
intervalfortheobservablemeans,thatisthemaximumoutput
sizeoftheconfidenceintervalsfortheparametersofthefit:in
uncertainty. These extrema can be found with an arbitrary
practical cases a detailed statistical analysis is recommended,
predefinedprecisionwithtwoMonteCarloruns,thefirstusing
since the knowledge of the expected PDF G(x) is affected as input Σ=Σ , the second using as input Σ=Σ . We
min max
by the uncertainties in the best-fit parameters for the relation
stress that this is not a full quantification of the uncertainty, if
x0(ΣS). nothing is known about the relation x (Σ); these two values,
0
This method can be used more generally to investigate
together with the statistical uncertainties on their determina-
possibleapproximationsofthex (Σ )functionalrelationship.
0 S tion, define an interval of variability for the output, not its
The resulting fitted function x (Σ) can be inserted directly in
0 probability distribution inside this interval.
equation (2) to obtain the output PDF.
A remark is here necessary, regarding the possibility of an
input PDF not supported on a bounded interval - this is the
VII. METHODFORESTIMATINGOUTPUTUNCERTAINTIES case for instance of a normal distribution which is supported
The theoretical findings of section IV and the heuristic on (−∞,+∞): we suggest in this case to first identify an
investigationofsectionVIinasimpleapplicationenvironment interval of interest for the input, for instance associated with
support the definition of a path for the calculation of the a given confidence interval for the input, and to proceed on
uncertainties of simulation observables in a generic scenario. with the analysis. In this case one will not obtain the output
We summarize it here from the practical perspective of a PDF, but the output PDF conditioned to the probability of the
Monte Carlo simulation user. This procedure assumes that the input to be in this interval.
uncertainties of input parameters are known. Approximate knowledge of the function x (Σ) can be
0
Equation (2) states that the form of the output PDF is achieved through statistical investigation, by performing a
assured by theory. As a consequence, the process of uncer- numberofMonteCarlosimulations,involvingdifferentvalues