MA 323 Monte Carlo Simulation Lecture 04

1 Generation from Univariate Normal Distribution

Recall that the PDF and CDF of univariate normal distribution with mean µ (µ ∈ R) and variance
σ 2 (σ > 0) are

1 x−µ 1 2 2
ϕµ,σ (x) = ϕ( )= √ e−(x−µ) /2σ for x ∈ R
σ σ 2πσ
and
  x
x−µ
Z
1 2 /2σ 2
Φµ,σ (x) = Φ =√ e−(t−µ) dt for x ∈ R,
σ 2πσ −∞

respectively, where
1 x2
ϕ(x) = √ e− 2 for x ∈ R
2π
and
Z x
1 t2
Φ(x) = √ e− 2 dt for x ∈ R.
2π −∞

Note that ϕ(·) and Φ(·) are PDF and CDF of standard normal distribution. The word “stan-
dard” indicates that the mean is 0 and the variance is 1. If Z ∼ N (0, 1), then µ + σZ ∼ N (µ, σ 2 ).
Thus, given a method for generating samples Z1 , Z2 , . . . from the standard normal distribution, we
can generate samples X1 , X2 , . . . from N (µ, σ 2 ) by setting Xi = µ + σZi . It, therefore, suffices to
consider methods for sampling from N (0, 1).
We now discuss algorithms for generating univariate normal distribution. We assume the avail-
ability of a sequence U1 , U2 , . . . of independent random variables uniformly distributed on the unit
interval (0, 1) and consider methods for transforming these uniform random variables to normally
distributed random variables.

1.1 Using Acceptance Rejection Method

x2
To obtain an upper bound of ϕ(x), we need to find an lower bound of 2
. Notice that

1 x2 1
(|x| − 1)2 = − |x| + ≥ 0.
2 2 2
x2 1
Thus, ≥ |x| − for all x ∈ R. This shows that, for all x ∈ R,
2 2
1 x2 1 1
ϕ(x) = √ e− 2 ≤ √ e−|x|+ 2 = cg(x),
2π 2π
 r
2e
where c = and g(x) = 12 e−|x| for all x ∈ R. Therefore, we can use acceptance rejection
π
method if we can generate from g easily. To generate from g(·) we can use the following result. The
PDF of X = ZY is g(·) if Y and Z are independent random variables such that Y ∼ Exp(1) and
P (Z = 1) = P (Z = −1) = 12 . To see it, notice that the CDF of X is

P (X ≤ x) = P (X ≤ x|Z = 1) P (Z = 1) + P (X ≤ x|Z = −1) P (Z = −1)

1 1
= P (ZY ≤ x|Z = 1) + P (ZY ≤ x|Z = −1)
2 2
1 1
= P (Y ≤ x) + P (Y ≥ −x)
2
( 2
1 x
e if x < 0
= 2 1 −x
1 − 2e if x ≥ 0.

This shows that the PDF of X is same as g(·). This result tells us the method of generation from
g. First, generate Y from exponential distribution with mean 1 and then assign a random
sign to it. Therefore, we have the Algorithm 1.

Algorithm 1 Generation of N (0, 1) by rejecting from Laplace PDF

1: repeat
2: generate U from U (0, 1) ▷ Will be used to generate from Exp(1)
3: X ← − ln(U ) ▷ X ∼ Exp(1)
4: generate V from U (0, 1) ▷ Will be used to generate random sign
1
5: if V < 2 then ▷ Assigning signs with equal probability
6: X ← −X
7: end if
8: generate W from U (0, 1) ▷ Will be used to implement acceptance rejection step
1 X2
9: until W e 2 −|X| ≤ e− 2

10: return X

Note that the condition at Line 9 of the Algorithm 1 does not depend on the sign of X.
1 X2
Thus, we may assign the sign if the candidate is accepted. Moreover, W e 2 −|X| ≤ e− 2 is equivalent
to −2 ln W ≥ (|X| − 1)2 . Thus, Algorithm 1 can be modified to obtain a better algorithm given
below.

Algorithm 2 Generation of N (0, 1) by rejecting from Laplace PDF

1: repeat
2: generate U from U (0, 1) ▷ Will be used to generate from Exp(1)
3: X ← − ln(U ) ▷ X ∼ Exp(1)
4: generate W from U (0, 1) ▷ Needed to implement acceptance rejection step and assign sign
2
5: until −2 ln (W ) ≥ (X − 1)
6: if W < 21 then
7: X ← −X
8: end if
9: return X
2 Box-Muller Method
Perhaps the simplest method to implement (though not the fastest or necessarily the most
convenient) is the one by Box-Muller. This algorithm generates a sample from a bivariate standard
normal distribution, each component of which is thus a univariate standard normal. The algorithm
is based on the following result.
i.i.d.
Theorem 1. Let U1 , U2 ∼ U (0, 1). Define
p p
Z1 = −2 ln U1 cos (2πU2 ) and Z2 = −2 ln U1 sin (2πU2 ) .
i.i.d.
Then Z1 , Z2 ∼ N (0, 1).
Proof. Note that

∂z1 cos (2πu2 )

=− √ ,
∂u1 u1 −2 ln u1
∂z2 sin (2πu2 )
=− √ ,
∂u1 u1 −2 ln u1
∂z1 p
= −2π sin (2πu2 ) −2 ln u1 ,
∂u2
∂z2 p
= 2π cos (2πu2 ) −2 ln u1 .
∂u2
Thus, the Jacobian of the transformation is
 ∂z1 ∂z2 
2π
J = det ∂u 1
∂z1
∂u1
∂z2 =− .
∂u2 ∂u2 u1

Therefore, the absolute value of the Jacobian of inverse transformation is

1 u1 1 − 12 (z12 +z22 )
= = e .
|J| 2π 2π

and hence the JPDF of (Z1 , Z2 ) is

1 − 21 (z12 +z22 )
fZ1 , Z2 (z1 , z2 ) = e = ϕ(z1 )ϕ(z2 ) for z1 ∈ R, z2 ∈ R.
2π
This proves the result.
Now, based on the theorem, we have the Algorithm 3 to generate a pair of independent standard
normal random numbers from a pair of independent U (0, 1) random numbers.

Algorithm 3 Box-Muller Method to generate N (0, 1) random numbers

1: generate
√ U1 and U2 from U (0, 1)
2: R ← −2 ln U1
3: θ ← 2πU2
4: Z1 ← R cos(θ)
5: Z2 ← R sin(θ)
6: return (Z1 , Z2 ).
3 Marsaglia and Bray Method
Marsaglia and Bray developed a modification of the Box-Muller method that reduces computing
time by avoiding evaluation of the “cos” and “sin” functions. The Marsaglia and Bray method
instead uses acceptance rejection method to sample points uniformly in the unit disc and
then transforms these points to normal variables. The algorithm is as follows:

Algorithm 4 Marsaglia and Bray Method to generate N (0, 1) random numbers

1: repeat
2: generate U1 and U2 from U (0, 1)
3: U1 ← 2U1 − 1 and U2 ← 2U2 − 1 ▷ Ui ∼ U (−1, 1)
4: until U12 + U22 ≤ 1 ▷ (U1 , U2 ) is uniformly distributed on the disc of radius 1 centered at the
origin
h i1 h i1
−2 ln(U12 +U22 ) 2 −2 ln(U12 +U22 ) 2
5: Z1 ← U1 U12 +U22
and Z 2 ← U2 U12 +U22
▷ Zi ∼ N (0, 1) and they are independent
6: return (Z1 , Z2 ).