Box-Muller Transform

Given a source of uniformly-distributed random numbers, generate normally-distributed random numbers

Programming
DUMP

This algorithm is named after George Edward Pelham Box and Mervin Edgar Muller, who published it on a two-page paper in 1958. The idea was not original, though: it appeared already in the 1934 book Fourier Transforms in the Complex Domain, by Raymond E. A. C. Paley and Norbert Wiener. Stigler’s law strikes again!

How exactly the Box-Muller transform works (in a deeper mathematical sense) is honestly beyond my understanding, but I hope this will at least provide a good intuition about it. The transform transforms two random numbers taken from a uniform distribution in the (0, 1) interval into two random numbers that follow a normal (Gaussian) distribution. Two numbers in, two numbers out:

⟨Box-Muller Transform⟩ =

function boxMullerTransform(u0: number, u1: number): [number, number] {
    ⟨Box-Muller Transform body⟩
}

For the transformation itself, we start by converting the input values u0 and u1 into a vector of length r and direction theta:

⟨Box-Muller Transform body⟩ =

let r = Math.sqrt(-2 * Math.log(u0));
let theta = 2 * Math.PI * u1;

We then transform this vector from the polar form to the rectangular form:

⟨Box-Muller Transform body⟩ +=

let z0 = r * Math.cos(theta);
let z1 = r * Math.sin(theta);

And both of these two numbers, z0 and z1, are random numbers following a standard normal distribution (that is, one with mean equals to 0 and standard deviation equals to 1). That’s exactly what we wanted, so we return these numbers and call it a day:

⟨Box-Muller Transform body⟩ +=

return [z0, z1];

Visualizing the Box-Muller Transform

Before anything else, credit where credit is due: this visualization is based on the very nice dynamic SVG visualization created by Wikipedia user Cmglee. I created my visualization from scratch, but the influence of Cmglee’s work is evident.

So, here’s the visualization. It generates pairs of random numbers (one, ten, or a hundred pairs at a time, depending on the button you click), and shows the Box-Muller transform doing its work. If you need some help to understand what is going on here, just keep reading after the diagram.

What’s happening here?

  • Each click in “One” generates a pair of random numbers (taken from a uniform distribution).
  • We use the boxMullerTransform() function to transform this pair of numbers. Recall from the explanation above that the transformation works in stages:
    • First, it generates a vector from u0 and u1; you can see this vector in the visualization.
    • Second, it converts the vector into its rectangular form, by projecting it into the horizontal and vertical axes; we show these projections as those cute falling circles. The point where the circles fall are the outputs of the transform (z0 and z1).
  • In this visualization, we also show histograms showing where the circles fall (in other words, the histograms show the distribution of values generated by the Box-Muller transform).
  • After you generate enough pairs of numbers (say, some hundreds of pairs), both histograms should resemble a normal distribution.
  • This shows that z0 and z1 (the outputs of boxMullerTransform()) follow a normal distribution.
  • And that’s exactly what the Box-Muller transform does: it transforms random numbers in a uniform distribution into random numbers in a normal distribution.
  • You may be thinking: a normal distribution is unbounded, so the Box-Muller transform could produce numbers from -∞ to +∞. How come in this visualization all the output values fit neatly in the [-3, +3] range? Simple: because my visualization discards any value that falls outside of this range. This is just for the sake of a clearer visualization. I am not really cheating.

Generating normally-distributed random numbers

What we saw so far was an implementation of the Box-Muller transform per se. But most people using it just want to generate normally-distributed random numbers. So, if you search the web for “Box-Muller transform”, perhaps most of the implementations you’ll find will include some extra machinery to make it behave more like a random number generator (RNG) than like the plain transform.

Let’s see one implementation like this. First, notice that a typical RNG generates one random number at a time, but our transform generates two, z0 and z1 at once! We don’t want to discard this extra number, so we’ll keep it around until we need it. Thus, here we define two global variables: one to tell if we have this extra number available, and the other to store the number itself:

⟨Box-Muller normal RNG globals⟩ =

let boxMullerNormalRNG_haveSpare = false;
let boxMullerNormalRNG_spare = 0.0;

Our RNG based on the Box-Muller transform will have the structure shown below: it makes use of these global variables and looks like a regular RNG: takes no parameters and returns a number.

⟨Box-Muller normal RNG⟩ =

⟨Box-Muller normal RNG globals⟩

function boxMullerNormalRNG(): number {
    ⟨Box-Muller normal RNG body⟩
}

In the implementation of boxMullerNormalRNG(), we first check if we have that spare number available. If we do, we return it immediately.

⟨Box-Muller normal RNG body⟩ =

if (boxMullerNormalRNG_haveSpare) {
    boxMullerNormalRNG_haveSpare = false;
    return boxMullerNormalRNG_spare;
}

If the extra random number is not available, though luck: we’ll need to do work to generate the random numbers. First, we generate our uniformly-distributed random numbers in the (0,1) range.

⟨Box-Muller normal RNG body⟩ +=

let u0: number;
do {
    u0 = Math.random();
} while (u0 <= 0.0 || u0 >= 1.0);

let u1: number;
do {
    u1 = Math.random();
} while (u1 <= 0.0 || u1 >= 1.0);

Now we are ready to use the Box-Muller transform itself to transform these uniformly-distributed numbers into normally-distributed ones.

⟨Box-Muller normal RNG body⟩ +=

let [z0, z1] = boxMullerTransform(u0, u1);

We got two numbers, but need only one for now. So, keep the spare one, z1 in our ugly-named globals for the future.

⟨Box-Muller normal RNG body⟩ +=

boxMullerNormalRNG_haveSpare = true;
boxMullerNormalRNG_spare = z1;

And the other normally-distributed number, z0, is the one we return:

⟨Box-Muller normal RNG body⟩ +=

return z0;

That’s it. An RNG that generates numbers following a standard normal distribution, based on the Box-Muller transform!

Nutrition Facts

AKA: Box-Muller Transformation.

See also: The Polar Method is a variation of this transform that doesn’t call the potentially expensive sine and cosine functions.

Keywords: Normal distribution, Gaussian distribution.

Source code: box-muller_transform.ts (the code generated from the snippets on this page); dump_box-muller_transform.ts (the interactive diagram above); dump.ts (helper utilities).

References

Changelog

  • 2020-03-17: Renamed variables u1, u2 to u0, u1 (for consistency).
  • 2020-01-28: Initial publication.

 Share!

 
Comments powered by Disqus