Archive for July, 2007

The Box-Muller Technique

Friday, July 20th, 2007

So you want a random number from the Normal (Gaussian) distribution, but all you have is a simple random number generator which gives you values from [0.0, 1.0). This came up for me when implementing a Monte Carlo simulation of electron emission from a hot metal surface into vacuum.

The brute-force method for this would be something like this:

  1. Calculate the value of the cumulative distribution function (CDF), i.e. the integral of the Normal distribution, at many points and normalize the whole CDF so that its values fit into the range [0.0, 1.0].
  2. Pick a random number in the range [0.0, 1.0]
  3. Map that random number onto the CDF and find the abscissa where that point is located. This is the random number you want.

This is in fact the general form of the solution to the problem of getting a random number out of any known continuous distribution, and is especially useful when the distribution is not known analytically. But I knew ahead of time that calculating the CDF would be a real bear so I searched around for something easier.

It turns out this problem was attacked in the 1950s in a nice way by Box and Muller. They came up with a nice algorithm for getting random numbers from the Normal distribution which is known as the Box-Muller technique. Several references, including Wikipedia, list more or less the same code for implementing the technique. I found this discussion to be useful. The following code listing gives the “polar form” which returns a random number from the Normal distribution with mean m and standard deviation s. Calls to ranf() yield random floats in [0.0, 1.0].

float box_muller(float m, float s)
float x1, x2, w, y1;
static float y2;
static int use_last = 0;

if (use_last) /* use value from previous call */
y1 = y2;
use_last = 0;
do {
x1 = 2.0 * ranf() - 1.0;
x2 = 2.0 * ranf() - 1.0;
w = x1 * x1 + x2 * x2;
} while ( w >= 1.0 );

w = sqrt( (-2.0 * log( w ) ) / w );
y1 = x1 * w;
y2 = x2 * w;
use_last = 1;

return( m + y1 * s );