Ziggurat Algorithm - Theory of Operation

Theory of Operation

The ziggurat algorithm is a rejection sampling algorithm; it randomly generates a point in a distribution slightly larger than the desired distribution, then tests whether the generated point is inside the desired distribution. If not, it tries again. Given a random point underneath a probability distribution curve, its x coordinate is a random number with the desired distribution.

The distribution the ziggurat algorithm chooses from is made up of n equal-area regions; n − 1 rectangles that cover the bulk of the desired distribution, on top of a non-rectangular base that includes the tail of the distribution.

Given a monotone decreasing probability distribution function f(x), defined for all x≥0, the base of the ziggurat is defined as all points inside the distribution and below y₁ = f(x₁). This consists of a rectangular region from (0, 0) to (x₁, y₁), and the (typically infinite) tail of the distribution, where x > x₁ (and y < y₁).

This layer (call it layer 0) has area A. On top of this, add a rectangular layer of width x₁ and height A/x₁, so it also has area A. The top of this layer is at height y₂ = y₁ + A/x₁, and intersects the distribution function at a point (x₂, y₂), where y₂ = f(x₂). This layer includes every point in the distribution function between y₁ and y₂, but (unlike the base layer) also includes points such as (x₁, y₂) which are not in the desired distribution.

Further layers are then stacked on top. To use a precomputed table of size n (n = 256 is typical), one chooses x₁ such that x_n=0, meaning that the top box, layer n−1, reaches the distribution's peak at (0, f(0)) exactly.

Layer i extends vertically from y_i to y_i+1, and can be divided into two regions horizontally: the (generally larger) portion from 0 to x_i+1 which is entirely contained within the desired distribution, and the (small) portion from x_i+1 to x_i, which is only partially contained.

Ignoring for a moment the problem of layer 0, and given uniform random variables U₀ and U₁ ∈ [0,1), the ziggurat algorithm can be described as:

1. Choose a random layer 0 ≤ i < n.
2. Let x = U₀x_i.
3. If x < x_i+1, return x.
4. Let y = y_i + U₁(y_i+1−y_i).
5. Compute f(x). If y < f(x), return x.
6. Otherwise, choose new random numbers and go back to step 1.

Step 1 amounts to choosing a low-resolution y coordinate. Step 3 tests if the x coordinate is clearly within the desired distribution function without knowing more about the y coordinate. If it is not, step 4 chooses a high-resolution y coordinate and step 5 does the rejection test.

With closely spaced layers, the algorithm terminates at step 3 a very large fraction of the time. Note that for the top layer n−1, however, this test always fails, because x_n = 0.

Layer 0 can also be divided into a central region and an edge, but the edge is an infinite tail. To use the same algorithm to check if the point is in the central region, generate a fictitious x₀ = A/y₁. This will generate points with x < x₁ with the correct frequency, and in the rare case that layer 0 is selected and x ≥ x₁, use a special fallback algorithm to select a point at random from the tail. Because the fallback algorithm is used less than one time in a thousand, speed is not essential.

Thus, the full ziggurat algorithm for one-sided distributions is:

1. Choose a random layer 0 ≤ i < n.
2. Let x = U₀x_i
3. If x < x_i+1, return x.
4. If i=0, generate a point from the tail using the fallback algorithm.
5. Let y = y_i + U₁(y_i+1−y_i).
6. Compute f(x). If y < f(x), return x.
7. Otherwise, choose new random numbers and go back to step 1.

For a two-sided distribution, of course, the result must be negated 50% of the time. This can often be done conveniently by choosing U₀ ∈ (−1,1) and, in step 3, testing if |x| < x_i+1.