Ziggurat algorithm

Remarks

Each layer has area $v$ by clauses 1 and 2, and the layers overlap only in boundary segments of measure zero. So a uniform choice of layer followed by a uniform point in it gives a point uniform on the union $Z = L_1 \cup \cdots \cup L_n$. Since $f$ is non-increasing, $Z$ contains the region $\set{(x, y) : x \ge 0,\ 0 \le y \le f(x)}$ under the graph. A uniform point of $Z$ conditioned to lie under the graph is uniform there, and its first coordinate has density proportional to $f$. This is rejection sampling with a uniform envelope on $Z$. The speed comes from a shortcut: in layer $i \le n-1$, if $X < x_{i-1}$ then $Y \le f(x_{i-1}) \le f(X)$, so the point is accepted without evaluating $f$. For large $n$, such as $n = 128$ or $256$, almost all trials end this way, after one random integer, one multiplication and one comparison. A uniform point of the base $L_n$ is drawn as follows. With probability $r f(r)/v$ it is a uniform point of the rectangle $[0, r] \times [0, f(r)]$. Otherwise $X$ is drawn from the tail density proportional to $f$ on $(r, \infty)$, which needs a separate tail sampler, and $Y$ is uniform on $[0, f(X)]$. For the standard normal distribution, one samples $\abs{Z}$ with $f(x) = e^{-x^2/2}$ and attaches a random sign. The numbers $x_i$, $r$ and $v$ are computed once, by solving the equations of clauses 1 and 2 numerically.

Going deeper

For a student

The algorithm, due to Marsaglia and Tsang, is a form of rejection sampling with a piecewise constant envelope that fits the density tightly. It is used for normal and exponential random numbers in many numerical libraries. For symmetric densities such as the normal, it samples the right half and attaches a random sign; often the sign comes from a spare bit of the same random integer. Its exactness depends on good uniform variates and an exact tail sampler; the tail of the normal is usually handled by a method of Marsaglia based on exponential variates.