Lecture Notes on Monte Carlo Methods

Fall Semester, 2005

Courant Institute of Mathematical Sciences, NYU
Jonathan Goodman, goodman@[Link]

Introduction
created October 15, 2005

This is a set of lecture notes for a graduate class on Monte Carlo methods
given at the Courant Institute of Mathematical Sciences at NYU in the fall of
2005. I am motivated by seeing many new people begin to use Monte Carlo.
Many new techniques and applications are emerging. I hope the class will cover
the basic principles quickly but with enough depth to be helpful to practition-
ers. Then we will cover a few of the many advanced areas, hopefully in a way
that people from diverse fields can understand. Given the variety of people us-
ing Monte Carlo and the range of application areas, it is surprising how much
common ground can be found.
Roughly speaking, Monte Carlo1 means computing using random numbers.
It is helpful to refine this definition and distinguish between true Monte Carlo
and simulation. We take simulation to mean generating individual random
objects faithfully according to some model. For example, we might want to see
what shape clouds come from a specific model of cloud formation that involves
randomness. The point of simulation might not be gather detailed statistics,
but just to see what a few random objects look like.
By contrast, Monte Carlo uses random numbers as a means to evaluate quan-
tities that themselves are not random. For example, suppose f (x) is the proba-
bility density for a one dimensional random variable, X. One way to evaluate2
A = E[X] is to generate many samples (random variables Xk with probability
density f ) and average them. More generally, we could generate a number of
random objects (e.g. clouds) and collect interesting statistics about them. The
difference is that the expected value of X or of some more complex statistic is a
property of a random variable but is not itself random. Therefore, it may be pos-
sible to evaluate A without generating random R samples with probability density
f (x). For example, we could estimate A = xf (x)dx by numerical quadrature.
Practitioners often find clever methods that are better (faster or more accurate)
than plain simulation. Simulation is mostly programming but Monte Carlo is
all about devising and understanding new computational strategies. Each piece
of creativity will improve your results.
When choosing a Monte Carlo method for a given problem, one should have
a strong bias against Monte Carlo at all. A practical deterministic method
almost always is better than Monte Carlo. If you have an integral in fewer
than, say, four variables, do it by deterministic quadrature. If you have a
1 The name comes from Monte Carlo, a traditional gambling center in Europe, another

place where random numbers are important.

2 See the next section for basic definitions and notation.

1
Markov chain with less than a thousand states or a diffusion process in less
than four dimensions, √ use a backward equation. The statistical error in Monte
Carlo is at least O(1/ L) (which isn’t very small), where L is the number of
samples3 It is a rare deterministic method that is worse than this.
Keep this in mind when reading Monte Carlo books. We authors often
use one dimensional examples to illustrate specific principles although we know
that this is not the fastest or most accurate way to compute a one dimensional
integral. The real world is full of problems where Monte Carlo is the method of
choice.

1 Overview
An example will illustrate the themes discussed in these notes. One part of a
physics problem reduces to calculating the integral

e−λ|x−y|
Z Z
A(λ) = dxdy , (1)
|x|<1 |y|<1 |x − y|

Here, the x and y variables both are in R3 . This is a six dimensional integral
(maybe reducible to five or even four dimensions), which makes direct quadra-
ture difficult. The problem is first to calculate A(λ) for various λ values and
next to find λ∗ that gives A(λ∗ ) = 21 .
The first step in evaluating A(λ) is to express in terms of the expected value
of a random variable. Here we will have

A(λ) = Const · B(λ) , B(λ) = Eλ [Z] . (2)

The Eλ [·] on the right means that the probability distribution of Z depends on
λ. For example, we could let X and  Y be independent random variables with
uniform probability density in B = x ∈ R3 with |x| < 1 . This probability
density is (recall that the volume of B is 34 π)
 3
4π if x ∈ B
f (x) =
0 otherwise.

Then multiplying and dividing by the normalizing 34 π gives

2 Z Z
e−λ|x−y|

4
A(λ) = f (x)f (y)dxdy
3π |x − y|
 2  −λ|X−Y | 
4 e
= E
3π |X − Y |
 2
4
= E[Z] ,
3π
3 We use L for the number of samples, the run length so that n always can be the number

of components of a multivariate random variable.

