1

42. Monte Carlo Techniques

42. Monte Carlo Techniques
Revised September 2021 by G. Cowan (RHUL).
Monte Carlo techniques are often the only practical way to evaluate difficult integrals or to
sample random variables governed by complicated probability density functions. Here we describe
an assortment of methods for sampling some commonly occurring probability density functions.

42.1 Sampling the uniform distribution
Most Monte Carlo sampling or integration techniques assume a “random number generator,”
which generates uniform statistically independent values on the half open interval [0, 1); for reviews
see, e.g., Refs. [1, 2].
Uniform random number generators are available in software libraries such as CLHEP [3], andROOT [4]. For example, in addition to a basic congruential generator TRandom (see below), ROOT
provides three more sophisticated routines: TRandom1 implements the RANLUX generator [5] basedon the method by Lüscher, and allows the user to select different quality levels, trading off quality
with speed; TRandom2 is based on the maximally equidistributed combined Tausworthe generator
by L’Ecuyer [6]; the TRandom3 generator implements the Mersenne twister algorithm of Matsumoto
and Nishimura [7]. All of the algorithms produce a periodic sequence of numbers, and to obtaineffectively random values, one must not use more than a small subset of a single period. The
Mersenne twister algorithm has an extremely long period of 219937 − 1. A recent review of highquality random number generators can be found in Ref. [8].
The performance of the generators can be investigated with tests such as DIEHARD [9] or
TestU01 [10]. Many commonly available congruential generators fail these tests and often havesequences (typically with periods less than 232 ), which can be easily exhausted on modern computers. A short period is a problem for the TRandom generator in ROOT, which, however, has theadvantage that its state is stored in a single 32-bit word. The generators TRandom1, TRandom2, orTRandom3 have much longer periods, with TRandom3 being recommended by the ROOT authors asproviding the best combination of speed and good random properties. For further information see,e.g., Ref. [11].

42.2 Inverse transform method
If the desired probability density function is f (x) on the range −∞ < x < ∞, its cumulative
distribution function (expressing the probability that x ≤ a) is given by Eq. (39.6). If a is chosenwith probability density f (a), then the integrated probability up to point a, F (a), is itself arandom variable which will occur with uniform probability density on [0, 1]. Suppose u is generated
according to a uniformly distributed in (0, 1). If x can take on any value, and ignoring the endpoints,
we can then find a unique x chosen from the p.d.f. f (x) for a given u if we set
u = F (x) ,(42.1)

provided we can find an inverse of F , defined by
x = F −1 (u) .

(42.2)

This method is shown in Fig. 42.1a. It is most convenient when one can calculate by hand theinverse function of the indefinite integral of f . This is the case for some common functions f (x)such as exp(x), (1 − x)n , and 1/(1 + x2 ) (Cauchy or Breit-Wigner), although it does not necessarily produce the fastest generator. Standard libraries contain software to implement this method
numerically, working from functions or histograms in one or more dimensions, e.g., the UNU.RAN
S. Navas et al. (Particle Data Group), Phys. Rev. D 110, 030001 (2024)
31st May, 2024 10:17am

42. Monte Carlo Techniques

2

1

(a)

Continuous
distribution

F(x)
u
0

x
x = F −1(u)

1

(b)

F(x)
u

Discrete
distribution} f(xk)0x

xk xk+1

Figure 42.1: Use of a random number u chosen from a uniform distribution (0,1) to find a random
number x from a distribution with cumulative distribution function F (x).
package [12], available in ROOT. For a discrete distribution, F (x) will have a discontinuous jumpof size f (xk ) at each allowed xk , k = 1, 2, · · · . Choose u from a uniform distribution on (0,1) as
before. Find xk such that
F (xk−1 ) < u ≤ F (xk ) ≡ Prob (x ≤ xk ) =kXf (xi ) ;(42.3)

i=1

then xk is the value we seek (note: F (x0 ) ≡ 0). This algorithm is illustrated in Fig. 42.1b.

