The Ziggurat Random Normal Generator
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.
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.
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.
댓글
댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.