Box–Muller transform

Remarks

Let $R = \sqrt{-2 \ln U_1}$ and $\Theta = 2\pi U_2$, so that $(Z_1, Z_2)$ has polar coordinates $(R, \Theta)$. For $t \ge 0$, $P(R^2 > t) = P(U_1 < e^{-t/2}) = e^{-t/2}$, so $R^2$ has the exponential distribution with mean $2$, and $R$ has density $r e^{-r^2/2}$ on $(0, \infty)$. Also $\Theta$ is uniform on $(0, 2\pi)$ and independent of $R$. By the change of variables formula with Jacobian $r$, the pair $(Z_1, Z_2)$ has density $\frac{1}{2\pi} e^{-r^2/2} = \frac{1}{\sqrt{2\pi}} e^{-z_1^2/2} \cdot \frac{1}{\sqrt{2\pi}} e^{-z_2^2/2}$, which is the product of two standard normal densities. Example: $U_1 = e^{-2}$ and $U_2 = 1/8$ give $R = 2$, $\Theta = \pi/4$, and $(Z_1, Z_2) = (\sqrt{2}, \sqrt{2})$. Independence of $U_1$ and $U_2$ is necessary: with $U_2 = U_1$ the point $(Z_1, Z_2)$ lies on one curve in the plane, so the pair is not a standard normal pair. To obtain a normal variable with mean $\mu$ and standard deviation $s$, use $\mu + s Z_1$.

Going deeper

For a student

One cannot invert the normal distribution function in closed form, so plain inverse transform sampling does not give normal variables cheaply; going to two dimensions and polar coordinates solves this. Marsaglia's polar method avoids the sine and cosine by sampling a uniform point in the unit disc by rejection and rescaling it. In floating point, $U_1$ must avoid $0$, and the tail of $Z_1$ is truncated by the smallest representable $U_1$. For speed, many libraries now prefer the ziggurat algorithm, which uses mostly one multiplication and one comparison per sample.