42.3 Acceptance-rejection method (Von Neumann)
Very commonly an analytic form for F (x) is unknown or too complex to work with, so thatobtaining an inverse as in Eq. (42.2) is impractical. We suppose that for any given value of x, theprobability density function f (x) can be computed, and further that enough is known about f (x)that we can enclose it entirely inside a shape which is C times an easily generated distribution
h(x), as illustrated in Fig. 42.2. That is, Ch(x) ≥ f (x) must hold for all x. Frequently h(x) is
uniform or is a normalized sum of uniform distributions. Note that both f (x) and h(x) must be
normalized to unit area, and therefore, the proportionality constant C > 1. To generate f (x), first
generate a candidate x according to h(x). Calculate f (x) and the height of the envelope C h(x);generate u and test if uC h(x) ≤ f (x). If so, accept x; if not reject x and try again. If we regard
x and uC h(x) as the abscissa and ordinate of a point in a two-dimensional plot, these points will
populate the entire area C h(x) in a smooth manner; then we accept those which fall under f (x).
31st May, 2024

42. Monte Carlo Techniques3

C h(x)

(a)

f(x)
(b)
C h(x)f(x)x
Figure 42.2: Illustration of the acceptance-rejection method. Random points are chosen insidethe upper bounding figure, and rejected if the ordinate exceeds f (x). The lower figure illustrates amethod to increase the efficiency (see text).

The efficiency is the ratio of areas, which must equal 1/C; therefore we must keep C as close as
possible to 1.0. Therefore, we try to choose C h(x) to be as close to f (x) as convenience dictates,
as in the lower part of Fig. 42.2.

42.4 Algorithms
Algorithms for generating random numbers belonging to many different distributions are given
for example by Press [13], Ahrens and Dieter [14], Rubinstein [15], Devroye [16], Walck [17] andGentle [18]. For many distributions, alternative algorithms exist, varying in complexity, speed,
and accuracy. For time-critical applications, these algorithms may be coded in-line to remove thesignificant overhead often encountered in making function calls.
In the examples given below, we use the notation for the variables and parameters given in Table 39.1. Variables named “u” are assumed to be independent and uniform on [0,1). Denominators
must be verified to be non-zero where relevant.
42.4.1 Exponential decay
This is a common application of the inverse transform method, and uses the fact that if u
is uniformly distributed in [0, 1], then (1 − u) is as well. Consider an exponential p.d.f. f (t) =(1/τ ) exp(−t/τ ) that is truncated so as to lie between two values, a and b, and renormalized to unit
area. To generate decay times t according to this p.d.f., first let α = exp(−a/τ ) and β = exp(−b/τ );31st May, 2024

42. Monte Carlo Techniques

4then generate u and lett = −τ ln(β + u(α − β)).

(42.4)For (a, b) = (0, ∞), we have simply t = −τ ln u. (See also Sec. 42.4.6.)
42.4.2 Isotropic direction in 3D
Isotropy means the density is proportional to solid angle, the differential element of which isdΩ = d(cos θ)dφ. Hence cos θ is uniform (2u1 − 1) and φ is uniform (2πu2 ). For alternative
generation of sin φ and cos φ, see the next subsection.
42.4.3 Sine and cosine of random angle in 2D
Generate u1 and u2 . Then v1 = 2u1 − 1 is uniform on (−1,1), and v2 = u2 is uniform on (0,1).
Calculate r2 = v12 + v22 . If r2 > 1, start over. Otherwise, the sine (S) and cosine (C) of a random
angle (i.e., uniformly distributed between zero and 2π) are given by
S = 2v1 v2 /r2and

C = (v12 − v22 )/r2 .

(42.5)

p(42.6)

42.4.4 Gaussian distribution
If u1 and u2 are uniform on (0,1), thenpz1 = sin(2πu1 ) −2 ln u2andz2 = cos(2πu1 ) −2 ln u2

are independent and Gaussian distributed with mean 0 and σ = 1.
There are many variants of this basic algorithm, which may be faster. For example, construct
v1 = 2u1 − 1 and v2 = 2u2 − 1, which are uniform on (−1,1). Calculate r2 = v12 + v22 , and if r2 > 1
start over. If r2 < 1, it is uniform on (0,1). Thensz1 = v 1−2 ln r2
r2

