Polar Method

We definitely don’t have a shortage of methods for generating normally-distributed random numbers from a source of uniformly-distributed random numbers. One of such methods is the so called Polar Method, a variation of the Box-Muller Transform, which I already described before. You might want to take a look at it before reading this.

The Polar Method was initially proposed, it seems, by the great George Marsaglia along with T. A. Bray, in the 1964 paper A Convenient Method for Generating Normal Variables. This paper must be really good, because it is available online for the whopping price of US$35.00, or seven bucks per page! Sorry, but no thanks.

Anyway, let’s see how it works. What we’ll do here is to construct a function that uses the Polar Method to generate normally-distributed random numbers from a source of uniformly-distributed random numbers (namely, Math.random()). As in the case of the Box-Muller Transform, our RNG (Random Number Generator) will generate two random numbers at once, but we want to return only one at a time. Therefore, here we’ll also use two global variables to let us keep that spare random number around until we need it:

⟨Polar Method normal RNG globals⟩ =

let polarMethodNormalRNG_haveSpare = false;
let polarMethodNormalRNG_spare = 0.0;

Our RNG has the structure shown below: it makes use of these global variables and looks like a regular RNG: takes no parameters and returns a number.

⟨Polar Method normal RNG⟩ =

⟨Polar Method normal RNG globals⟩

function polarMethodNormalRNG(): number {
    ⟨Polar Method normal RNG body⟩
}

On to the implementation of polarMethodNormalRNG(): we first check if we have that spare number available. If we do, we return it immediately.

⟨Polar Method normal RNG body⟩ =

if (polarMethodNormalRNG_haveSpare) {
    polarMethodNormalRNG_haveSpare = false;
    return polarMethodNormalRNG_spare;
}

Otherwise, we reached the interesting part of this explanation: the point where we’ll actually use the Polar Method!

We first create two random numbers, u1 and u2, taken from a uniform distribution, so that the point (u1, u2) lies in the unit circle. Notice that we do so by trial and error, rejecting points that fall outside it and trying again:

⟨Polar Method normal RNG body⟩ +=

let u0, u1, s: number;

do {
    u0 = Math.random() * 2.0 - 1.0;
    u1 = Math.random() * 2.0 - 1.0;
    s = u0*u0 + u1*u1;
} while (s >= 1 || s == 0);

Then we scale this point by a magic value m¹ that will leave the resulting point (z0, z1) in exactly the same position as the tip of the vector we generated in the Box-Muller Transform:

⟨Polar Method normal RNG body⟩ +=

let m = Math.sqrt(-2.0 * Math.log(s) / s);
let z0 = u0 * m;
let z1 = u1 * m;

With the Box-Muller Transform, the generated vector was in polar coordinates and therefore we had to use sine and cosine to project it into the Cartesian axes and thus obtain the desired normally-distributed numbers. But here things are different: the point (z0, z1) is already born in rectangular coordinates, so we can just use it directly²: z0 and z1 are the numbers we want. So we return z0 while keeping z1 at hand for future invocations:

⟨Polar Method normal RNG body⟩ +=

polarMethodNormalRNG_haveSpare = true;
polarMethodNormalRNG_spare = z1;

return z0;

Again, my limited mathematical knowledge is not enough to explain all the details about how the Polar Method works.³ But, really, the thing about the Polar Method is not mathematical (the math here is pretty much the same as in the Box-Muller transform). The important idea here is computational: even though the Polar Method discards and re-generates more than 20% of all generated points⁴ in that while loop, it will still be faster than Box-Muller because it avoids calling sin() and cos(). (Ahem, usually it will be faster. Don’t miss the performance considerations below.)

Visualization

This visualization shows the Polar Method doing its work. Depending on the button you click, it will generate one, ten, or one hundred pairs of random numbers. I suggest you to first generate one point at a time to get a feeling of what the algorithm does in a micro-level; then, generate lots of points at once to see how it behaves statistically. There are some notes after the diagram that may help you understand it.

What’s happening here?

• Each click in “One” shows the steps polarMethodNormalRNG() goes through to generate the random numbers:
  • First it shows the point (u0, u1). If the point lies outside the unit circle, the point is drawn in red and the visualization ends. (The actual algorithm, as we’ve seen, would continue trying until getting a point that lies inside the unit circle.)
  • If the point is inside the unit circle, it is drawn in blue and translated to its final position: (z0, z1). This movement depicts the multiplication by m. z0 and z1 are the two random numbers generated by polarMethodNormalRNG(), so the Polar Method is complete at this point.
• For clarity sake, we show one additional animated step: the point (z0, z1) being projected into the vertical and horizontal axes, and the creation of histograms showing where these projections fall (in other words, the histograms show the distribution of values generated by the Polar Method).
• 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 follow a normal distribution.
• And that’s exactly what the Polar Method 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 Polar Method 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.

Performance Considerations

Descriptions of the Polar Method often praise it for its improved performance over the Box-Muller Transform. After all, given that the Polar Method avoids calling those expensive trigonometric functions it must be faster, right?

Well, I guess this was close an absolute truth when the Polar Method was introduced in the 1960s, but nowadays in addition to the CPUs we all love, we also have things like GPUs and FPGAs. In a nice 2009 paper called (hold you breath) A Comparison of CPUs, GPUs, FPGAs, and Massively Parallel Processor Arrays for Random Number Generation, the authors investigate different algorithms related with random number generation and how they perform in different categories of hardware.

An interesting conclusion of the paper is that for each of the four hardware categories analyzed, a different algorithm was deemed the best. For GPUs, in particular, the Box-Muller Transform was the top choice, despite all that trigonometry. This happened for two reasons. First, GPUs can compute sines and cosines really fast. Second, to make a long story short, GPUs don’t deal well with branching and thus that while loop in the Polar Method was a performance disadvantage.

So, I guess it is fair to say that the Polar Method is faster than the Box-Muller transform on most of the hardware we use most of the time⁵ – but keep your mind open.

Nutrition Facts

AKA: Marsaglia Polar Method.

See also: This method is a variation of the Box-Muller Transform.

Keywords: Normal distribution, Gaussian distribution.

Source code: polar_method.ts (the code generated from the snippets on this page); dump_polar_method.ts (the interactive diagram above); dump.ts (helper utilities).

References

• David B. Thomas, Lee Howes, Wayne Luk. A Comparison of CPUs, GPUs, FPGAs, and Massively Parallel Processor Arrays for Random Number Generation. FPGA ‘09: Proceedings of the ACM/SIGDA International Symposium on Field Programmable Gate Arrays, 2009.
• Wikipedia, Marsaglia polar method.

 Share!

• Share on Facebook
   
• Tweet
   
• Post to Tumblr
   
• Pin it
   
• Submit to Reddit