Sampling From the Normal Distribution Using the Box-Muller Transform
The Normal Distribution is the workhorse of many common statistical analyses and being able to draw samples from this distribution lies at the heart of many statistical/machine learning algorithms. There have been a number of methods developed to sample from the Normal distribution including Inverse Transform Sampling, the Ziggurat Algorithm, and the Ratio Method (a rejection sampling technique). In this post we will focus on an elegant method called the Box-Muller transform.
A quick review of Cartesian and polar coordinates.
Before we can talk about using the Box-Muller transform, let’s refresh our understanding of the relationship between Cartesian and polar coordinates. You may remember from geometry that if x and y are two points in the Cartesian plane they can be represented in polar coordinates with a radius and an angle using the following relationships:
, and therefore
Notice that if and , then we map out values contained in the unit circle, as shown in the figure below. Also note that random variables in such a circle can be generated by transforming values sampled from the uniform distribution. Specifically, radii can be sampled from and angle can be sampled from . A similar mechanism (i.e. drawing points in a circle using uniform variables) is at the heart of the Box-Muller transform for sampling Normal random variables.
Drawing Normally-distributed samples with the Box-Muller transform
Ok, now that we’ve discussed how Cartesian coordinates are represented in polar coordinates, let’s move on to how we can use this relationship to generate random variables. Box-Muller sampling is based on representing the joint distribution of two independent standard Normal random Cartesian variables and
in polar coordinates. The joint distribution (which is circular-symmetric) is:
If we notice that the term in the numerator of the exponent is equal to (as above) we can make the connection between the the Cartesian representation of the joint Normal distribution and its polar representation:
which is the product of two density functions, an exponential distribution over squared radii:
and a uniform distribution over angles:
just like those mentioned above when generating points on the unit circle. Now, if we make another connection between the exponential distribution and the uniform distribution, namely that:
This gives us a way to generate points from the joint Gaussian distribution by sampling from two independent uniform distributions, one for and another for , and transforming them into Cartesian coordinates via the relationships above. In detail, the procedure goes as follows:
- Transform the variables into radius and angle representation , and
- Transform radius and angle into Cartesian coordinates:
What results are two independent Normal random variables, and . A MATLAB implementation of the Box-Muller algorithm is shown below:
% NORMAL SAMPLES USING BOX-MUELLER METHOD % DRAW SAMPLES FROM PROPOSAL DISTRIBUTION u = rand(2,100000); r = sqrt(-2*log(u(1,:))); theta = 2*pi*u(2,:); x = r.*cos(theta); y = r.*sin(theta); % DISPLAY BOX-MULLER SAMPLES figure % X SAMPLES subplot(121); hist(x,100); colormap hot;axis square title(sprintf('Box-Muller Samples Y\n Mean = %1.2f\n Variance = %1.2f\n Kurtosis = %1.2f',mean(x),var(x),3-kurtosis(x))) xlim([-6 6]) % Y SAMPLES subplot(122); hist(y,100); colormap hot;axis square title(sprintf('Box-Muller Samples X\n Mean = %1.2f\n Variance = %1.2f\n Kurtosis = %1.2f',mean(y),var(y),3-kurtosis(y))) xlim([-6 6])
The output of the MATLAB code is shown above. Notice the first, second, and fourth central moments (mean, variance, and kurtosis) of the generated samples are consistent with the standard normal. The Box-Muller transform is another example of of how uniform variables on the interval (0,1) and can be transformed in order to sample from a more complicated distribution.