s

and

z2 = v 2−2 ln r2
r2

(42.7)are independent numbers chosen from a normal distribution with mean 0 and variance 1. zi0 = µ+σzi
distributes with mean µ and variance σ 2 .
For a multivariate Gaussian with an n × n covariance matrix V , one can start by generatingn independent Gaussian variables, {ηj }, with mean 0 and variance 1 as above. Then the new set
P{xi } is obtained as xi = µi + j Lij ηj , where µi is the mean of xi , and Lij are the components ofL, the unique lower triangular matrix that fulfils V = LLT . The matrix L can be easily computedby the following recursive relation (Cholesky’s method):

Ljj = Vjj −

j−1
X

1/2

L2jk 

,(42.8a)

k=1

Lij =

Vij −

Pj−1

k=1 Lik Ljk

Ljj, j = 1, ..., n ; i = j + 1, ..., n,(42.8b)where Vij = ρij σi σj are the components of V . For n = 2 one hasL=

σ1 p 0ρσ21 − ρ2 σ2!

,(42.9)

andp therefore the correlated Gaussian variables are generated as x1 = µ1 + σ1 η1 , x2 = µ2 + ρσ2 η1 +1 − ρ 2 σ2 η2 .
31st May, 2024

42. Monte Carlo Techniques

542.4.5 χ2 (n) distribution
To generate a variable following the χ2 distribution for n degrees of freedom, use the Gamma
distribution with k = n/2 and λ = 1/2 using the method of Sec. 42.4.6.
42.4.6 Gamma distribution
All of the following algorithms are given for λ = 1. For λ 6= 1, divide the resulting random
number x by λ.
• If k = 1 (the exponential distribution), accept x = − ln u. (See also Sec. 42.4.1.)• If 0 < k < 1, initialize with v1 = (e + k)/e (with e = 2.71828... being the natural log base).
Generate u1 , u2 . Define v2 = v1 u1 .
1/kCase 1: v2 ≤ 1. Define x = v2 . If u2 ≤ e−x , accept x and stop, else restart by generating
new u1 , u2 .
Case 2: v2 > 1. Define x = −ln([v1 − v2 ]/k). If u2 ≤ xk−1 , accept x and stop, else restartby generating new u1 , u2 . Note that, for k < 1, the probability density has a pole at
x = 0, so that return values of zero due to underflow must be accepted or otherwise
dealt with.
• Otherwise, if k p> 1, initialize with c = 3k−0.75. Generate u1 and compute v1 = u1 (1−u1 ) and
v2 = (u1 − 0.5) c/v1 . If x = k + v2 − 1 ≤ 0, go back and generate new u1 ; otherwise generate
u2 and compute v3 = 64v13 u22 . If v3 ≤ 1 − 2v22 /x or if ln v3 ≤ 2{[k − 1] ln[x/(k − 1)] − v2 },accept x and stop; otherwise go back and generate new u1 .
42.4.7 Binomial distribution
Begin with k = 0 and generate u uniform in [0, 1). Compute Pk = (1 − p)n and store Pk into
B. If u ≤ B accept rk = k and stop. Otherwise, increment k by one; compute the next Pk as
Pk · (p/(1 − p)) · (n − k)/(k + 1); add this to B. Again, if u ≤ B, accept rk = k and stop, otherwise
iterate until a value is accepted. If p > 1/2, it will be more efficient to generate r from f (r; n, q),i.e., with p and q interchanged, and then set rk = n − r.
42.4.8 Poisson distribution
Iterate until a successful choice is made: Begin with k = 1 and set A = 1 to start. Generate u.
Replace A with uA; if now A < exp(−µ), where µ is the Poisson parameter, accept nk = k − 1 andstop. Otherwise increment k by 1, generate a new u and repeat, always starting with the value ofA left from the previous try.
Note that the Poisson generator used in ROOT’s TRandom classes before version 5.12 (including the derived classes TRandom1, TRandom2, TRandom3) uses a Gaussian approximation whenµ exceeds a given threshold. This may be satisfactory (and much faster) for some applications.
To do this, generate z from a Gaussian with zero mean and unit standard deviation; then use√
x = max(0, [µ + z µ + 0.5]) where [ ] signifies the greatest integer ≤ the expression. The routinesfrom Numerical Recipes [13] and CLHEP’s routine RandPoisson do not make this approximation(see, e.g., Ref. [11]).
42.4.9 Student’s t distribution−2/nGenerate u1 and u2 uniform in (0, 1); then t = sin(2πu1 )[n(u2− 1)]1/2 follows the Student’s
t distribution for n > 0 degrees of freedom (n not necessarily an integer).
31st May, 2024