2
where
e−λ|X−Y |
Z= . (3)
|X − Y |
We do not have a formula for the probability density of Z, but we can generate
a random sample by choosing X and Y independently in B then applying (3).
To evaluate B(λ), we generate L independent Z samples, Z1 , . . ., ZL , and
use the estimator4
L
1X
E[Z] = B(λ) ≈ B(λ)
b = Zk . (4)
L
k=1

The law of large numbers states that B b → B as L → ∞. It may take quite a

large L, and lots of computer time, to get B b close enough to B to satisfy us.
Computational error bars tell us how far B is likely to be from B. We will
b
follow the statistician’s habit of expressing this with a confidence interval. For
example, we might say that the probability that B is not between B b − r and
B + r is about 5%. A Monte Carlo practitioner who neglects error bars deserves
b
to get the wrong answer, and will. It is not always necessary to present the
error bars to the consumers of your results, just as you don’t show them how
you debugged and tested your computer code.
Applications demand more than just B(λ b j ) for a few λj . We want to know
the function B(λ), the shape of its graph, its derivatives, the inverse function,
etc. The statistical noise in the estimators B(λ
b j ) may be amplified in finding
this other properties of B(λ). There often are better ways5 . Sensitivity analysis
provides estimates for derivatives of B(λ) having less noise than the naive finite
difference (see below). Stochastic approximation provides more sophisticated

2
ways to find solutions of equations such as B(λ) = 12 pi 4 (solve for λ so that
A(λ) = 12 ).
There are other estimators of B(λ) besides (3),(4). Variance reduction means
searching for alternatives with (hopefully) less statistical error. For example, if
we could sample (X, Y ) pairs from the 6 dimensional probability density

f (x)f (y)
g(x, y) = Const · , (5)
|x − y|

then B(λ) = Eg [e−λ|X−Y | ] would give more accurate estimates of B for the
same number of samples. This is an example of importance sampling. It is
harder to create samples from the density (5) than to sample X and Y from f
independently. Powerful sampling methods such as rejection and Markov chain
Monte Carlo (MCMC) will be very handy.
The simple estimator (3),(4) is unbiased, which means that B(λ) = E[B(λ)].
b
Many estimators have bias as well as statistical error. For example, since B(λ) is
4 Statisticians put a “hat” on a quantity to indicate a statistical estimate of it. In this way,

the statistical estimate of B is B.

b
5 Ask A.C. for a joke with this punch line.

3
 1 pi 2


a nonlinear function of λ, an estimator, λ
c∗ , of λ∗ with B(λ∗ ) =
2 4 probably
has bias  
bias = E[λ∗ ] − E[λc∗ ] = O 1 .
L
There may be a tradeoff between statistical
 error and bias. For example, (3),(4)
has statistical error B(λ) − B(λ) = O √1L . The finite difference estimator
b

c0 (λ) = B(λ + ∆λ) − B(λ)

b b
B (6)
∆λ
has bias

c0 − B 0 = B(λ + ∆λ) − B(λ) − B 0 (λ) = O(∆λ) .

h i
E B
∆λ
Assuming B(λ+∆λ)
b and B(λ)
b have independent statistical error, the statistical
√
0
error of B is on the order of 1/∆λ L. Altogether, we minimize the total error
c
by (setting all constants equal to one)
 
1 1
min (bias + noise) = min ∆λ + √ = 1/4 .
∆λ ∆λ ∆λ L L
Sensitivity analysis and variance reduction will improve this.

2 Background
The prerequisites for the class are a course on stochastic processes (Stochastic
Calculus at the Courant Institute) and some experience with scientific comput-
ing at the level of our Scientific Computing class. We will use linear algebra
and multivariate calculus at the level of the Courant Institute beginning classes.
The discussion will be as informal mathematically as possible.
Random variables generally are denoted by capitol letters, X, Y , T , etc.,
with specific values denoted by lower case: x, y, t, etc. A random element of Rn
would be called a random vector or a multivariate random variable and could
be written X = (X1 , . . . , Xn ). We have a scalar random variable when n = 1.
The probability density for X might be called f (x), so that
Z
P (A) = P (X ∈ A) = f (x)dx . (7)
A

Here P (A) is the probability of the event A ⊆ Rn . The expected value of X

