Fall Semester, 2005

Courant Institute of Mathematical Sciences, NYU

Jonathan Goodman, goodman@cims.nyu.edu

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 indeﬁnite 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, Θ) deﬁned by 0 ≤ Θ < 2π and X = R cos(Θ), and Y = R sin(Θ). Clearly Θ is uniformly distributed in the interval [0, 2π ] and may be sampled using

Θ = 2πU1 .

Unlike the original distribution function N (x), there is a simple expression for the R distribution function:

2π

r

G(R) = P (R ≤ r) =

r =0

Θ=0

r

1 −r 2 /2

e

rdrdθ =

2π

e−r

2

/2

rdr .

r =0

The same change of variable r 2 /2 = s, r dr = ds (so that r = r when s = r2 /2) allows us to calculate

r 2 /2

e−s dx = 1 − e−r

G(r) =

2

/2

.

s=0

Therefore, we may sample R by solving the distribution function equation1 G(R) = 1 − e−R

2

/2

= 1 − U2 ,

whose solution is R = −2 ln(U2 ). Altogether, the Box Muller method takes independent standard uniform random variables U1 and U2 and produces independent standard normals X and Y using the formulas Θ = 2πU1 , R =

−2 ln(U2 ) , X = R cos(Θ) , Y = R sin(Θ) .

(1)

It may seem odd that X and Y in (13) are independent given that they use the same R and Θ. Not only does our algebra shows that this is true, but we can test the independence computationally, and it will be conﬁrmed. Part of this method was generating a point “at random” on the unit circle. We suggested doing this by choosing Θ uniformly in the interval [0, 2π ] then taking the point on the circle to be (cos(Θ), sin(Θ)). This has the possible drawback that the computer must evaluate the sine and cosine functions. Another way to do this2 is to choose a point uniformly in the 2 × 2 square −1 ≤ x ≤ 1, 1 ≤ y ≤ 1 then rejecting it if it falls outside the unit circle. The ﬁrst accepted point will be uniformly distributed in the unit disk x2 + y 2 ≤ 1, so its angle will be random and uniformly distributed. The ﬁnal step is to get a point on the unit circle x2 + y 2 = 1 by dividing by the length. The methods have equal accuracy (both are exact in exact arithmetic). What distinguishes them is computer performance (a topic discussed more in a later lecture, hopefully). The rejection method, with an...