42. Monte Carlo Techniques

6

Alternatively, generate x from a Gaussian with mean 0 and σ 2 = 1 according to the method of
42.4.4. Next generate y, anpindependent gamma random variate, according to 42.4.6 with λ = 1/2
and k = n/2. Then z = x/ y/n is distributed as a t with n degrees of freedom.
For the special case n = 1, the Breit-Wigner distribution, generate u1 and u2 ; set v1 = 2u1 − 1and v2 = 2u2 − 1. If v12 + v22 ≤ 1 accept z = v1 /v2 as a Breit-Wigner distribution with unit
area, center at 0.0, and FWHM 2.0. Otherwise start over. For center M0 and FWHM Γ , use
W = zΓ/2 + M0 .
42.4.10 Beta distribution
The choice of an appropriate algorithm for generation of beta distributed random numbers
depends on the values of the parameters α and β. For, e.g., α = 1, one can use the transformation
method to find x = 1 − u1/β , and similarly if β = 1 one has x = u1/α . For more general cases see,e.g., Refs. [17, 18] and references therein.

42.5 Importance sampling and weighted Monte Carlo
Often the goal of a Monte Carlo calculation is to determine an expectation value of a function
h(x), where x is a (single or vector) random variable that follows a pdf f (x),Z

Ef [h(x)] =h(x)f (x) dx ≡ µ .

(42.10)A Monte-Carlo estimator µ̂MC for µ is the average of N values of h(x) where x is sampled (generated)Pfrom f (x), i.e., µ̂MC = N1 N
i=1 h(xi ). This has a varianceV [µ̂MC ] =11 
Vf [h(x)] =Ef [h2 (x)] − µ2 .
NN(42.11)

By using the method of importance sampling, one can achieve a reduction in this variance andthus a more accurate determination of the expectation value for a given number of random values
generated. The key idea is to rewrite the expectation value in Eq. (42.10) asZµ=h(x)f (x)h(x)f (x)g(x) dx = Eg
g(x)
g(x)


Z

h(x)f (x) dx =,(42.12)where g(x) is any other pdf of x with the same support as f (x) (i.e., nonzero for the same regionof x). Thus the desired quantity µ is the expectation value with respect to g of h(x)f (x)/g(x). It
can be estimated by generating N values of x sampled from g and computingµ̂IS =N
1 X
h(xi )f (xi ).
N i=1 g(xi )(42.13)The variance of µ̂IS is given by
1
h(x)f (x)1=V [µ̂IS ] = Vg
N
g(x)
N


"



Eg

#

h2 (x)f 2 (x)− µ2
g 2 (x)!

,(42.14)

By choosing g(x) such that h(x)f (x)/g(x) is as constant as possible, the variance of µ̂IS can besubstantially reduced. One can show (see, e.g., Refs. [19, 20]) that the variance is minimized wheng(x) ∝ |h(x)|f (x).
31st May, 2024

742. Monte Carlo Techniques

An alternative importance sampling estimator can be constructed by replacing the number ofP
generated values N in Eq. (42.13) by the sum N
i=1 f (xi )/g(xi ). This can given an even smaller
variance at the price of a small bias. It can be of further advantage in problems where the pdf f (x)is known only up to a normalization constant, which then cancels (see Refs. [19, 20]).
A closely related application of importance sampling is the use of weighted Monte Carlo tocompute the probability Pf (x ∈ A) for x to be in a specified region A:
Pf (x ∈ A) =Z

f (x) dx .

(42.15)A

