Ziggurat algorithm

<h2 id="theory-of-operation">Theory of operation</h2>
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 density 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 density function f(x), defined for all x ≥ 0, the base of the ziggurat is defined as all points inside the distribution and below y1 = f(x1). This consists of a rectangular region from (0, 0) to (x1, y1), and the (typically infinite) tail of the distribution, where x > x1 (and y < y1).
This layer (call it layer 0) has area A. On top of this, add a rectangular layer of width x1 and height A/x1, so it also has area A. The top of this layer is at height y2 = y1 + A/x1, and intersects the density function at a point (x2, y2), where y2 = f(x2). This layer includes every point in the density function between y1 and y2, but (unlike the base layer) also includes points such as (x1, y2) 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 x1 such that xn = 0, meaning that the top box, layer n − 1, reaches the distribution's peak at (0, f(0)) exactly.
Layer i extends vertically from yi to yi +1, and can be divided into two regions horizontally: the (generally larger) portion from 0 to xi +1 which is entirely contained within the desired distribution, and the (small) portion from xi +1 to xi, which is only partially contained.
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:

<ol><li>Choose a random layer 0 ≤ i < n.</li>
<li>Let x = U0xi.</li>
<li>If x < xi +1, return x.</li>
<li>Let y = yi + U1(yi +1 − yi).</li>
<li>Compute f(x). If y < f(x), return x.</li>
<li>Otherwise, choose new random numbers and go back to step 1.</li></ol>
Step 1 amounts to choosing a low-resolution y coordinate. Step 3 tests if the x coordinate is clearly within the desired density 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. For the top layer n − 1, however, this test always fails, because xn = 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 x0 = A/y1. This will generate points with x < x1 with the correct frequency, and in the rare case that layer 0 is selected and x ≥ x1, 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:

<ol><li>Choose a random layer 0 ≤ i < n.</li>
<li>Let x = U0xi</li>
<li>If x < xi +1, return x.</li>
<li>If i = 0, generate a point from the tail using the fallback algorithm.</li>
<li>Let y = yi + U1(yi +1 − yi).</li>
<li>Compute f(x). If y < f(x), return x.</li>
<li>Otherwise, choose new random numbers and go back to step 1.</li></ol>
For a two-sided distribution, the result must be negated 50% of the time. This can often be done conveniently by choosing U0 ∈ (−1,1) and, in step 3, testing if |x| < xi +1.

<h3>Fallback algorithms for the tail</h3>
Because the ziggurat algorithm only generates most outputs very rapidly, and requires a fallback algorithm whenever x > x1, it is always more complex than a more direct implementation. The specific fallback algorithm depends on the distribution.
For an exponential distribution, the tail looks just like the body of the distribution. One way is to fall back to the most elementary algorithm E = −ln(U1) and let x = x1 − ln(U1). Another is to call the ziggurat algorithm <a href="/facts/Recursion/CxSKxJoW">recursively</a> and add x1 to the result.
For a normal distribution, Marsaglia suggests a compact algorithm:

<ol><li>Let x = −ln(U1)/x1.</li>
<li>Let y = −ln(U2).</li>
<li>If 2y > x2, return x + x1.</li>
<li>Otherwise, go back to step 1.</li></ol>
Since x1 ≈ 3.5 for typical table sizes, the test in step 3 is almost always successful. Since −ln(U1) is an exponentially distributed variate, an implementation of the exponential distribution may be used.

<h3>Optimizations</h3>
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:

<ul><li>Nothing in the ziggurat algorithm depends on the probability distribution function being normalized (integral under the curve equal to 1), removing <a href="/facts/Normalizing_constant/KKTfybqr">normalizing constants</a> can speed up the computation of f(x).</li>
<li>Most uniform random number generators are based on integer random number generators which return an integer in the range [0,  232 − 1]. A table of 2−32xi lets you use such numbers directly for U0.</li>
<li>When computing two-sided distributions using a two-sided U0 as described earlier, the random integer can be interpreted as a signed number in the range [−231, 231 − 1], and a scale factor of 2−31 can be used.</li>
<li>Rather than comparing U0xi to xi +1 in step 3, it is possible to precompute xi +1/xi and compare U0 with that directly. If U0 is an integer random number generator, these limits may be premultiplied by 232 (or 231, as appropriate) so an integer comparison can be used.</li>
<li>With the above two changes, the table of unmodified xi values is no longer needed and may be deleted.</li>
<li>When generating <a href="/facts/IEEE_754/dD3e7zAl">IEEE 754</a> single-precision <a href="/facts/Floating-point_arithmetic/eIckahxe">floating point values</a>, which only have a 24-bit mantissa (including the implicit leading 1), the least-significant bits of a 32-bit integer random number are not used. These bits may be used to select the layer number. (See the references below for a detailed discussion of this.)</li>
<li>The first three steps may be put into an <a href="/facts/Inline_function/foqU4fpQ">inline function</a>, which can call an out-of-line implementation of the less frequently needed steps.</li></ul>
<h3>Generating the tables</h3>
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.
As previously described, you can find xi = f −1(yi) and yi +1 = yi + A/xi. Repeat n − 1 times for the layers of the ziggurat. At the end, you should have yn = f(0). There will be some <a href="/facts/Round-off_error/aF9N8X4I">round-off error</a>, but it is a useful <a href="/facts/Sanity_test/4l0uVCdR">sanity test</a> to see that it is acceptably small.
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.

