Computational Efficiency of Polar and Box Muller Method Using Monte Carlo

Lecture Notes on Monte Carlo Methods
Fall Semester, 2005
Courant Institute of Mathematical Sciences, NYU
Jonathan Goodman, [e-mail]

Chapter 2: Simple Sampling of Gaussians. created August 26, 2005
Generating univariate or multivariate Gaussian random variables is simple and fast. There should be no reason ever to use approximate methods based, for example, on the Central limit theorem.

1

Box Muller

It would be nice to get a standard normal from a standard uniform by inverting the distribution function, but there is no closed form formula for this distribution
2
x function N (x) = P (X < x) = √1 π −∞ e−x /2 dx . The Box Muller method is a
2
brilliant trick to overcome this by producing two independent standard normals from two independent uniforms. It is based on the familiar trick for calculating
∞

2

e−x

I=

/2

dx .

−∞

This cannot be calculated by “integration” – the indefinite integral does not have an algebraic expression in terms of elementary functions (exponentials, logs, trig functions). However,
∞

2

e−x

I2 =

∞
/2

e−y

dx

−∞

2

∞
/2

∞

−∞

2

e−(x

dy =
−∞

+y 2 )/2

dxdy .

−∞

The last integral can be calculated using polar coordinates x = r cos(θ), y = r sin(θ) with area element dxdy = rdrdθ, so that
2π

I2 =

r = 0∞

e−r

2

/2

rdrdθ = 2π

r = 0∞ e−r

2

/2

rdr .

θ =0

Unlike the original x integral, this r integral is elementary. The substitution s = r2 /2 gives ds = rdr and
∞

e−s ds = 2π .

I 2 = 2π s=0 The Box Muller algorithm is a probabilistic interpretation of this trick. If
(X, Y ) is a pair of independent standard normals, then the probability density is a product:
2
2
1
1 −(x2 +y2 )/2
1
e
.
f (x, y ) = √ e−x /2 · √ e−y /2 =
2π
2π
2π

1

Since this density is radially symmetric, it is natural to consider the polar coordinate random variables (R, Θ) defined by 0 ≤ Θ < 2π and X = R