Chapter 2: Generating Random Variates

ST4231: Computer Intensive Statistical Methods

Semester 1, AY 2023/2024

Department of Statistics and Data Science

National University of Singapore

LI Cheng

stalic@[Link]

1 / 57
 Outline

Random Number Generators

Basics of Monte Carlo Methods

Inversion Method

Rejection Sampling

Polar Method for Bivariate Normal

2 / 57
 Random Number Generators

• The building block of simulation is the ability to generate random

numbers, where a random number represents the value of a random
variable uniformly distributed on (0,1), i.e. X ∼ Uniform(0, 1).
• Although not an emphasis of this module, we briefly introduce the
congruential method to generate pseudo random numbers.
• The idea is to construct one system using number theory which can
generate a sequence of values, which, although deterministically
generated, have all the appearances of being independent uniform
(0,1) random variables.

3 / 57
 Multiplicative Congruential Method

• Starting with an initial value x0 , called the seed, and we recursively

compute successive values xn , n ≥ 1, by letting

xn = axn−1 mod m (1)

where a and m are some given positive integers. Thus, each xn takes
value in 0, 1, · · · , m − 1.
• Here a mod m means the remainder after the integer a is divided by
the integer m. Here “mod” is a shorthand of “modulus”. For
example, 5 mod 2 = 1, 48 mod 7 = 6, 2022 mod 100 = 22.
• The quantity xn /m (called a pseudo random number) is taken as an
approximation to the value of a uniform (0,1) random variable.

4 / 57
 Three Criteria for Pseudo-random Number

• In general the constants a and m should be chosen to satisfy three

criteria:
1. For any initial seed, the resultant sequence has the “appearance” of
being a sequence of independent uniform (0,1) random variables.
2. For any initial seed, the number of variables that can be generated
before repetition begins is large.
3. The values can be computed efficiently on a digital computer.

5 / 57
 Guideline

• A guideline that appears to be of help in satisfying the above three

conditions is that m should be chosen to be a large prime number
that can be fitted to the computer word size.
• For a 32-bit word machine, it has been shown that the choice of
m = 231 − 1 and a = 75 = 16807 result in desirable properties.
• Three main generators
• Mixed Congruential xn = (axn−1 + c) mod m
−1
• Inversive Congruential xn = (axn−1 + c) mod m
• Nonlinear Congruential xn = (dxn−1
2
+ axn−1 + c) mod m

6 / 57
 Other Generators

• Other random number generators include Feedback Shift Register

methods, Coupled Generators, Portable Generators, and the
methods based on Chaotic Systems.
• Most computer languages already have a built-in (standard uniform)
random number generator.
• Based on these random number generators, we can further simulate
non-uniform random variates.

7 / 57
 Our Starting Point

(i) Throughout this course, we will always assume that we can generate
random numbers (or sample) from Uniform(0, 1).

(ii) For the rest of this chapter, we assume that we can only sample
from Uniform(0, 1).
In other words, all random variates from other distributions (e.g.
Exp(λ), N(0, 1), Gamma(a, b), . . .) must be some transformed
version of Uniform(0, 1) random variates.

8 / 57
 Outline

Random Number Generators

Basics of Monte Carlo Methods

Inversion Method

Rejection Sampling

Polar Method for Bivariate Normal

9 / 57
 Evaluating an Integral

• Numerical integration methods based on interpolation, e.g. using

Riemann sums.1

• Problem: Curse of dimensionality!

1
Figures from Wolfram MathWorld [Link] and
Wolfram Mathematica [Link]
10 / 57
 Monte Carlo Integration

• This is one of the earliest applications of random numbers.

• Suppose we want to compute
Z 1 Z 1
θ= g (x)dx = g (x) · 1dx.
0 0

• If U ∼ Uniform(0, 1), then we can write θ as

θ = E[g (U)]

• If U1 , U2 , · · · , Uk are independent Uniform(0, 1) random variables,