<h3>Finding x1 and A</h3>
Given an initial (guess at) x1, you need a way to compute the area t of the tail for which x > x1. For the exponential distribution, this is just e−x1, while for the normal distribution, assuming you are using the unnormalized f(x) = e−x2/2, this is √π/2 <a href="/facts/Error_function/Ca2H7pGr">erfc</a>(x/√2). For more awkward distributions, <a href="/facts/Numerical_integration/MwJnvcDV">numerical integration</a> may be required.
With this in hand, from x1, you can find y1 = f(x1), the area t in the tail, and the area of the base layer A = x1y1 + t.
Then compute the series yi and xi as above. If yi > f(0) for any i < n, then the initial estimate x1 was too low, leading to too large an area A. If yn < f(0), then the initial estimate x1 was too high.
Given this, use a <a href="/facts/Root-finding_algorithm/IxTDkhBq">root-finding algorithm</a> (such as the <a href="/facts/Bisection_method/3tWmtoIp">bisection method</a>) to find the value x1 which produces yn−1 as close to f(0) as possible. Alternatively, look for the value which makes the area of the topmost layer, xn−1(f(0) − yn−1), as close to the desired value A as possible. This saves one evaluation of f −1(x) and is actually the condition of greatest interest.

<h2 id="mcfarlands-variation">McFarland's variation</h2>
Christopher D. McFarland has proposed a further-optimized version.<a class="footnote-ref" id="fnref:1" href="#fn:1">1</a> This applies three algorithmic changes, at the expense of slightly larger tables.
First, the common case considers only the rectangular portions, from (0, yi −1) to (xi, yi) The odd-shaped regions to the right of these (mostly almost triangular, plus the tail) are handled separately. This simplifies and speeds up the algorithm's <a href="/facts/Fast_path/xAf8p4Nv">fast path</a>.
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). This increases the probability that the fast path will be used.
One major consequence of this is that the number of layers is slightly less than n. Even though the area of the odd-shaped portions is taken exactly, the total adds up to more than one layer's worth. The area per layer is adjusted so that the number of rectangular layers is an integer. If the initial 0 ≤ i < n exceeds the number of rectangular layers, phase 2 proceeds.
If the value sought lies in any of the odd-shaped regions, the <a href="/facts/Alias_method/XVYcurtk">alias method</a> is used to choose one, based on its true area. This is a small amount of additional work, and requires additional alias tables, but chooses one of the layers' right-hand sides.
The chosen odd-shaped region is rejection sampled, but if a sample is rejected, the algorithm does not return to the beginning. The true area of each odd-shaped region was used to choose a layer, so the rejection-sampling loop stays in that layer until a point is chosen.
Third, the almost-triangular shape of most odd-shaped portions is exploited, although this must be divided into three cases depending on the <a href="/facts/Second_derivative/3aA6j4EY">second derivative</a> of the probability distribution function in the selected layer.
If the function is convex (as the exponential distribution is everywhere, and the normal distribution is for |x| > 1), then the function is strictly contained within the lower triangle. 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. If U1 + U2 > 1, the point is in the upper triangle and can be reflected to (1−U1, 1−U2). Then, if U1 + U2 < 1−ε, for some suitable tolerance ε, the point is definitely below the curve and can immediately be accepted. Only for points very close to the diagonal is it necessary to compute the distribution function f(x) to perform an exact rejection test. (The tolerance ε should in theory depend on the layer, but a single maximum value can be used on all layers with little loss.)
If the function is concave (as the normal distribution is for |x| < 1), it includes a small portion of the upper triangle so reflection is impossible, but points whose normalized coordinates satisfy U1 + U2 ≤ 1 can be immediately accepted, and points for which U1 + U2 > 1+ε can be immediately rejected.
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−ε <U1 + U2 < 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.

