Ziggurat Revisited

Introduction

Generating random number for use in simulation is a classic topic in scientific computing and about as old as the field itself. Most algorithms concentrate on the uniform distribution. Its values can be used to generate randomly distributed values from other distributions simply by using the relevant inverse function.

Regarding terminology, all computational algorithms for generation of numbers are by definition deterministic. Here, we follow standard convention and refer to such numbers as as they can always be recreated given the seed value for a given sequence. We are not concerned in this note with numbers (also called low-discrepnancy sequences). We consider this topic to be subset of pseudo-random numbers subject to distributional constraints. Without lack of generality, we will henceforth drop the qualifier when refering to random numbers.

Due to its importance for many modeling tasks, the Normal distribution has also been studied extensively in order to find suitable direct algorithms which generate stream of normally distributed (pseudo) random numbers. A well-known examples for such algorithms includes the method by Box and Muller . A useful recent survey of the field is provided by .

In an important paper, introduced the Ziggurat method. Since its initial publication, this algorithm has become reasonably popular as it provides a very useful combination of both excellent statistical properties and execution speed. conclude their survey by saying that

This paper reexamines the Ziggurat method, provides a new and portable C++ implementation, and compares it to several other Open Source implementations of the underlying algorithm by applying three different statistical tests.

Ziggurat

This sections briefly discusses the key papers related to Ziggurat.

Marsaglia and Tsang

introduced the Ziggurat method. It provides a fast algorithm for generating both normally and exponentially distributed random numbers. The original paper also contains a corresponding implementation in the C language.

The listing in Figure shows this initial implementation. We have removed the code for generating exponentially distributed random numbers as well a comment header from the listing to keep the display more compact. The full listing is also included in the package for reference.

As can be seen from Figure , the Ziggurat algorithm is implemented in bare-bones C code using a number of macros, and two helper functions. The helper functions allow setting a seed, and deal with parameter updates needed in about 2.5% of cases. The actual core component—the function to draw a random number distributed according to thstandard normal distribution—is provided by the macro . Needless to say, using C macros is no longer cosidered . Possible side effects include inadvertent changes in globally visible variables, as well as possible bugs from the macro evaluation.

A more important concern is that this implementation uses types, and explicit bit mapping operations. The code was originally developed for 32-bit operating systems where and are typically four bytes (or 32 bits) wide. Hence the code does not produce correct results on a 64-bit operating system as (signed and unsigned) types are typically eight bytes (or 64 bits) wide (whereas is still 32 bits).

Our modified version introduced below overcomes both issues.

Leong, Zhang et al

show in a comment that the Ziggurat method as introduced by suffers from another weakness due to the generator (by Marsaglia). The authors show via a χ2-test that the generator has a short period of about 232 − 1, or the four byte limit. Replacing it with the generator (also by Marsaglia) improves the performance beyond this limit.

#define MWC  ((znew<<16)+wnew )
#define SHR3 (jz=jsr, jsr^=(jsr<<13), \
        jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define CONG (jcong=69069*jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)

#define RNOR (hz=KISS, iz=hz&127, \
        (fabs(hz)<kn[iz]) ? hz*wn[iz] : nfix())

Following , Ziggurat code should use an improved uniform generator. Choices are either KISS as suggested initially, or another trusted (and fast) uniform generator such as the Mersenne Twister . Other Open Source implementations (such as the ones discussed below) frequently use the Mersenne Twister as the souce of uniformly distributed random numbers.

Voss’s implementation in GNU GSL

provided another Ziggurat implementation for use in the GNU Scientific Library or GSL . It uses the Mersenne Twister generator by which also avoids the issue identified by in which the originally-used uniform generator had too short a cycle.

notes two differences between his implementation and the original work by . First, he uses only 128 instead of 256 steps which reduces the memory requirements and computational cost at a possible (though presumably minor) loss of precision. Second, he uses an exponential distribution with tail rejection for the base strip, which is motivated by a simpler implementations. Both of these aspects could have implications for the statistical properties of the generator. Voss also appears to be unaware of the work by , yet sidesteps the issue raised by these author by relying on the Mersenne Twister generator.

Here, the implementation from the current GSL sources and file is used, and adapted to the class strucuture detailed in section .

Gretl

The Gretl econometrics program contains another Open Source implementation of the Ziggurat algorithm. The code credits the implementation by described above. review the Gretl implementation and performance of Ziggurat and find it to be satisfactory.

We use the implementation from the file from the current Gretl sources and adapt to the class strucuture detailed in section .

QuantLib

The QuantLib library for quantitative finance contains another open source implementation of Ziggurat. It is provided in the files and . As part of the experimental section, it is made available for further study and use, but not yet part of the default build. As before, we integrate this source into the class structure used here.

Speed

R Generators

Before comparing the speed of the different Ziggurat implementations, it is also illustrative to compare the different R generators. Figure provides a comparison.

We see that the Box-Muller generator is the slowest by some margin. However, both Kinderman-Ramage and Ahrens-Dieter are faster than the Inversion method chosen as the default in R. So even before considering Ziggurat generators, R users could reap a speed benefit simply by calling or .

Ziggurat Generators

All Ziggurat generators are significantly faster the the default generator in R which uses a inversion method. Among the Ziggurat generators, we notice that approaches using an external uniform number generator (GNU GSL, GNU Gretl, QuantLib) are all slower than our compact and self-contained implementation which is seen as the fastest method.