then we have that g (U1 ), · · · , g (Uk ) are i.i.d. random variables.

11 / 57
 Monte Carlo Integration
• By strong law of large numbers, it follows that: if
Z 1
|g (x)|dx < ∞, then with probability 1,
0

k
1X
g (Ui ) → E[g (U)] = θ as k → ∞.
k
i=1

• Hence, we can generate a large number of random numbers ui from

Uniform(0, 1), and approximate θ by the average value of g (ui ), i.e.,
k
1X
θ̂ = g (ui )
k
i=1

• This approach to approximating integrals is called the Monte Carlo

integration.

12 / 57
 Multi-dimensional Integral
Z 1Z 1 Z 1
θ= ··· g (x1 , x2 , · · · , xn )dx1 dx2 · · · dxn
0 0 0
Generate k independent sets, each consisting of n independent
Uniform(0, 1) random variates
(1) (1)
u1 ··· un
(2) (2)
u1 ··· un
..
.
(k) (k)
u1 ··· un

(i) (i)
Since g (U1 , · · · , Un ), i = 1, 2, · · · , k are i.i.d, we have that with
probability 1,
k
1X (i) (i)
θ̂ = g (U1 , · · · , Un ) → θ
k
i=1

13 / 57
 Example: Estimation of π

• Suppose the random vector (X , Y ) is uniformly distributed in the

square of area 4 centered at the origin (0,0).
• We have the following relationship
π
P(X 2 + Y 2 ≤ 1) =
4
• Hence, if we generate a large number of random points in the
square, the proportion of points that fall within the circle will be
approximated π/4.
• Question: What is the g function and what is θ here?

14 / 57
 Example: Estimation of π

• We will come back to Monte Carlo integration in the next chapter.

• In this chapter, we first introduce several methods to draw samples

from a distribution. After we address this problem, we can use the
samples drawn to do Monte Carlo integration.

15 / 57
 Outline

Random Number Generators

Basics of Monte Carlo Methods

Inversion Method

Rejection Sampling

Polar Method for Bivariate Normal

16 / 57
 Overview

• Goal: Generate a random variable X (or a random sample of

{X1 , . . . , Xn }) from a given distribution with cdf F .
• Tool: One can generate U ∼ Uniform(0, 1).

• Methods: We will introduce several methods to achieve this goal,

including inversion method, rejection sampling (accept-reject
method), as well as methods for multivariate distributions.

17 / 57
 The Idea of Inversion Method
• Goal: Generate a random variable X (or a random sample of
{X1 , . . . , Xn }) from a given distribution with cdf F .
• Tool: One can generate U ∼ Uniform(0, 1).

• Inversion Method: Find a transformation T : [0, 1] → R such that

X = T (U).
• If we can find T , we are done!

• For Monte Carlo integration: for any function φ that satisfies

E[|φ(X )|] < ∞, we can use the fact that (why?) with probability 1,
as n → ∞,
φ(T (U1 )) + · · · + φ(T (Un ))
→ E[φ(X )].
n

18 / 57
 Inversion Method

• What is the transform T ?

• Suppose X = T (U), then assume that T is a bijection, we have

F (x) = P{X ≤ x} = P{T (U) ≤ x}

= P{U ≤ T −1 (x)} = T −1 (x).

Therefore, we can set T (x) = F −1 (x).

Theorem
Let F be the cdf of a random variable, and let U be a random variable
with the Uniform(0, 1) distribution. Then F −1 (U) ∼ F , i.e. the cdf of
the random variable X = F −1 (U) is F .

19 / 57
 Inversion Method

About the quantile function F −1 :

• When F : R → [0, 1] is continuous and strictly increasing (as for
most continuous random variables), then F −1 : [0, 1] → R is also
continuous and strictly increasing.
• In general, we only know that F is right-continuous and
non-decreasing. F could be discontinuous (as for most discrete
random variables). Then we can define F −1 as the generalized
inverse:

F −1 (u) = inf{z ∈ R : F (z) ≥ u}, u ∈ [0, 1].

20 / 57
 Discrete Random Number Generators

• For discrete distributions, the inversion method by using the quantile

function is equivalent to the following method.
• Suppose that we want to generate the value of a discrete random
variable X having probability mass function
X
P(X = xj ) = pj , j = 0, 1, . . . , pj = 1
j

• Algorithm
• Generate U ∼ Uniform(0, 1).
• if U < p0 , set X = x0 and stop.
• if U < p0 + p1 , set X = x1 and stop.
• if U < p0 + p1 + p2 , set X = x2 and stop.
..
• .

21 / 57
 Discrete Random Number Generators

• If the xi , i ≥ 0, are ordered so that x0 < x1 < · · · and if we let F

Xk
denote the distribution function of X , then F (xk ) = pi and so
i=0

X will be equal to xj if F (xj−1 ) ≤ U < F (xj )

• In other words, after generating U ∼ Uniform(0, 1), we determine

the value of X by finding the interval [F (xj−1 ), F (xj )) in which U
lies (or, equivalently, by finding the inverse of F (U) and see which
[xj−1 , xj )-interval U lies in).

22 / 57
 Discrete Random Number Generators

Example: A simple discrete random number generator

• If we want to simulate a random variable X such that

p1 = 0.2, p2 = 0.15, p3 = 0.25, p4 = 0.4,

where pj = P(X = j).

• We could generate U ∼ Uniform(0, 1) and do the following:
• if U < 0.2 set X = 1 and stop.
• if U < 0.35 set X = 2 and stop.
• if U < 0.6 set X = 3 and stop.
• Otherwise set X = 4.

23 / 57
 Discrete Random Number Generators

• However, a more efficient procedure is the following (why?)

• if U < 0.4 set X = 4 and stop.
• if U < 0.65 set X = 3 and stop.
• if U < 0.85 set X = 1 and stop.
• Otherwise set X = 2.

24 / 57
 Continuous Random Variates

Example 1. A simple continuous random number generator:

• Suppose we want to generate X from the distribution

F (x) = x n , 0<x <1

• Let x = F −1 (u), then

x = u 1/n , or, equivalently, u = F (x) = x n

• Hence we have the following algorithm for generating a random
variable from the cdf F .
• Generate a random number U ∼ Uniform(0, 1).
• Set X = U 1/n .

25 / 57
 Continuous Random Variates

Example 2. Exponential random number generator:

• If X ∼ Exp(λ), then its distribution function is given by

F (x) = 1 − e −λx

• Let x = F −1 (u), we have

u = F (x) = 1 − e −λx

or, taking logarithms,

1
x = − log(1 − u)
λ

26 / 57
 Continuous Random Variates
• Hence, we can generate X ∼ Exp(λ) by first generating a random
number U ∼ Uniform(0, 1) and then setting
1
X = F −1 (U) = − log(1 − U)
λ
• Since 1 − U ∼ Uniform(0, 1), we can also simply use
1
X = − log(U)
λ

2
2
Figure from wikipedia. 27 / 57
 Continuous Random Variates

Example 3. Gamma random number generator:

• Suppose we want to generate the value of a Gamma(α, λ) random
variable.
Z x α α−1 −λy Z λx α−1 −t
λ y e t e
F (x) = dy = dt
0 Γ(α) 0 Γ(α)

• The inverse F −1 does not have a closed form, so we can not use the
inversion method here.
• When α is an integer α = n, using the result that a Gamma(n, λ)
random variable X can be regarded as being the sum of n
independent Exp(λ) r.v.’s, we can make use of Example 2 to
generate X .

28 / 57
 Continuous Random Variates

• Based on the above idea, we can generate a Gamma(n, λ) random

variable by first generating n random numbers
U1 , . . . , Un ∼ Uniform(0, 1), and then setting

