Technical Articles and Newsletters |

By Cleve Moler, MathWorks

MATLAB's function `rand(m,n)` generates an ` m`-by-

Almost all algorithms for generating normally distributed random numbers are based on transformations of uniform distributions. The simplest way to generate an ` m`-by-

`sum(rand(m,n,12),3) - 6`

This works because `R = rand(m,n,p)` generates a three-dimensional uniformly distributed array and `sum(R,3)` sums along the third dimension. The result is a two-dimensional array with elements drawn from a distribution with mean *p/2* and variance *p/12* that approaches a normal distribution as *p* increases. If we take *p =* 12, we get a pretty good approximation to the normal distribution and we get the variance to be equal to one without any additional scaling. There are two difficulties with this approach. It requires twelve uniforms to generate one normal, so it is slow. And the finite *p* approximation causes it to have poor behavior in the tails of the distribution.

Older versions of MATLAB - before MATLAB 5 - used the *polar* algorithm. This generates two values at a time. It involves finding a random point in the unit circle by generating uniformly distributed points in the `[0,1]x[0,1]` square and rejecting any points outside of the circle. Points in the square are represented by vectors with two components. The rejection portion of the code is:

`r = 2;`

while r > 1

u = 2*rand(2,1)-1

r = u'*u

end

For each point accepted, the polar transformation

`v = sqrt(-2*log(r)/r)*u`

produces a vector with two independent normally distributed elements. This algorithm does not involve any approximations, so it has the proper behavior in the tails of the distribution. But it is moderately expensive. Over 21% of the uniform numbers are rejected when they fall outside of the circle and the square root and logarithm calculations contribute significantly to the cost.

Beginning with MATLAB 5, and continuing with MATLAB 6, `randn` uses a sophisticated table lookup algorithm developed by George Marsaglia of Florida State University. Marsaglia calls his approach the "ziggurat" algorithm. Ziggurats are ancient Mesopotamian terraced temple mounds that, mathematically, are two-dimensional step functions. A one-dimensional ziggurat underlies Marsaglia's algorithm.

Marsaglia has refined his ziggurat algorithm over the years. An early version is described in Knuth's classic *The Art of Computer Programming,* volume II. MATLAB's version is described by Marsaglia and W. W. Tsang in the *SIAM Journal of Scientific and Statistical Programming*, volume 5, 1984. A Fortran subroutine is discussed briefly in the Kahaner, Moler, and Nash text book, *Numerical Methods and Software*, and is available from the MathWorks FTP site. A recent version is available in the online electronic journal, *Journal of Statistical Software*. We describe this recent version here because it is the most elegant. The version actually used in MATLAB is more complicated, but is based on the same ideas and is just as effective.

The probability density function, or *pdf*, of the normal distribution is the bell-shaped curve*F(x) = c exp(-x ^{2}/2)*

where *c=1/(2π) ^{1/2}* is a normalizing constant that we can ignore. If we generate random points

The ziggurat algorithm covers the area under the pdf by a slightly larger area with *n* sections. The accompanying picture has ` n` = 8; the actual code uses

For a specified number, ` n`, of sections, it is possible to solve a transcendental equation to find

After the initialization, normally distributed random numbers can be computed very quickly. The key portion of the code computes a single random integer, ` k`, between 1 and

`k = ceil(128*rand);`

u = 2*rand-1;

if abs(u) < sigma(k)

v = u*x(k);

return

end

Most of the *σ _{k}'s* are greater than 0.98, and the test is true over 97% of the time. One normal random number can usually be computed from one random integer, one random uniform, an if-test and a multiplication. No square roots or logarithms are required.

The point determined by ` k` and

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

With this algorithm, MATLAB 6 can generate normally distributed random numbers as fast as it can generate uniformly distributed ones. In fact, MATLAB on my 800 MHz Pentium III laptop can generate over 10 million random numbers from either distribution in less than one second.

MATLAB's built-in function `randn` uses two unsigned 32-bit integers, one for the seed in the uniform generator and one as the shift register for the integer generator. The code that processes points outside the core of the ziggurat accesses these two generators independently, so the period of the overall generator is 2^{64}. At 10 million normal deviates per second, `randn` will take over 58,000 years before it repeats.

Published 2001