Marsaglia polar method

The underlying theory may be summarized as follows:

If u is uniformly distributed in the interval

0 ≤ u < 1, then the point

(cos(2πu), sin(2πu))

is uniformly distributed on the unit circumference

x² + y² = 1, and multiplying that point by an independent

random variable ρ whose distribution is

:

will produce a point

:$\backslash left(\backslash rho\backslash cos(2\backslash pi\; u),\backslash rho\backslash sin(2\backslash pi\; u)\backslash right)$

whose coordinates are jointly distributed as two independent standard

normal random variables.

This idea dates back to Laplace, whom Gauss credits with finding the above

:$I=\backslash int\_\{-\backslash infty\}^\backslash infty\; e^\{-x^2/2\}\backslash ,dx$

by taking the square root of

:$I^2\; =\; \backslash int\_\{-\backslash infty\}^\backslash infty\backslash int\_\{-\backslash infty\}^\backslash infty\; e^\{-(x^2+y^2)/2\}\backslash ,dx\backslash ,dy$

=\int_0^{2\pi}\int_0^\infty re^{-r^2/2} \, dr \, d\theta.

The transformation to polar coordinates makes evident that θ is

uniformly distributed (constant density) from 0 to 2π, and that the

radial distance r has density

:$re^\{-r^2/2\}.\; \backslash ,$

(r² has the appropriate chi square distribution.)

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,

A direct application of this idea,

:$x=\backslash sqrt\{-2\backslash ln(u\_1)\}\backslash cos(2\backslash pi\; u\_2),\backslash quad\; y=\backslash sqrt\{-2\backslash ln(u\_1)\}\backslash sin(2\backslash pi\; u\_2)$

is called the Box–Muller transform, in which the chi variate is usually

generated as

:$\backslash sqrt\{-2\backslash ln(u\_1)\},$

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.{{cite web|last=Kanter|first=David|title=Intel's Ivy Bridge Graphics Architecture|url=http://www.realworldtech.com/ivy-bridge-gpu/5/|work=Real World Tech|accessdate=8 April 2013}} Notably for Intel-based machines, one can use fsincos assembler instruction or the expi instruction (available e.g. in D), to calculate complex

: $\backslash operatorname\{expi\}(z)\; =\; e^\{i\; z\}\; =\; \backslash cos(z)\; +\; i\; \backslash sin(z),\; \backslash ,$

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 $r\; =\; \backslash sqrt\{-2\; \backslash ln(u\_1)\}$ and $z\; =\; 2\; \backslash pi\; u\_2.$ Then

: $\backslash \; re^\{i\; z\}\; =\; \backslash sqrt\{-2\; \backslash ln(u\_1)\}\; e^\{i\; 2\; \backslash pi\; u\_2\}\; =\backslash sqrt\{-2\; \backslash ln(u\_1)\}\backslash left[\; \backslash cos(2\; \backslash pi\; u\_2)\; +\; i\; \backslash sin(2\; \backslash pi\; u\_2)\backslash right].$

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 $\backslash sqrt\{x^2\; +\; y^2\}$ radius. In particular, a random point (x, y) inside the unit circle is projected onto the unit circumference by setting $s=x^2+y^2$ and forming the point

:$\backslash left(\; \backslash frac\{x\}\{\backslash sqrt\{s\}\},\; \backslash frac\{y\}\{\backslash sqrt\{s\}\}\; \backslash right),\; \backslash ,$

which is a faster procedure than calculating the cosine and sine. Some researchers argue that the conditional if instruction (for rejecting a point outside of the unit circle), can make programs slower on modern processors equipped with pipelining and branch prediction.This effect can be heightened in a GPU generating many variates in parallel, where a rejection on one processor can slow down many other processors. See section 7 of {{citation

| last1 = Thomas | first1 = David B.

| last2 = Howes | first2 = Lee W.

| last3 = Luk | first3 = Wayne

| editor1-last = Chow | editor1-first = Paul

| editor2-last = Cheung | editor2-first = Peter Y. K.

| contribution = A comparison of CPUs, GPUs, FPGAs, and massively parallel processor arrays for random number generation

| doi = 10.1145/1508128.1508139

| pages = 63–72

| publisher = Association for Computing Machinery

| title = Proceedings of the ACM/SIGDA 17th International Symposium on Field Programmable Gate Arrays, FPGA 2009, Monterey, California, USA, February 22–24, 2009

| year = 2009| isbn = 9781605584102

| citeseerx = 10.1.1.149.6066

| s2cid = 465785

}}. Also this procedure requires about 27% more evaluations of the underlying random number generator (only $\backslash pi/4\; \backslash approx\; 79\backslash \%$ of generated points lie inside of unit circle).

That random point on the circumference is then radially projected the required random distance by means of

:$\backslash sqrt\{-2\backslash ln(s)\},\; \backslash ,$

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.

A simple implementation in Python:

import math

import random

def marsaglia_sample():

while True:

u1 = random.uniform(-1, 1)

u2 = random.uniform(-1, 1)

if (w := u1**2 + u2**2) < 1:

break

z1 = u1 * math.sqrt(-2 * math.log(w) / w)

z2 = u2 * math.sqrt(-2 * math.log(w) / w)

return z1, z2

Simple implementation in Java using the mean and standard deviation:

private static double spare;

private static boolean hasSpare = false;

public static synchronized double generateGaussian(double mean, double stdDev) {

if (hasSpare) {

hasSpare = false;

return spare * stdDev + mean;

} else {

double u, v, s;

do {

u = Math.random() * 2 - 1;

v = Math.random() * 2 - 1;

s = u * u + v * v;

} while (s >= 1 || s == 0);

s = Math.sqrt(-2.0 * Math.log(s) / s);

spare = v * s;

hasSpare = true;

return mean + stdDev * u * s;

}

}

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

double generateGaussian(double mean, double stdDev) {

static double spare;

static bool hasSpare = false;

if (hasSpare) {

hasSpare = false;

return spare * stdDev + mean;

} else {

double u, v, s;

do {

u = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;

v = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;

s = u * u + v * v;

} while (s >= 1.0 || s == 0.0);

s = sqrt(-2.0 * log(s) / s);

spare = v * s;

hasSpare = true;

return mean + stdDev * u * s;

}

}

C++11 GNU GCC libstdc++'s implementation of std::normal_distribution [https://stackoverflow.com/a/25684539 uses] the Marsaglia polar method, as quoted from [https://gcc.gnu.org/onlinedocs/gcc-4.6.0/libstdc++/api/a01001_source.html#l01649 herein].

A simple Julia implementation:

"""

marsagliasample(N)

Generate `2N` samples from the standard normal distribution using the Marsaglia method.

"""

function marsagliasample(N)

z = Array{Float64}(undef,N,2);

for i in axes(z,1)

s = Inf;

while s > 1

z[i,:] .= 2rand(2) .- 1;

s = sum(abs2.(z[i,:]))

end

z[i,:] .*= sqrt(-2log(s)/s);

end

vec(z)

end

"""

marsagliasample(n,μ,σ)

Generate `n` samples from the normal distribution with mean `μ` and standard deviation `σ` using the Marsaglia method.

"""

function marsagliasample(n,μ,σ)

μ .+ σ*marsagliasample(cld(n,2))[1:n];

end

The for loop can be parallelized by using the Threads.@threads macro.

{{Reflist}}

Category:Monte Carlo methods

Category:Pseudorandom number generators

Category:Non-uniform random numbers