1 1
X = − log U1 − · · · − log Un
λ λ
1
= − log(U1 · · · Un )
λ
n
X
• Where the use of the identity log(xi ) = log(x1 · · · xn ) is
i=1
computationally time saving in that it requires only one rather than
n logarithm computations.

29 / 57
 Outline

Random Number Generators

Basics of Monte Carlo Methods

Inversion Method

Rejection Sampling

Polar Method for Bivariate Normal

30 / 57
 Fundamental Theorem of Simulation

• We consider the general problem of sampling X ∼ F or X ∼ f ,

where F is a cdf and f is a pdf.
• The following theorem (Theorem 2.15, Chapter 2 of Robert and
Casella) is fundamentally important.

Theorem
If X is a random variable with pdf f (x), then simulating X is equivalent
to simulating a pair of random variables (X , U) jointly from

(X , U) ∼ Uniform {(x, u) : 0 < u < f (x)} .

31 / 57
 Fundamental Theorem of Simulation
Proof (not examinable): Let S = {(x, u) : 0 < u < f (x)} and
C = the area of S (normalizing constant). If (X , U) is uniformly
distributed in S, then the cdf of X is

1 x
Z f (t)
1 x
Z Z Z
1
P(X ≤ x) = dudt = dudt = f (t)dt.
S C C −∞ 0 C −∞

1
Taking x → +∞, we have 1 = or C = 1. This implies that f (x) is the
C
density of X . 2
• In practice, it is often difficult to simulate uniform r.v. on S, since
the shape of the region S depends on f (x), which is often difficult
to sample from.
• Instead, we can use a “proxy” distribution g (x) that is easy to
sample from, and then perform rejection sampling.

32 / 57
 Rejection Sampling
The following plot illustrates rejection sampling for a mixture distribution
0.5Beta(2, 6) + 0.5Beta(4, 2).

33 / 57
 Rejection Sampling
• Goal: To generate a d-dimensional random vector X from a known
target pdf f (x ) on Rd .
• Idea (Fundamental Theorem of Simulation): If one can generate
(X , U) uniformly distributed on the domain
S = {(x , u) : 0 ≤ u ≤ f (x )}, then X has density f (x ).

34 / 57
 Rejection Sampling

When we can use Rejection Sampling:

• Assume that g (x ) is another pdf on Rd satisfying two conditions:
• We already know how to generate random vector Y from g (x );
• There is a positive constant M such that f (x ) ≤ Mg (x ) for every
x ∈ Rd .
• Question: How do we find M?.
• Example: We use rejection sampling to sample from N(0, 1) (this is
f (x); the plot is in the previous slide).
• Assume that we only know how to sample from a Cauchy distribution;
• We can set M = 3, if we use Cauchy(0, 2) as g (x).

35 / 57
 Rejection Sampling Algorithm

(1) Generate Y ∼ g ;
(2) Generate U ∼ Uniform[0, Mg (Y )], where M satisfies that
f (x ) ≤ Mg (x ) for all x ;
(3) If U ≤ f (Y ), then accept: set X = Y and stop. Otherwise, reject:
return to step (1).

In practice, the following version is more commonly used with the

constant M satisfying f (x ) ≤ Mg (x ) for all x .
(1) Generate Y ∼ g ;
(2) Generate U ∼ Uniform[0, 1];
f (Y )
(3) If U ≤ , then accept: set X = Y and stop. Otherwise,
Mg (Y )
reject: return to step (1).

36 / 57
 Rejection Sampling Algorithm

Theorem
If X is generated via Steps 1-3 of the rejection sampling method in the
previous slide, then X has pdf f (x ).

Proof (not examinable): The pair (Y , U) generated from steps (1)-(2) satisfies Y ∼g
and U | Y ∼ Uniform[0, Mg (Y )]. Therefore, by step (3), for any set A,

