Publication: Some Fundamental Properties of a Multivariate von Mises Distribution. Multivariate von Mises dis­tri­bu­tions can be used to describe components in mixture models for angular data. Here, we derive sufficient criteria for von Mises dis­tri­bu­tions to be unimodal. … read more

Blog: Strange random number generator. The default random number generator of Go seems to be of unclear provenance. … read more

Übersicht über die gängigen Wahr­schein­lich­keits­ver­tei­lun­gen. Definition, Er­war­tungs­wert, Varianz und Mo­men­ten­er­zeu­gen­de Funktion für verschiedene Wahr­schein­lich­keits­ver­tei­lun­gen. … read more

The Ziggurat Method for Generating Gaussian Random Numbers

By Jochen Voss, last updated 2014-06-11

Introduction

The Ziggurat method for generating random numbers is a special case of the rejection method, where the density of the proposal distribution looks like the silhouette of a Ziggurat. The method can be used to produce very fast random number generators for distributions like the normal distribution or the exponential distribution. The main advantage comes from the fact that for a high percentage of the generated numbers no slow floating point operations are necessary.

The technical details of the Ziggurat method are, for example, contained in the article The Ziggurat Method for Generating Random Variables by G. Marsaglia and W. W. Tsang.

On this page I provide an implementation of the Ziggurat method to generate Gaussian random numbers. I hope to see this function to be included in the GSL library at some point of time. Therefore I compare the speed of my implementation with the two functions for Gaussian random number generation, which currently are included in the GSL library.

Implementation

My implementation of the Ziggurat method can be downloaded below. The file gauss.c contains the implementation, the other files are only there for illustration.

There are some differences between my implementation and the one in the article cited above: 1) I use only 128 steps instead of 256 to decrease the amount of required tabulated data. 2) I use an exponential distribution together with rejection for the tails of the base strip to simplify the implementation. Both changes seem to have no significant performance impact.

Speed Measurements

Since the Ziggurat method needs 32-bit random integers, it works best for random number generators with gsl_rng_min=0 and gsl_rng_max=232-1. The following table shows the speed improvement for the different random number generators. For the default generator mt19937 the speed up factor is 3.9.

RNG Ziggurat Box-Muller ratio method improvement
borosh13 8.89 2.31 2.39 3.7
cmrg 2.84 1.34 1.3 2.1
coveyou 8.88 2.39 2.07 3.7
fishman18 1.23 0.882 0.863 1.4
fishman20 5.88 2.21 2.13 2.7
fishman2x 4.35 1.69 1.51 2.6
gfsr4 9.03 2.19 2.29 4.0
knuthran 6.47 2.21 2.27 2.8
knuthran2 0.662 0.527 0.499 1.3
lecuyer21 6.61 2.22 2.17 3.0
minstd 5.61 2.16 2.1 2.6
mrg 4.38 1.73 1.66 2.5
mt19937 8.28 2.01 2.14 3.9
mt19937_1999 8.28 2.02 2.14 3.9
mt19937_1998 8.27 2.01 2.12 3.9
r250 9.1 2.2 2.05 4.1
ran0 5.61 2.17 2.1 2.6
ran1 5.94 2.1 2.04 2.8
ran2 4.4 1.68 1.63 2.6
ran3 6.69 2.05 2.12 3.2
rand 8.47 2.37 2.37 3.6
rand48 8.14 0.887 0.84 9.2
random128-bsd 7.91 2.2 2.3 3.4
random128-glibc2 7.72 2.18 2.31 3.3
random128-libc5 7.94 2.2 2.31 3.4
random256-bsd 8.02 2.18 2.29 3.5
random256-glibc2 8 2.2 2.28 3.5
random256-libc5 7.72 2.18 2.29 3.4
random32-bsd 8.14 2.24 2.34 3.5
random32-glibc2 7.93 2.25 2.34 3.4
random32-libc5 8.19 2.24 2.34 3.5
random64-bsd 7.71 2.2 2.32 3.3
random64-glibc2 7.71 2.21 2.33 3.3
random64-libc5 7.68 2.2 2.32 3.3
random8-bsd 8.45 2.37 2.36 3.6
random8-glibc2 8.45 2.37 2.36 3.6
random8-libc5 8.47 2.37 2.37 3.6

Table 1. This table compares the speed of Gaussian random number generation for different methods and different random number generators. The first column gives the GSL name for the underlying random number generator. Columns 2, 3, and 4 give the speed (in million generated numbers per second) of the Ziggurat method gsl_ran_gaussian_ziggurat, the Box-Muller method gsl_ran_gaussian and the ratio method gsl_ran_gaussian_ratio_method. All measurements were done on a Pentium 4 with 2 GHz, the code was compiled with gcc -O2. The last column gives the ratio of the values for the Ziggurat method and for the better one of the other two methods.

Statistical Tests

In order to validate my implementation I performed several statistical tests on my random number generator. The program to perform these tests can be downloaded below. The (a little bit arbitrary) list of tests used is:

All three methods, gsl_ran_gaussian_ziggurat, gsl_ran_gaussian and gsl_ran_gaussian_ratio_method, pass these tests.

Download

References

Copyright © 2014, Jochen Voss. All content on this website (including text, pictures, and any other original works), unless otherwise noted, is licensed under a Creative Commons Attribution-Share Alike 3.0 License.