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(u1: number, u2: number): [number, number] {
⟨Box-Muller Transform body⟩
}
```

For the transformation itself, we start by converting the input values `u1`

and
`u2`

into a vector of length `r`

and direction `theta`

:

⟨Box-Muller Transform body⟩ =

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

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`

).

- First, it generates a vector from
- 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.

**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 the interactive diagram above); `dump.ts`

(helper utilities).

## References

- Wikipedia,
*Box-Muller Transform*. - Eric W. Weisstein,
*Box-Muller Transformation*. Wolfram's MathWorld. - George Edward Pelham Box and Mervin Edgar Muller,
*A Note on the Generation of Random Normal Deviates*. The Annals of Mathematical Statistics, 1958.