Proof of Marsaglia polar method

Question

I studied Polar method and I can use it very well to simulate to Standard Normal Variable.

But I can't figure it out that how it works!

So is there any proof/theorem to learn reasoning behind Polar method?

Thanks! Regards,

Answer 1

First I'll quickly go through the Box-Muller method and then show how the Polar method takes advantage of it.

Suppose that $$ \left( \begin{array}{c} X\\ Y\\ \end{array} \right) \sim N(0, I_2). $$ Note that by multiplying by $\Sigma^{1/2}$ and adding $\mu$ we can get any normal distribution at all, so the important part is getting iid normals.

Define $R^2 = X^2 + Y^2$ and $\Theta = \arctan(Y/X)$. It can be shown (via Jacobian theorem) that $R^2 \sim \textrm{Exp}(2)$, $\Theta \sim \textrm{Unif}(0, 2\pi)$, and $R^2 \perp \Theta$. Using polar coordinates we have that

$$ \left( \begin{array}{c} X\\ Y\\ \end{array} \right) = \left( \begin{array}{c} R \cos\Theta \\ R \sin \Theta\\ \end{array} \right) $$

so if we can sample $R^2$ and $\Theta$ we're good to go. This turns out to be easy: Let $V_1, V_2 \sim \ \textrm{iid} \ \textrm{Unif}(0,1)$. Then by the inverse CDF method $-2\log V_1 := R^2 \sim \textrm{Exp}(2)$ and $2\pi V_2 := \Theta \sim \textrm{Unif}(0, 2\pi)$ with both independent. Plugging these in to the above equation gives us our $X$ and $Y$.

Now I'll show how the Polar method uses this.

Let $\mathcal U = \{(u, v) : u^2 + v^2 \leq 1\}$, i.e. $\mathcal U$ is the unit disk. Take $(X, Y) \sim \textrm{Unif}(\mathcal U)$. We can get this by taking $U_1, U_2 \sim \ \textrm{iid} \ \textrm{Unif}[-1, 1]$ and rejecting until we're in $\mathcal U$.

Now define $T = X^2 + Y^2$, i.e. $T$ is the squared distance from the origin to our point $(X, Y)$.

$$ P(T \leq t) = P(X^2 + Y^2 \leq t) = \frac{\pi (\sqrt t)^2}{\pi (1)^2} = t, $$

i.e. the probability that $T \leq t$ is the ratio of the area of the disk with radius $\sqrt t$ over the area of the entire disk, which is just $t$. This means that $T \sim \textrm{Unif}(0, 1)$.

From the Box-Muller bit we know that

$$ \left( \begin{array}{c} R \cos\Theta \\ R \sin \Theta\\ \end{array} \right) \sim N(0, I_2) $$

when $R^2 \sim \textrm{Exp}(2) \perp \Theta \sim \textrm{Unif}(0, 2\pi)$. $T \sim \textrm{Unif}(0, 1)$ so we can take $R^2 := -2 \log T$.

Now note that for a particular $(X, Y)$ we'll have a triangle like so:

This means that $\sin \Theta = \frac{Y}{\sqrt T}$ and $\cos \Theta = \frac{X}{\sqrt T}$. We can just plug these into the Box-Muller equations and there we go. It also turns out that the independence assumption is satisfied. This frees us from computing any trig functions which can be time consuming.

Answer 2

I think there are 2 points to complete above proof.

1. $\Theta$ of the point uniformly generated in the unit circle is uniform and $\Theta \sim \mathrm{Unif}[0,2\pi]$.

2. $T$ is independent of $\Theta$.

This is not a hard work. For the first part, $$ \int^{2\pi}_{0}\int^1_0 r f_{(X,Y)}(r \cos\theta,r \sin\theta)~\mathrm dr\mathrm d\theta = 1 \\ \Longrightarrow f_{(X,Y)}(r \cos\theta,r \sin\theta) = \frac{1}{\pi} $$ since $X$ and $Y$ generated uniformly(constant joint density). Then we have $$ f_{(\Theta,R)}(\theta,r) = \frac{r}{\pi} $$ by transformation of pdf.

$$f_{\Theta}(\theta) = \int^{1}_{0}\frac{r}{\pi}~\mathrm dr = \frac{1}{2 \pi} $$ Thus, $\Theta \sim \textrm{Unif}[0,2\pi]$. For the second part, $$ f_{R}(r) = \int^{2\pi}_{0}\frac{r}{\pi}~\mathrm dr = 2r. $$ Thus, r is independent of $\Theta$ also $T = r^2$ is independent of $\Theta$.