<ul><li><a href="/facts/George_Marsaglia/Jlu1Defc">Marsaglia, George</a>; Tsang, Wai Wan (2 October 2000). <a href="http://www.jstatsoft.org/v05/i08/paper">"The Ziggurat Method for Generating Random Variables"</a>. Journal of Statistical Software. 5 (8). Retrieved 2007-06-20. This paper numbers the layers from 1 starting at the top, and makes layer 0 at the bottom a special case, while the explanation above numbers layers from 0 at the bottom.</li>
<li><a href="https://www.jstatsoft.org/article/downloadSuppFile/v005i08/rnorrexp.c">C implementation of the ziggurat method for the normal density function and the exponential density function</a>, that is essentially a copy of the code in the paper. (Potential users should be aware that this C code assumes 32-bit integers.)</li>
<li><a href="http://heliosphan.org/zigguratalgorithm/zigguratalgorithm.html">A C# implementation</a> of the ziggurat algorithm and overview of the method.</li>
<li>Jurgen A. Doornik (2005). <a href="http://www.doornik.com/research/ziggurat.pdf">"An Improved Ziggurat Method to Generate Normal Random Samples"</a> (PDF). Nuffield College, Oxford. Retrieved 2007-06-20. Describes the hazards of using the least-significant bits of the integer random number generator to choose the layer number.</li>
<li><a href="https://au.mathworks.com/company/newsletters/articles/normal-behavior.html">Normal Behavior</a> By Cleve Moler, MathWorks, describing the ziggurat algorithm introduced in <a href="/facts/MATLAB/qPjLISCk">MATLAB</a> version 5, 2001.</li>
<li><a href="https://blogs.mathworks.com/cleve/2015/05/18/the-ziggurat-random-normal-generator/">The Ziggurat Random Normal Generator</a> Blogs of MathWorks, posted by Cleve Moler, May 18, 2015.</li>
<li>David B. Thomas; Philip H.W. Leong; Wayne Luk; John D. Villasenor (October 2007). <a href="http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf">"Gaussian Random Number Generators"</a> (PDF). ACM Computing Surveys. 39 (4): 11:1–38. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1145%2F1287620.1287622">10.1145/1287620.1287622</a>. <a href="/facts/ISSN_(identifier)/DPAflDvU">ISSN</a> <a href="https://search.worldcat.org/issn/0360-0300">0360-0300</a>. <a href="/facts/S2CID_(identifier)/ldJsHa2Y">S2CID</a> <a href="https://api.semanticscholar.org/CorpusID:10948255">10948255</a>. Retrieved 2009-07-27. [W]hen maintaining extremely high statistical quality is the first priority, and subject to that constraint, speed is also desired, the Ziggurat method will often be the most appropriate choice. Comparison of several algorithms for generating <a href="/facts/Gaussian_distribution/UapjjPyQ">Gaussian</a> random numbers.</li>
<li>Nadler, Boaz (2006). "Design Flaws in the Implementation of the Ziggurat and Monty Python methods (And some remarks on Matlab randn)". <a href="/facts/ArXiv_(identifier)/H6EtgnBe">arXiv</a>:<a href="https://arxiv.org/abs/math/0603058">math/0603058</a>.. Illustrates problems with underlying uniform pseudo-random number generators and how those problems affect the ziggurat algorithm's output.</li>
<li>Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian (13–16 July 2009). <a href="http://www.ee.cooper.edu/~nummey/ersa2009.pdf">Hardware-Optimized Ziggurat Algorithm for High-Speed Gaussian Random Number Generators</a> (PDF). 2009 International Conference on Engineering of Reconfigurable Systems & Algorithms. Las Vegas.</li>
<li><a href="/facts/George_Marsaglia/Jlu1Defc">Marsaglia, George</a> (September 1963). <a href="https://apps.dtic.mil/sti/citations/AD0423993">Generating a Variable from the Tail of the Normal Distribution</a> (Technical report). Boeing Scientific Research Labs. Mathematical Note No. 322, DTIC accession number AD0423993. <a href="https://web.archive.org/web/20140910003249/http://www.dtic.mil/docs/citations/AD0423993">Archived</a> from the original on September 10, 2014 – via <a href="/facts/Defense_Technical_Information_Center/E7Vq7hhD">Defense Technical Information Center</a>.</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1">McFarland, Christopher D. (24 June 2015). "A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers". Journal of Statistical Computation and Simulation. 86 (7): 1281–1294. arXiv:1403.6870. doi:10.1080/00949655.2015.1060234. Note that the Bitbucket repository mentioned in the paper is no longer available and the code is now at https://github.com/cd-mcfarland/fast_prng <a href="https://mcfarlandlab.org/publication/mcfarland-2016/" target="_blank">https://mcfarlandlab.org/publication/mcfarland-2016/</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></li>
</ol>

Ziggurat algorithm open-in-new

Ziggurat algorithm