It may be, however, that one does not have an MC model capable of generating x ∼ f (x), butrather one can generate x according to a different density g(x). The probability Pf (x ∈ A) can be
writtenR f (x)Z
A g(x) g(x) dx
g(x) dx = Eg [w(x)|x ∈ A] Pg (x ∈ A) ,Pf (x ∈ A) = R
A g(x) dx(42.16)A

R

where w(x) = f (x)/g(x) is the weight function and Pg (x ∈ A) = A g(x) dx is the probability tofind x ∈ A assuming x ∼ g(x). That is, Pf (x ∈ A) is the conditional expectation value of w(x)with respect to g(x) given x ∈ A multiplied by the probability to find x ∈ A under assumption ofg(x).
Suppose N values of x are generated according to g(x) and m of them are found in the region
A. Then the probability to be in A for x ∼ g(x) can be estimated by m/N , and the expectation
value above can be obtained from the average of the weights in A. Therefore the desired probability
Pf (x ∈ A) can be estimated usingP̂f (x ∈ A) =

m
m
1 X
m
1 X
wi ×
wi ,
=
m i=1
N
N i=1

(42.17)

where wi = w(xi ) and the sum includes only the m values of x found in A. That is, when generating
the x values according to g(x) instead of f (x), the number of events m found in A is replaced by
the corresponding sum of weights. The variance of P̂f (x ∈ A) can be found fromVb [P̂f (x ∈ A)] =

m
1 X
w2 .
N 2 i=1 i(42.18)

By choosing g(x) so that a larger fraction of x values are sampled in the “important” region A, one
can reduce the variance of the estimated probability for a given total number of generated values.

42.6 Markov Chain Monte Carlo
In applications involving generation of random numbers following a multivariate distribution
with a high number of dimensions, the transformation method may not be possible and theacceptance-rejection technique may have too low of an efficiency to be practical. If it is not required
to have independent random values, but only that they follow a certain distribution, then Markov
Chain Monte Carlo (MCMC) methods can be used. In depth treatments of MCMC can be found,e.g., in the texts by Robert and Casella [19], Liu [20], and the review by Neal [21]. HEP-oriented
software for MCMC is available from the Bayesian Analysis Toolkit (BAT) [22, 23].
31st May, 2024

842. Monte Carlo Techniques

MCMC is particularly useful in connection with Bayesian statistics, where a p.d.f. p(θ) for
an n-dimensional vector of parameters θ = (θ1 , . . . , θn ) is obtained, and one needs the marginal
distribution of a subset of the components. Here one samples θ from p(θ) and simply records themarginal distribution for the components of interest.
A simple and broadly applicable MCMC method is the Metropolis-Hastings algorithm, whichallows one to generate multidimensional points θ distributed according to a target p.d.f. that is
proportional to a given function p(θ). It is not necessary to have p(θ) normalized to unit area,
which is useful in Bayesian statistics, as posterior probability densities are often determined only
up to an unknown normalization constant.
To generate points that follow p(θ), one first needs a proposal p.d.f. q(θ; θ0 ), which can be(almost) any p.d.f. from which independent random values θ can be generated, and which contains
as a parameter another point in the same space θ0 . For example, a multivariate Gaussian centeredabout θ0 can be used. Beginning at an arbitrary starting point θ0 , the Hastings algorithm iterates
the following steps:
1. Generate a value θ using the proposal density q(θ; θ0 );hip(θ)q(θ0 ;θ)2. Form the Hastings test ratio, α = min 1, p(θ
;0 )q(θ;θ0 )3. Generate a value u uniformly distributed in [0, 1];
4. If u ≤ α, take θ1 = θ. Otherwise, repeat the old point, i.e., θ1 = θ0 .
5. Set θ0 = θ1 and return to step 1.
If one takes the proposal density to be symmetric in θ and θ0 , then this is the Metropolis -Hastings
algorithm, and the test ratio becomes α = min[1, p(θ)/p(θ0 )]. That is, if the proposed θ is at avalue of probability higher than θ0 , the step is taken. If the proposed step is rejected, the old pointis repeated.
Methods for assessing and optimizing the performance of the algorithm are discussed in, e.g.,
Refs. [19–21]. One can, for example, examine the autocorrelation as a function of the lag k, i.e.,the correlation of a sampled point with that k steps removed. This should decrease as quickly as
possible for increasing k.
Generally one chooses the proposal density so as to optimize some quality measure such as the
autocorrelation. For certain problems it has been shown that one achieves optimal performance
when the acceptance fraction, that is, the fraction of points with u ≤ α, is around 40%. This can be
adjusted by varying the width of the proposal density. For example, one can use for the proposal
p.d.f. a multivariate Gaussian with the same covariance matrix as that of the target p.d.f., butscaled by a constant.
References[1] F. James, Comput. Phys. Commun. 60, 329 (1990).
[2] P. L’Ecuyer, Proc. 1997 Winter Simulation Conference, IEEE Press, Dec. 1997, 127–134.
[3] L. Lonnblad, Comput. Phys. Commun. 84, 307 (1994).
[4] R. Brun and F. Rademakers, Nucl. Instrum. Meth. A389, 81 (1997); See also root.cern.ch.
[5] F. James, Comput. Phys. Commun. 79, 111 (1994), [Erratum: Comput. Phys. Commun.97,357(1996)]; M. Luscher, Comput. Phys. Commun. 79, 100 (1994), [hep-lat/9309020].
[6] P. L’Ecuyer, Mathematics of Computation, 65, 213 (1996) and 65, 225 (1999).
[7] M. Matsumoto and T. Nishimura, ACM Transactions on Modeling and Computer Simulation,
Vol. 8, No. 1, January 1998, 3–30.
31st May, 2024