with the law f is Z
Ef [X] = xf (x)dx .
Rn
We write E[·] for Ef [·] when the f is clear. In one dimension it may be clearer
to use old fashion differential notation
P (x ≤ X ≤ x + dx) = f (x)dx . (8)

4
Here the event A is the small interval (x, x + dx). The numbers P (A) form a
probability measure even if they are not given in terms of a density as in (1.1).
The law of a random variable either is its probability density or its probability
measure. We write X ∼ f or X ∼ P to indicate that X has the law given by
f or P . The vector notation above notwithstanding, we often write Xk ∼ f
to indicate a sequence random variables each with the law f . If Xk ∈ Rn ,
we can indicate vector components by Xk = (Xk1 , . . . , Xkn ). We say the Xk
are samples of the law f or samples of the random variable X. Monte Carlo
computations involve thousands or millions of samples.
A standard uniform random variable is a scalar, U , with probability density
f (u) = 1 for 0 ≤ u ≤ 1 and f (u) = 0 otherwise. A pseudo random number
generator is a piece of computer code that produces a sequence Uk , k = 1, 2, . . .,
that resembles a sequence of independent samples of a standard uniform. This
is not to be confused with quasi random numbers which are deliberately not
random, in an attempt to be more uniform than actual random numbers. The
term “pseudo random” indicates that the output of the pseudo random num-
ber generator is not actually random. If you run the pseudo random number
generator twice, with the same seed (see a later lecture for a more technical
discussion), you will get the same Uk .
The central limit theorem underlies much Monte Carlo error analysis. The
simplest case is Y1 , . . ., YL , independent samples of a scalar random variable,
Y with E[Y ] = 0 and E[Y 2 ] = σY2 . The theorem is that the law of ZL =
PL
L−1/2 k=1 Yk converges to a Gaussian law with mean zero and variance σY2 as
L → ∞. The next case is a multivariate random variable, Y with mean zero and
covariance CY = E[Y Y t ]. The law of ZL converges to a multivariate normal
with mean zero and covariance CY .
Wick’s theorem is a recipe that evaluates any moment of a multivariate
normal. Suppose l1 (Y ), . . ., l2m (Y ) are an even number of linear functionals.6
A pairing is a collection of m pairs {{j1 , k1 } , · · · , {jm , km }}, so that each of the
numbers 1, . . . , 2m appears exactly once. Reordering the pairs or exchanging
the numbers in a pair gives the same pairing. The pairings of {1, 2, 3, 4} are
{{1, 2} , {3, 4}}, {{1, 3} , {2, 4}}, and {{1, 4} , {2, 3}}. The number of pairings of
2m numbers is7 (2m − 1)(2m − 3) · · · 3. Wick’s theorem states that if Y is a
mean zero multivariate normal, then
X
E [l1 (Y ) · · · · · l2m (Y )] = E [lj1 (Y )lk1 (Y )] · · · · · E [ljm (Y )lkm (Y )] (9)
pairings

An example illustrates many features of this formula. Let (Y1 , Y2 ) be a mean

zero bivariate normal with var(Y1 ) = 2, var(Y2 ) = 5, and cov(Y1 , Y2 ) = 3, then
6 A linear functional is a scalar linear function of Y . If Y is a column vector, l(Y ) can be

represented as l · Y , for some row vector, also called l. Scalar functions, particularly of high
or infinite dimensional variables, often are called functionals.
7 Number one chooses a partner from among the 2m − 1 other numbers. This is the factor

2m − 1. Then the lowest as yet unpaired number chooses a partner from the 2m − 3 remaining
numbers, and so on.

5
corresponding to the three pairings above, we have the three terms

E Y12 Y22 = E [Y1 · Y1 · Y2 · Y2 ]

 

= E [Y1 Y1 ] E [Y2 Y2 ] + E [Y1 Y2 ] E [Y1 Y2 ] + E [Y1 Y2 ] E [Y1 Y2 ]

= 2·5+3·3+3·3
= 33 .

The last two terms in the sum are the same but correspond to the distinct
pairings {{1, 3} , {2, 4}} and {{1, 4} , {2, 3}}. All the terms on the right are given
by covariances, which illustrates the general fact that a Gaussian is determined
by its mean and covariance matrix. This example used row vectors l1 = l2 =
(1, 0) and l3 = l4 = (0, 1). For more complicated row vectors, we have the
formula (which the reader should verify)

