This is the third in a multi-part series on the MATLAB random number generators. MATLAB has used variants of George Marsaglia's ziggurat algorithm to generate normally distributed random numbers for almost twenty years.

Contents

The Ziggurat

It's important to have a memorable name for an algorithm. Ziggurats are ancient Mesopotamian terraced temple mounds that, mathematically, are two-dimensional step functions. A one-dimensional ziggurat underlies Marsaglia's algorithm. I'm not sure if we would still be using the algorithm today if Marsaglia had called it the "step function algorithm". The probability density function, or pdf, of the normal distribution is the bell-shaped curve $$ f(x) = c e^{-x^2/2} $$ where $c = 1/(2\pi)^{1/2}$ is a normalizing constant that we can ignore. If we were to generate random points $(x,y)$, uniformly distributed in the plane, and reject any of them that do not fall under this curve, the remaining $x$'s form our desired normal distribution. The ziggurat algorithm covers the area under the pdf by a slightly larger area with $n$ sections. The following figure has $n = 6$; the actual code we use today has $n = 256$. The choice of $n$ affects the speed, but not the accuracy of the code.

   z = zigplot(6)

z =

         0
    0.8288
    1.1713
    1.4696
    1.7819
    2.1761

The top $n-1$ sections of the ziggurat are rectangles. The bottom section is a rectangle together with an infinite tail under the graph of $f(x)$. The right-hand edges of the rectangles are at the points $z_k$ shown with circles in the figure. With $f(z_1) = 1$ and $f(z_{n+1}) = 0$, the height of the $k$ th section is $f(z_k) - f(z_{k+1})$. The key idea is to choose the $z_k$'s so that all $n$ sections, including the unbounded one on the bottom, have the same area. There are other algorithms that approximate the area under the pdf with rectangles. The distinguishing features of Marsaglia's algorithm are the facts that the rectangles are horizontal and have equal areas.

Initialization

For a specified number, $n$, of sections, it is possible to solve a transcendental equation to find $z_n$, the point where the infinite tail meets the last rectangular section. In our picture with $n = 6$, it turns out that $z_n = 2.18$. In an actual code with $n = 256$, $z_n = 3.6542$. Once $z_n$ is known, it is easy to compute the common area of the sections and the other right-hand endpoints $z_k$. It is also possible to compute $\sigma_k = z_{k-1}/z_k$, which is the fraction of each section that lies underneath the section above it. Let's call these fractional sections the core of the ziggurat. The right-hand edge of the core is the dotted line in our figure. The remaining portions of the rectangles, to the right of the dotted lines in the area covered by the graph of $f(x)$, are the tips. The computation of the $z_k$ 's and $\sigma_k$ 's is done only once and the values included in the header of the source code.

Central algorithm

With this initialization, normally distributed random numbers can be computed very quickly. The key portion of the code computes a single random integer, $j$, between $1$ and $n$ and a single uniformly distributed random number, $u$, between $-1$ and $1$. A check is then made to see if $u$ falls in the core of the $j$ th section. If it does, then we know that $u z_j$ is the $x$-coordinate of a point under the pdf and this value can be returned as one sample from the normal distribution. The code looks something like this.

  j = randi(256);
  u = 2*rand-1;
  if abs(u) < sigma(j)
     r = u*z(j);
  else
     r = tip_computation
  end

Most of the $\sigma_j$'s are greater than 0.98, and the test is true over 98.5% of the time. One normal random number can usually be computed from one random integer, one random uniform, an if-test, and a multiplication. The point determined by $j$ and $u$ will fall in a tip region less than 1.5% of the time. This happens if $j = 1$ because the top section has no core, if $j$ is between $2$ and $n-1$ and the random point is in one of the tips, or if $j = n$ and the point is in the infinite tail. In these cases, the more costly tip computation is required.

Accuracy

It is important to realize that, even though the ziggurat step function only approximates the probability density function, the resulting distribution is exactly normal. Decreasing $n$ decreases the amount of storage required for the tables and increases the fraction of time that extra computation is required, but does not affect the accuracy. Even with $n = 6$, we would have to do the more costly corrections over 25% of the time, instead of less than 1.5%, but we would still get an exact normal distribution.

Variations

Details of the ziggurat implementation have varied over the years. The original 1984 paper by Marsaglia and Tsang included a fairly elaborate transformation algorithm for handling the tips. We used this algorithm for several releases of MATLAB and I reproduced its behavior in the program randntx in Numerical Computing with MATLAB. That method and that code are now obsolete. The 2000 paper by Marsaglia and Tsang available at the jstatsoft link given below has a simpler rejection algorithm for use in the tips. We have been using this in more recent releases of MATLAB, including current ones.

Underlying Uniform Generator

Marsaglia and Tsang were anxious to make their normal generator as fast as their uniform generator. But they were a little too anxious. Their original code made one call to a 32-bit uniform generator. They used the high order 7 bits for the vertical index $j$ into the ziggurat and then reused all 32 bits to get the horizontal abscissa $u$. Later investigators, including Jurgen Doornik, found this correlation led to failures of certain statistical tests. We now make two calls to the 32-bit Mersenne Twister generator (during sequential computation). We take 8 bits to get $j$ and then 52 of the remaining 56 bits to get a double precision $u$. How does this affect the timing? Allocate a long vector.

   clear
   m = 25e6
   x = zeros(m,1);

m =

    25000000

Generate a random uniform vector and a random normal vector. Then compare the two execution times.

   tic, x = rand(m,1); tu = toc
   tic, x = randn(m,1); tn = toc
   ratio = tu/tn

tu =

    0.3416


tn =

    0.4520


ratio =

    0.7558

So, random uniform execution time is about three-quarters of the random normal execution time.

Acknowledgements

Thanks to Peter Perkins for his help on the entire series about the MATLAB random number generators. This post on the ziggurat is adapted from section 9.3 of Numerical Computing with MATLAB.

References

George Marsaglia and W. W. Tsang, "A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions." SIAM Journal on Scientific and Statistical Computing 5, 349-359, 1984. <http://epubs.siam.org/doi/pdf/10.1137/0905026> George Marsaglia and W. W. Tsang, "The ziggurat method for generating random variables." Journal of Statistical Software 5, 1-7, 2000. <http://www.jstatsoft.org/v05/i08> Jurgen A. Doornik, "An Improved Ziggurat Method to Generate Normal Random Samples." PDF, Nutfield College, Oxford, 2005. <http://www.doornik.com/research/ziggurat.pdf> Cleve Moler, Numerical Computing with MATLAB, SIAM 2004, <http://dx.doi.org/10.1137/1.9780898717952> Also available at: <https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/random.pdf>

Get the MATLAB code Published with MATLAB® R2014a