iconEuler Examples

The Box-Muller Method

by R. Grothmann

Euler simulates the normal distribution with the following trick due
to Box and Muller. We like to explain that trick here.

Assume, we have two [0,1]-equi-distributed variables, which we call x
and y.
>x:=linspace(0.01,1,20); y:=linspace(0,1,80)';
Now we map this into the plane, using r=sqrt(-2*log(x)) and phi=2*pi*y.
The image is the complete plane. Since we take x>0.01, we only
get part of the plane.
>z:=sqrt(-2*log(x))*exp(2pi*I*y); plot2d(z); insimg;

Box Muller Method

Now, the Box Muller method says, that both coordinates of z are
0-1-normal distributed. We check this for our chosen grid. So
we have to look at the distribution of all re(z).
>t:=redim(re(z),1,cols(z)*rows(z)); ...
>plot2d(t,distribution=10); plot2d("qnormal(x,0,1)",add=1,thickness=3); insimg;

Box Muller Method

But why is that so?

Let us take a tiny square Q in the plane at z=(t,s) with side length
(dt,ds). We want to know, how many (x,y) fall into that square under
the mapping

Box Muller Method

This is the area of the inverse of the Q under f, which again is
assymtotically equal to

Box Muller Method

So the densitity function of z in the plane is 1/det(Jf(x,y)), where
f(x,y)=z. Let us compute this.
>f1 &= sqrt(-2*log(x))*sin(2*pi*y), f2 &= sqrt(-2*log(x))*cos(2*pi*y)
                  sqrt(2) sqrt(- log(x)) sin(2 pi y)


                  sqrt(2) sqrt(- log(x)) cos(2 pi y)

>&1/trigsimp(determinant(jacobian([f1,f2],[x,y])))
                                  x
                                 ----
                                 2 pi

Scaled times 2*pi this is the following densitity distribution.
>plot3d(re(z),im(z),x); insimg;

Box Muller Method

If we start with z=(t,s), we see that

Box Muller Method

so

Box Muller Method

The density is thus equal to the two dimensional 0-1-normal
distribution. And this distribution has the property that the
coordinates are 0-1-normal distributed.

We can check this with Maxima, but it is also easy to see by hand. To
see, how many z fall into a strip of width ds at s, we integrate the
density for t=-inf to t=inf.
>&integrate(exp(-(t^2+s^2)/2)/(2*%pi),t,-inf,inf)
                                    2
                                   s
                                 - --
                                   2
                                E
                           ----------------
                           sqrt(2) sqrt(pi)

>

Examples Homepage