E[l1 (Y )l2 (Y )] = l1 Cl2t ,

where C is the covariance matrix of Y . For a scalar mean zero Gaussian, all
pairings give the same contribution, so

E[X 4 ] = 3σX
4
, E[X 6 ] = 15σX
6
, E[X 8 ] = 105σX
8
, etc.

One of the proofs of the central limit theorem also explains Wick’s theorem.
If a mean zero random variable satisfies Wick’s theorem, then all its moments
are the same as the Gaussian moments so it is Gaussian8 . We take a concrete
case:
E[l1 (Zl ) · · · · · l4 (ZL )] → (Wick formula) as L → ∞.
To evaluate the sum, we use a different summation index for each factor j =
1, 2, 3, 4:
L
1 X
lj (ZL ) = √ lj (Ykj ) .
L kj =1

so that

E [l1 (ZL )l2 (ZL )l3 (ZL )l4 (ZL )]

L L L L
1 X X X X
= E[l1 (Yk1 )l2 (Yk2 )l3 (Yk3 )l4 (Yk4 )]
L2
k1 =1 k2 =1 k3 =1 k4 =1

Most of the expectations on the right side are zero. For example, if k1 6= k2 ,
k1 6= k3 , and k1 6= k4 then l1 (Yk1 ) has mean zero and is independent of the
other l(Y ) factors so E[l1 (Yk1 )l2 (Yk2 )l3 (Yk3 )l4 (Yk4 )] = 0. To get a nonzero
expectation, each term must be paired with at least one other. For example,
the pairing {{1, 3} , {2, 4}} corresponds to possibly nonzero terms with k1 = k3
and k2 = k4 : E[l1 (Yk1 )l2 (Yk2 )l3 (Yk1 )l4 (Yk2 )].
8 The justification of this statement is an easy case of the moment problem, showing that

two probability laws with the same moments are the same.

6
 If k1 6= k2 then Yk1 and Yk2 are independent, and we get a Wick type
contribution:

E[l1 (Yk1 )l2 (Yk2 )l3 (Yk1 )l4 (Yk2 )] = E[l1 (Yk1 )l3 (Yk1 )]E[l2 (Yk2 )l4 (Yk2 )]
= l1 Cl3t l2 Cl4t .

There are L(L − 1) such terms. The L remaining terms have k1 = k2 = k3 = k4 .

Altogether,

E [l1 (ZL )l2 (ZL )l3 (ZL )l4 (ZL )]

L(L − 1) n t

t

t

t

t

t

o
= l 1 Cl 2 l 3 Cl 4 + l 1 Cl 3 l 2 Cl 4 + l 1 Cl 4 l 2 Cl 3
L2
L
+ 2 E[l1 (Y )l2 (Y )l3 (Y )l4 (Y )] .
L
The first term on the right converges to the Wick formula as L → ∞ while the
second converges to zero.
This argument has drawbacks as a proof of the central limit theorem. It
requires Y to have finite moments (e.g. so that the last term above is finite).
It relies on the moment problem, which, at a minimum, is more sophisticated
than this. It does not answer questions about convergence of the distribution
of ZL other than moments. However, if you believe the central limit theorem
already and just want Wick’s formula, take Y to be Gaussian.
Here is an example of the multivariate central limit theorem that comes
up in the theory of time stepping methods for stochastic differential equations.
Suppose X ∼ N (0, 1) and we define a three component random variable Y
as Yk = (X 2 − 1, X 3 , X 4 − 3)t . The Y variances (handily computed using
Wick’s theorem) are C11 = 2, C22 = 15 and C33 = 96. The off diagonal
covariances all are zero except C13 = 12. Note that although Y1 and Y2 are
uncorrelated√(C12 = 0), they are far from independent. In fact, if we know Y1 ,
then X = ± Y1 + 1: If we know Y1 then we know everything about Y2 except
the sign. Nevertheless, if we have L independent X samples then the random
variables
L L
1 X 1 X 3
Xk2 − 1 , ZL,2 = √


ZL,1 = √ Xk
L k=1 L k=1
are nearly independent Gaussians for large L.