P(X ∈ A, U ≤ f (Y )) P(Y ∈ A, U ≤ f (Y ))
P(X ∈ A | U ≤ f (Y )) = =
P(U ≤ f (Y )) P(U ≤ f (Y ))
du g (y )dy · g (y )dy
R  R f (y ) 1  R f (y )
A 0 Mg (y ) A Mg (y )
= R R = R
f (y )
du g (y )dy
f (y )
· g (y )dy

1
Rd 0 Mg (y ) Rd Mg (y )

f (y )dy
R Z Z
= RA = f (y )dy (since f (y )dy = 1)
Rd
f (y )dy A Rd

Since this relation holds for any set A ⊆ Rd , f is the density of X. 2

37 / 57
 Efficiency of Rejection Sampling
• For each proposal (Y , U) obtained via steps (1) and (2).

area under f (x ) 1
P{(Y , U) is accepted} = = .
area under Mg (x ) M
• Therefore, the expected number of proposals needed to sample one
single X from f is M. In fact, the number of proposals needed has
the geometric distribution with parameter 1/M.
• Thus to improve the efficiency of rejection sampling, we would like
to choose g (x ) so that M is small.
• Theoretically, if we use a proposal g ≡ f , then M = 1. However,
making M = 1 is practically impossible as it defeats the purpose of
rejection sampling.
f (x )
• Clearly, taking M = sup gives the optimal M.
x ∈Rd g (x )

38 / 57
 Example 1
Use the rejection sampling method to generate a random variable with
pdf
f (x) = 20x(1 − x)3 , 0 < x < 1

• Step 1: Specify the proposal pdf g (x). Since f (x) is the pdf of
Beta(2, 4) (why?), X is defined in the interval (0,1). Consider
Uniform(0, 1) with pdf
g (x) = 1, 0 < x < 1.
• Step 2: Find M. To determine the constant M , we maximize the
ratio
f (x)
= 20x(1 − x)3 , for 0 < x < 1.
g (x)
d
(20x(1 − x)3 ) = 20[(1 − x)3 − 3x(1 − x)2 ] = 0
dx
=⇒ x = 1/4 or 1.
39 / 57
 Example 1 (cont’d)

• A quick check shows x = 1/4 attains the maximum of f (x)/g (x).

We set M = gf (x)
(x) = 20(0.25)(1 − 0.25)3 = 135
64 .
x=1/4

• Step 3: Specify the rejection function. From Step 2, we have

f (x) ≤ Mg (x) = 135
64 g (x). Hence,

f (x) 256
= x(1 − x)3 .
Mg (x) 27

40 / 57
 Example 1 (cont’d)
• Step 4: Write down the algorithm:
• Generate random number Y and U from Uniform(0, 1).
256
• If U ≤ Y (1 − Y )3 , stop and set X = Y . Otherwise, return to
27
step 1.
• The average number of times that step 1 will be performed is
c = 135/64 ≈ 2.11.

41 / 57
 Example 2

• Suppose we want to generate a random variable from

Gamma(1.5, 1), with pdf

f (x) = Kx 1/2 e −x , x > 0

√
where K = 1/Γ(1.5) = 2/ π.
• Because such a random variable is concentrated on the positive axis
and has mean 1.5, it is natural to try the rejection technique with an
exponential random variable with the same mean.
• We can let
2
g (x) = e −2x/3 , x > 0.
3

42 / 57
 Example 2 (cont’d)
• We have
f (x) 3
= Kx 1/2 e −x/3
g (x) 2
f (x) 33/2
• Maximize this ratio. We get M = sup = , and
x>0 g (x) (2πe)1/2

f (x)
= (2e/3)1/2 x 1/2 e −x/3 .
Mg (x)
• The rejection sampling algorithm is as follows:
1. Generate a random number U1 ∼ Uniform(0, 1), and set
3
Y = − log U1 ; (why?)
2
2. Generate a random number U2 ∼ Uniform(0, 1);
3. If U2 < (2eY /3)1/2 e −Y /3 , then set X = Y and stop; otherwise,
return to step 1.