942. Monte Carlo Techniques[8] F. James and L. Moneta, Comput. Softw. Big Sci. 4, 2. 12 p (2019), 4 figures,[arXiv:1903.01247], URL http://cds.cern.ch/record/2678858.
[9] Much of DIEHARD is described in: G. Marsaglia, A Current View of Random Number Generators, keynote address, Computer Science and Statistics: 16th Symposium on the Interface,
Elsevier (1985).

[10] P. L’Ecuyer and R. Simard, ACM Transactions on Mathematical Software 33, 4, Article 1,
December 2007.
[11] J. Heinrich, CDF Note CDF/MEMO/STATISTICS/PUBLIC /8032, 2006.
[12] UNU.RAN is described at statmath.wu.ac.at/software/ unuran; See also W. Hörmann, J.
Leydold, and G. Derflinger, Automatic Nonuniform Random Variate Generation, (Springer,
New York, 2004).
[13] W.H. Press et al., Numerical Recipes, 3rd edition, (Cambridge University Press, New York,
2007).
[14] J.H. Ahrens and U. Dieter, Computing 12, 223 (1974).
[15] R.Y. Rubinstein, Simulation and the Monte Carlo Method, (John Wiley and Sons, Inc., New
York, 1981).
[16] L. Devroye, Non-Uniform Random Variate Generation, (Springer-Verlag, New York, 1986);
Available online at luc.devroye.org/rnbookindex.html.
[17] C. Walck, Handbook on Statistical Distributions for Experimentalists, University of Stockholm
Report SUF-PFY/96-01, available from www.fysik.su.se/˜walck.
[18] J.E. Gentle, Random Number Generation and Monte Carlo Methods, 2nd ed., (Springer, New
York, 2003).
[19] C.P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., (Springer, New York,
2004).
[20] J.S. Liu, Monte Carlo Strategies in Scientific Computing, (Springer, New York, 2001).
[21] R.M. Neal, Probabilistic Inference Using Markov Chain Monte Carlo Methods, Technical
Report CRG-TR-93-1, Dept. of Computer Science, University of Toronto, available fromwww.cs.toronto.edu/˜radford/res-mcmc.html .
[22] A. Caldwell, D. Kollar, K. Kröninger, Comput. Phys. Commun. 180 (2009) pages 2197-2209.
[23] Schulz, O., Beaujean, F., Caldwell, A. et al., SN Computer Science 2, 210 (2021).

31st May, 2024