Ziggurat algorithm

Sometimes (2.5% of the time, in the case of a normal or exponential distribution when using typical table sizes)[citation needed] more computations are required.

Nevertheless, the algorithm is computationally much faster[citation needed] than the two most commonly used methods of generating normally distributed random numbers, the Marsaglia polar method and the Box–Muller transform, which require at least one logarithm and one square root calculation for each pair of generated values.

The term ziggurat algorithm dates from Marsaglia's paper with Wai Wan Tsang in 2000; it is so named because it is conceptually based on covering the probability distribution with rectangular segments stacked in decreasing order of size, resulting in a figure that resembles a ziggurat.

To use a precomputed table of size n (n = 256 is typical), one chooses x1 such that xn = 0, meaning that the top box, layer n − 1, reaches the distribution's peak at (0, f(0)) exactly.

Ignoring for a moment the problem of layer 0, and given uniform random variables U0 and U1 ∈ [0,1), the ziggurat algorithm can be described as: Step 1 amounts to choosing a low-resolution y coordinate.

For a normal distribution, Marsaglia suggests a compact algorithm: Since x1 ≈ 3.5 for typical table sizes, the test in step 3 is almost always successful.

The algorithm can be performed efficiently with precomputed tables of xi and yi = f(xi), but there are some modifications to make it even faster: It is possible to store the entire table precomputed, or just include the values n, y1, A, and an implementation of f −1(y) in the source code, and compute the remaining values when initializing the random number generator.

When actually filling in the table values, just assume that xn = 0 and yn = f(0), and accept the slight difference in layer n − 1's area as rounding error.

Second, the exact area of the odd-shaped regions is used; they are not rounded up to include the entire rectangle to (xi −1, yi).

If the value sought lies in any of the odd-shaped regions, the alias method is used to choose one, based on its true area.

Third, the almost-triangular shape of most odd-shaped portions is exploited, although this must be divided into three cases depending on the second derivative of the probability distribution function in the selected layer.

Two unit uniform deviates U1 and U2 are chosen, and before they are scaled to the rectangle enclosing the odd-shaped region, their sum is tested.

Only for points very close to the diagonal is it necessary to compute the distribution function f(x) to perform an exact rejection test.

In the one layer which straddles |x| = 1, the normal distribution has an inflection point, and the exact rejection test must be applied if 1−ε 

The tail is handled as in the original Ziggurat algorithm, and can be thought of as a fourth case for the shape of the odd-shaped region to the right.

The Ziggurat algorithm used to generate sample values with a normal distribution . (Only positive values are shown for simplicity.) The pink dots are initially uniform-distributed random numbers. The desired distribution function is first segmented into equal areas "A". One layer i is selected at random by the uniform source at the left. Then a random value from the top source is multiplied by the width of the chosen layer, and the result is x tested to see which region of the layer it falls into with 3 possible outcomes: 1) (left, solid black region) the sample is clearly under the curve and may be output immediately, 2) (right, vertically striped region) the sample value may lie under the curve, and must be tested further. In that case, a random y value within the chosen layer is generated and compared to f(x) . If less, the point is under the curve and the value x is output. If not, (the third case), the chosen point x is rejected and the algorithm is restarted from the beginning.