43 / 57
Rejection Sampling: Unknown Normalizing Constant

• Question: Can we use rejection sampling when f (x ) is known only

up to a normalizing constant?
• For example,  we
pwould like to sample from
f (x) ∝ exp − x 2 + 1/3 on x > 0. This means that
 p 
f (x) = c · exp − x 2 + 1/3 on x > 0. However, c is unknown
and usually difficult to calculate.
• In general, we can assume that f (x ) = c fe(x ), where fe(x ) is known,
and c normalizes fe(x ) into a proper density.
• Does rejection sampling still work?

• The answer is yes. In fact, we do not even need to know c!

44 / 57
Rejection Sampling: Unknown Normalizing Constant
• Recall that in rejection sampling, the key step is to find a bound M,
such that f (x )/g (x ) ≤ M for all x .
• In the algorithm, the key steps are
• Generate Y∼ g (Y ), U ∼ Uniform[0, 1];
f (Y )
• Decide if U ≤
Mg (Y )
or not.

• Now with only fe(x ) known, we cannot find M. However, we can

still find a constant Me to upper bound the ratio of fe(x )/g (x ), e.g.
M = sup f (x )/g (x ).
e e
x

• Then we can simply modify the two steps above into

• Generate Y∼ g (Y ), U ∼ Uniform[0, 1];
fe(Y )
• Decide if U ≤ or not.
e (Y )
Mg

45 / 57
 Rejection Sampling: Unknown Normalizing Constant

Suppose f (x ) = c fe(x ) where fe(x ) is known and c is unknown. M e

satisfies that f (x ) ≤ Mg (x ) for all x . Then the rejection sampling
e e
algorithm is given as follows:

(1) Generate Y ∼ g ;
(2) Generate U ∼ Uniform[0, 1];
fe(Y )
(3) If U ≤ , then accept: set X = Y and stop. Otherwise,
e (Y )
Mg
reject: return to step (1).

Corollary
If X is generated via steps 1-3 of the rejection sampling method above,
then X has pdf f (x ).

46 / 57
 Example 3

• Suppose we want to generate a random variable from the pdf

 p 
f (x) ∝ exp − x 2 + 1/3 , x > 0.

• Since the domain is (0, +∞), we can use the exponential

distribution as proposal, e.g. Exp(1/3).
• We let
1
g (x) = e −x/3 , x > 0.
3

47 / 57
 Example 3 (cont’d)

 p 
• Clearly, f˜(x)/g (x) = 3 exp − x 2 + 1/3 + x/3 . The function is
maximized when x → +∞ (explain), and f˜(x)/g (x) ≤ 3.
 p 
• Therefore, we can set M e = 3, fe(x) = exp − x 2 + 1/3 + x/3 .
Mg (x)
e

• The rejection sampling algorithm is as follows:

1. Generate a random number U1 ∼ Uniform(0, 1), and set
Y = −3 log U1 ; (why?)
 p number U2 ∼ Uniform(0,
2. Generate a random  1);
3. If U2 < exp − Y 2 + 1/3 + Y /3 , then set X = Y and stop;
otherwise, return to step 1.

48 / 57
 Outline

Random Number Generators

Basics of Monte Carlo Methods

Inversion Method

Rejection Sampling

Polar Method for Bivariate Normal

49 / 57
 The Polar Method for Bivariate Normal
We first consider generating independent bivariate normal.
• Let X and Y be independent standard normal random variables. Let
R and θ denote the polar coordinates of the vector (X , Y ). That is

Y
S = R 2 = X 2 + Y 2, tan θ = .
X
• Since X and Y are independent, we have the joint density

1 2 1 2 1 −(x 2 +y 2 )/2
f (x, y ) = √ e −x /2 √ e −y /2 = e
2π 2π 2π
• Using the change-of-variable formula from (X , Y ) to (S, θ), we have
that
1 1 −s/2
f (s, θ) = e , 0 < s < ∞, 0 < θ < 2π.
2 2π

