Marsaglia polar method


The Marsaglia polar method is a pseudo-random number sampling method for generating a pair of independent standard normal random variables. While it is superior to the Box–Muller transform, the Ziggurat algorithm is even more efficient.
Standard normal random variables are frequently used in computer science, computational statistics, and in particular, in applications of the Monte Carlo method.
The polar method works by choosing random points in the square −1 < x < 1, −1 < y < 1 until
and then returning the required pair of normal random variables as
or, equivalently,
where and represent the cosine and sine of the angle that the vector makes with x axis.

Theoretical basis

The underlying theory may be summarized as follows:
If u is uniformly distributed in the interval
0 ≤ u < 1, then the point
, sin)
is uniformly distributed on the unit circumference
x2 + y2 = 1, and multiplying that point by an independent
random variable ρ whose distribution is
will produce a point
whose coordinates are jointly distributed as two independent standard
normal random variables.

History

This idea dates back to Laplace, whom Gauss credits with finding the above
by taking the square root of
The transformation to polar coordinates makes evident that θ is
uniformly distributed from 0 to 2π, and that the
radial distance r has density
This method of producing a pair of independent standard normal variates by radially projecting a random point on the unit circumference to a distance given by the square root of a chi-square-2 variate is called the polar method for generating a pair of normal random variables,

Practical considerations

A direct application of this idea,
is called the Box–Muller transform, in which the chi variate is usually
generated as
but that transform requires logarithm, square root, sine and cosine functions. On some processors, the cosine and sine of the same argument can be calculated in parallel using a single instruction. Notably for Intel-based machines, one can use fsincos assembler instruction or the expi instruction, to calculate complex
and just separate the real and imaginary parts.
Note:
To explicitly calculate the complex-polar form use the following substitutions in the general form,
Let and Then
In contrast, the polar method here removes the need to calculate a cosine and sine. Instead, by solving for a point on the unit circle, these two functions can be replaced with the x and y coordinates normalized to the radius. In particular, a random point inside the unit circle is projected onto the unit circumference by setting and forming the point
which is a faster procedure than calculating the cosine and sine. Some researchers argue that the conditional if instruction, can make programs slower on modern processors equipped with pipelining and branch prediction. Also this procedure requires about 27% more evaluations of the underlying random number generator.
That random point on the circumference is then radially projected the required random distance by means of
using the same s because that s is independent of the random point on the circumference and is itself uniformly distributed from 0 to 1.

Implementation

Simple implementation in Java using the mean and standard deviation:

private static double spare;
private static boolean hasSpare = false;
public static synchronized double generateGaussian

A non-thread safe implementation in C++ using the mean and standard deviation:

double generateGaussian

C++11 GNU GCC libstdc++'s implementation of std::normal_distribution the Marsaglia polar method, as quoted from .