Accuracy

Standard Test for Uniform RNG draws

Test for random number generators are often focussed on the case of uniform generators which are the most common type of generators. As detailed for example in , a test proceeds as follow:

Here, we can apply this test for first converting the N(0, 1) distributed values produced by the given Ziggurat implementation to U(0, 1) distributed values by using the inverse of the normal distribution. We are then ready to simulate and test. Figure below displays Q-Q plots for the empirical distribution against the uniform, and displays the p-values of a Kolmogorov-Smirnov as well as a Wilcoxon test.

We see that five of the six generators pass the test. In the case of the original generator, we can see the departure from the expected diagonal clearly once we draw more than 4.2 × 109 numbers (which is the limit of the representation of an unsigned four-byte nunber). However, only the Kolmogorov-Smirnov test can formally reject; the Wilcoxon test appear to lack sufficient power in this setting. The QuantLib generator is seen as suspicious which p-value just below a conventional rejection level.

Normal Test for Normal RNG draws

We can propose a simpler variant of the test outlined in the previous section. As the random numbers we are drawing are following a N(0, 1) distribution, the sum of their values follows a N(0, n) distribution. This allows us to skip one inversion step:

Results, shown in Figure , are qualitatively similar to the result discussed above. Kolmogorov-Smirnov rejects for the generator. However, we note that the Wilcoxon test now has a lower p-value—we would now reject at conventional test levels. The QuantLib implementation is now rejected by the Wilcoxon test but not the Kolmogorov-Smirnov.

χ2 test

Another test variant is the χ2 test which was also used by . The basic idea is as follow:

As can be seen in Figure , the original proposal by fails as was shown by . All other tests pass again.

C++ Implementation

Preceding work by and also contained implementations in the C language. These versions were implemented in just a few lines, and used idioms common to C programmers such as macros and global variables.

C++ programming style permits encapsulation in order to avoid possible collisions and side-effects. Moreover, by using a modest amount of object-oriented programming we can use a class structure with a common base class to express commonalities between the implementations. The following code segment shows the virtual base class used here.


#ifndef RcppZiggurat__Zigg_h
#define RcppZiggurat__Zigg_h

#include <cmath>
#include <stdint.h>         // or cstdint (C++11)

namespace Ziggurat {

   class Zigg {
   public:
      virtual ~Zigg() {};
      virtual void setSeed(const uint32_t s) = 0;
      // no getSeed() as GSL has none
      virtual double norm() = 0;
   };
}

#endif

As shown in the preceding code segment, we provide two user-accessible functions to obtain a normal random deviate, and to set the seed, respectively. The actual implementation uses portable types such as , an unsigned 32-bit integer provided by the C header file , which provides correct and identical results on both 32-bit and 64-bit operating systems.

Each actual implementation can then encapsulate its state variable as a private variable inaccessible to other functions. Such a small core to each class also makes it feasible to provide a Ziggurat generator in each thread in a parallel execution framework.

Our Ziggurat implementation has no external dependencies and can therefore be used in other projects. The testing framework used for this note has a single dependency on the GNU GSL as the generator by is used via its GSL implementations. The generator and testing framework in the corresponding R package have a build-dependency on R, and are of course accessed by R. But the generator discussed here could equally well be used in standalone programs or with other scripting languages.

In the package, the Rcpp Attributes feature of the C++/R integration package are used to access instances of the corresponding generator class.


#include <Rcpp.h>

#include <Ziggurat.h>

static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
Rcpp::NumericVector zrnorm(int n) {
    Rcpp::NumericVector x(n);
    for (int i=0; i<n; i++) {
        x[i] = zigg.norm();
    }
    return x;
}

// [[Rcpp::export]]
void zsetseed(unsigned long int s) {
    zigg.setSeed(s);
    return;
}

In this particular reference implementation, we have chosen a namespace for the entire project. Within this namespace, we opted to provided an extra namespace layer for each generator as some of these generators still use global variables—which are therefore shielded in their own namespace. For example, for the generator, we use . Next is the name of the actual class—which in the case of the reference implementation shown above is once again leading to the triple use of the term. Actual deployment of the Ziggurat generator (without comparison to other implementations and concerns about interaction between variables, particularly for the older implementations having global variables) can of course be used with a single namespace.

The remainder shows how two functions and are provided via the attribute as described in the Rcpp Attributes vignette of the package .

Conclusion

This note describes the package and its new implementation of the Ziggurat generator for normally distributed random numbers. The package is implemented in a way which is portable so that it can be used on 32-bit and 64-bit operating systems, filling a gap left by the original implementations .

By embedding the code in a simple C++ class structure, we can ease testing and comparison of several variants of the algorithm. Our note reconfirmed the findings by of a short cycle due to to the use of an inferior uniform generator. Replacing the generator leads to better performance.

We suggest a new test for random number generators producing deviates distributed according to the standard normal distribution by adapting and simplifying an existing test framework for uniform deviates. Both tests confirm the (previously documented) failure of the original Ziggurat proposal but do not find a problem with any of the new implementations (apart from the still-experimental QuantLib generator).

A key motivation for this work has been a desire to improve the speed of creating standard-normally distributed random numbers in R. We find Ziggurat to be faster than the existing implementations, and hope that this generator will be of use to those generating large numbers of draws.