50 / 57
 The Polar Method for Bivariate Normal

• It can be seen that this joint density of (S, θ) is equal to the product
1
of an exponential density e −s/2 and the uniform density on (0, 2π).
2

• Therefore, S and θ are independent, with S = R 2 ∼ Exp(1/2) and

θ ∼ Uniform(0, 2π). (By definition, S has a χ22 distribution, which is
Exp(1/2).)
• To generate (X , Y ), we can to first generate their polar coordinates
(R, θ) and then transform back.

51 / 57
 The Polar Method for Bivariate Normal

Box-Muller Algorithm v1:

1. Generate random numbers U1 ∼ Uniform(0, 1) and
U2 ∼ Uniform(0, 1).
p
2. Set R = −2 log(U1 ) (why?), θ = 2πU2 .
3. Set
p
X = −2 log U1 cos(2πU2 )
p
Y = −2 log U1 sin(2πU2 )

This algorithm is not the most efficient, because evaluating the

trigonometric functions can be expensive.

52 / 57
 The Polar Method for Bivariate Normal
• Suppose that (V1 , V2 ) is uniformly distributed in the disk centered
at (0, 0) with radius 1. Then for the random angle
θ ∼ Uniform(0, 2π), we can write

V1 V2
cos θ = q , sin θ = q .
V12 + V22 V12 + V22

53 / 57
 The Polar Method for Bivariate Normal

• Following the Box-Muller transformation, we can generate

independent unit normals (X , Y ) as follows,

V1
X = (−2 log(U))1/2
(V12 + V22 )1/2
V2
Y = (−2 log(U))1/2 2 (2)
(V1 + V22 )1/2

where U ∼ Uniform(0, 1) and (V1 , V2 ) is uniformly distributed in the

disk centered at (0, 0) with radius 1.
• We can use the rejection sampling method to draw (V1 , V2 ) (see
slide 16).

54 / 57
 The Polar Method for Bivariate Normal
• We need the following important result:
S = V12 + V22 is uniformly distributed on (0,1) and is independent of
the random angle θ (associated to the point (V1 , V2 )).
• Therefore, we can use S as the random number U needed in
Equation (2).
√
• By viewing S as the “radial coordinate”, we obtain that
p
X = −2 log(S)/S · V1
p
Y = −2 log(S)/S · V2

are independent unit normals (explain).

• We just need to ensure that (V1 , V2 ) is uniformly distributed in the
disk centered at (0, 0) with radius 1, and set S = V12 + V22 .

55 / 57
 The Polar Method for Bivariate Normal

Box-Muller Algorithm v2:

1. Generate random numbers U1 and U2 from Uniform(0, 1).
2. Set V1 = 2U1 − 1, V2 = 2U2 − 1, S = V12 + V22 .
3. If S > 1 return to Step 1 (rejection sampling).
4. Return the independent unit normals
p
X = −2 log(S)/S · V1
p
Y = −2 log(S)/S · V2

Step 3 is to ensure that (V1 , V2 ) is uniformly distributed in the disk

centered at (0, 0) with radius 1.

56 / 57
 General Multivariate Normal

• d−dimensional normal with mean µ and covariance matrix Σ:

 
1 1
f (x ; µ, Σ) = exp − (x − µ) Σ (x − µ)
⊤ −1
(2π)d/2 |Σ|1/2 2

where |Σ| is the determinant of Σ.

• The most commonly used method to simulate from this distribution
is to first generate Z = (Z1 , . . . , Zd )⊤ with Z1 , . . . , Zd i.i.d. N(0, 1),
and then set
X = LZ + µ,
where L satisfies LL⊤ = Σ. Usually L is taken as the Cholesky
factor, a lower triangular matrix with positive diagonal entries.
• Then X ∼ N(µ, Σ).

57 / 57