ARVAG - Automatic nonuniform Random VAriate Generation

W. Hörmann, J. Leydold, and G. Derflinger
Automatic Nonuniform Random Variate Generation
Series: Statistics and Computing
2004, X, 442pp, Hardcover
ISBN: 3-540-40652-2
Springer, Berlin

About this book   |   Table of Contents   |   Introduction   |   Errata   |   Demo of Algorithms


Non-uniform random variate generation is a small field of research somewhere between mathematics, statistics and computer science. It started in the fifties of the last century in the "stone-age" of computers. Since then its development has mainly been driven by the wish of researchers to solve generation problems necessary to run their simulation models. Also the need for fast generators has been an important driving force. The main mathematical problems that have to be solved concern the distribution of transformed random variates and finding tight inequalities. Also implementing and testing the proposed algorithms has been an important part of the research work. A large number of research papers in this field has been published in the seventies and early eighties. The main bibliographical landmark of this development is the book of Devroye (1986a), that is commonly addressed as the "bible" of random variate generation. We can certainly say that random variate generation has become an accepted research area considered as a subarea of statistical computing and simulation methodology. Practically all text-books on discrete event simulation or Monte Carlo methods include at least one chapter on random variate generation; within simulation courses it is taught even to undergraduate students.

More important is the fact that random variate generation is used by lots of more or less educated users of stochastic simulation. Random variate generation code is found in spreadsheets and in expensive discrete event simulation software and of course in a variety of programming languages. Probably many of these users do not bother about the methods and ideas of the generation algorithms. They just want to generate random variates with the desired distribution. "The problems of random variate generation are solved" these people may say. And they are right as long as they are only interested in popular standard distributions like normal, gamma, beta, or Weibull distributions.

Why Universal Random Variate Generation?

The situation changes considerably if the user is interested in non standard distributions. For example, she wants to simulate models that include the generalized inverse Gaussian distribution, the hyperbolic distribution, or any other distribution that is less common or even newly defined for her purpose. Then the user had two possibilities: She could either find some code (or a paper) dealing with the generation of random variates from her distribution, or she needed some knowledge on random variate generation (or find an expert) to design her own algorithm. Today the user has a third possibility: She can find a universal generator suitable for her distribution, perhaps in our C library UNU.RAN (Universal Non-Uniform RANdom variate generation). Then she can generate variates without designing a new generator; a function that evaluates e.g. the density of the distribution or the hazard rate is sufficient.

This book is concentrating on the third possibility. We present ideas and unifying concepts of universal random variate generation, and demonstrate how they can be used to obtain fast and robust algorithms. The book is presenting the first step of random variate generation, the design of the algorithms. The second step, i.e. the implementation of most of the presented algorithms, can be found in our C library UNU.RAN (Universal Non-Uniform RANdom number generators, Leydold, Hömann, Janka, and Tirler, 2000). As some of these algorithms are rather long it is not possible, and probably not desirable, to present all implementation details in this book as they would hide the main ideas.

What Is Non-Uniform Random Variate Generation?

Usually random variates are generated by transforming a sequence of independent uniform random numbers on (0,1) into a sequence of independent random variates of the desired distribution. This transformation needs not be one-to-one. We assume here, as it is generally done in the literature, that we have an ideal source of uniform random numbers available. An assumption which is not too unrealistic if we think of fast, modern uniform random number generators that have cycle lengths in the order 2 raised to the power of several hundreds or even thousands and equidistribution property up to several hundred dimensions, e.g., Matsumoto and Nishimura (1998) or L'Ecuyer (1999); see also the pLab website maintained by Hellekalek (2002) for further links. A collection of many published uniform random number generators -- good ones and bad ones -- is compiled by Entacher (2000).

Given that (ideal) source of uniform random numbers, the well known inversion, (acceptance-) rejection and decomposition methods can be used to obtain exact random variate generation algorithms for standard distributions. We do not want to give a historical overview here but it is remarkable that the rejection method dates back to von Neumann (1951). Later refinements of the general methods were developed mainly to design fast algorithms, if possible with short code and small memory requirements. For the normal distribution compare e.g. Box and Muller (1958), Marsaglia, MacLaren, and Bray (1964), Ahrens and Dieter (1972,1973,1988a), and Kinderman and Ramage (1976). For the gamma distribution see e.g. Cheng (1977), Schmeiser and Lal (1980), Cheng and Feast (1980), and Ahrens and Dieter (1982). In the books of Devroye (1986a) and Dagpunar (1988) you can find an impressive number of references for articles dealing with random variate generation for standard distributions. The techniques and general methods we describe in this book are very closely related to those developed for standard distributions. We know what we owe to these "pioneers" in random variate generation. Even more as we have started our research in this field with generators for standard distributions as well (see e.g. Hö:rmann and Derflinger, 1990a). However, in this book we try to demonstrate that several of these main ideas can be applied to build universal generators for fairly large distribution families. Thus we do not concentrate on standard distributions but try to identify large classes of distributions that allow for universal generation. Ideally these classes should contain most important standard distributions. So the reader will come across quite a few standard distributions as we use them as examples for figures and timings. They are best suited for this purpose as the reader is familiar with their properties. We could have included fairly exotic distributions as well, as we have done it in some of the empirical comparisons.

What Is a Universal Generator?

A universal (also called automatic or black-box) generator is a computer program that can sample from a large family of distributions. The distributions are characterized by a program that evaluates (e.g.) the density, the cumulative distribution function, or the hazard rate; often some other information like the mode of the distribution is required. A universal generator typically starts with a setup that computes all constants necessary for the generation, e.g. the hat function for a rejection algorithm. In the sampling part of the program these constants are used to generate random variates. If we want to generate a large sample from a single distribution, the setup part is executed only once whereas the sampling is repeated very often. The average execution time of the sampling algorithm to generate one random variate (without taking the setup into account) is called marginal execution time. Clearly the setup time is less important than the marginal execution time if we want to generate a large sample from a single distribution.

Fast universal generators for discrete distributions are well known and frequently used. The indexed search method (Chen and Asau, 1974) and the alias method (Walker, 1974, 1977) can generate from any discrete distribution with known probability vector and bounded domain. For continuous distributions the design of universal generators started in the eighties and is mainly linked with the name of Luc Devroye. During the last decade research in this direction was intensified (see e.g. Gilks and Wild, 1992; Hömann, 1995; Ahrens, 1995; Leydold, 2000) and generalized to random vectors (see e.g. Devroye 1997; Leydold and Hö:rmann, 1998). However it seems that these developments are little known and hardly used. When we have started our UNU.RAN project in 1999 we checked several well known scientific libraries: IMSL, Cern, NAG, Crand, Ranlib, Numerical recipes, and GSL (Gnu Scientific Library). Although all of them include quite a few random variate generation algorithms for continuous standard distributions we were not able to find a single universal generator for continuous distributions in any of them. And the situation is similar for random-vector generation methods. This has been a main motivation for us to start our project with the aim to write both a book and a C library to describe and realize theory and practice of universal random variate and random vector generation.

Why We Have Written this Book?

We are convinced that universal random variate generation is a concept of greatest practical importance. It also leads to nice mathematical theory and results. Neither the mathematics nor the implemented algorithms have been easily accessible up to now, as there is -- up to our knowledge -- no book and no software library available yet that includes the main concepts of universal random variate generation.

Implementation of Universal Algorithms

One main problem when implementing universal algorithms for continuous distributions in a software library is certainly the application programming interface (API). Considering generators for a fixed distribution everything is simple. Any programmer can easily guess that the following C statements are assigning a realization from a uniform, a standard normal, and a gamma(5) random variate to the respective variables.

  x  = random();
  xn = randnormal();
  xg = randgamma(5.);
For a universal method the situation is clearly more complicated. The setup often is very expensive compared to the marginal generation time and thus has to be separated from the sampling part of a universal algorithm. Moreover, to run such a setup routine we need a lot more parameters. In a routine like randgamma(5.) all required information is used to build the algorithm and is thus contained implicitly in the algorithm. For black-box algorithms we of course have to provide this information explicitly. This may include -- depending on the chosen algorithm -- the density, its derivative and its mode, its cumulative distribution function, its hazard rate, or similar data and functions to describe the desired distribution. Furthermore, all the black-box algorithms have their own parameters to adjust the algorithm to the given distribution. All of this information is needed in the setup routine to construct a generator for the distribution and to store all necessary constants. The sampling program then uses these constants to generate random variates.

There are two very different general solutions for the implementation of automatic algorithms. First, we can make a library that contains both a setup routine and a sampling routine using an object-oriented design. The setup routine creates an instance of a generator object that contains all necessary constants. The sampling routine is using this generator object for generation. If we want to sample from several different distributions in one simulation we can create instances of such generator objects for all of them and can use them for sampling as required. In our library UNU.RAN we have implemented this idea. For further details see Sect. 8.1.

Our experiences with the implementation of universal algorithms in a flexible, reliable, and robust way results in rather large computer code. As the reader will find out herself, the complexity of such a library arises from the setup step, from parts performing adaptive steps, and (especially) from checking the data given by the user, since not every method can be used for every distribution. The sampling routine itself, however, is very simple and consists only of a few lines of code. Installing and using such a library might seem too tedious for "just a random number generator" at a first glance, especially when only a generator for a particular distribution is required. As a solution to this problem we can use universal methods to realize the concept of an automatic code generator for random variate generation. In this second approach we use the constants that are computed in the setup to produce a single piece of code in a high level language for a generator of the desired distribution. Such a code generator has the advantage that it is also comparatively simple to generate code for different programming languages. Moreover, we can use a graphical user interface (GUI) to simplify the task of obtaining a generator for the desired distribution for a practitioner or researcher with little background in random variate generation. We also have implemented a proof of concept study using a web based interface. It can be found at Currently, program code in C, FORTRAN, and Java can be generated and downloaded. For more details see Leydold, Derflinger, Tirler, Hörmann (2003).

It should be clear that the algorithm design of the setup and the sampling routine remain the same for both possible implementation concepts. Thus we will not consider the differences between these two approaches throughout the book. The only important difference to remember is that the speed of the setup is of little or no concern for the code generator whereas it may become important if we use generator objects.

Theoretical Concepts in the Book

We have spoken quite a lot about algorithms and implementation so far but most of the pages of the book are devoted to the theoretical concepts we need for designing random variate generation algorithms. Typically we are concerned with the following problems:

Chapters of the Book

Chapter 2 presents the most important basic concepts of random variate generation: Inversion, rejection, and composition for continuous random variates. Thus it is crucial to the rest of the book as practical all of the algorithms presented in the book depend on one or several of these principles. Chapter 3 continues with the basic methods for generating from discrete distributions, among them inversion by sequential search, the indexed search and the alias method.

Part II of the book deals with continuous univariate distributions. Chapters 4, 5, and 6 present three quite different approaches to design universal algorithms by utilizing the rejection method. Chapter 7 realizes the same task using numerical inversion. Chapter 8 compares different aspects of the universal algorithms presented so far including our computational experiences when using the UNU.RAN library. It also describes the main design of the UNU.RAN programming interface that can be used for generating variates from discrete distributions and from random vectors as well. Part II closes with Chap. 9 that collects different special algorithms for the case that the density or cumulative distribution function of the distribution is not known.

Part III consists only of Chap. 10. It explains recent universal algorithms for discrete distributions, among them the indexed search method for distributions with unbounded domain and different universal rejection methods.

Part IV, that only consists of Chap. 11, presents general methods to generate random vectors. It also demonstrates how the rejection method can be utilized to obtain universal algorithms for multivariate distributions. The practical application of these methods are restricted to dimensions up to about ten.

Part V contains different methods and applications that are closely related to random variate generation. Chapter 12 collects random variate generation procedures for different situations where no full characterization of the distribution is available. Thus the decision about the generation procedures is implicitly also including a modeling decision. Chapter 13 (co-authored by M. Hauser) presents very efficient algorithms to sample Gaussian time series and time series with non-Gaussian one-dimensional marginals. They work for time series of length up to one million. Markov Chain Monte Carlo (MCMC) algorithms have become frequently used in the last years. Chapter 14 gives a short introduction into these methods. It compares them with the random vector generation algorithms of Chap. 11 and discusses how MCMC can be used to generate iid. random vectors. The final Chap. 15 presents some simulation examples for financial engineering and Bayesian statistics to demonstrate, at least briefly, how some of the algorithms presented in the book can be used in practice.

Reader Guidelines

This book is a research monograph as we tried to cover -- at least shortly -- all relevant universal random variate generation methods found in the literature. On the other hand the necessary mathematical and statistical tools are fairly basic which should make most of the concepts and algorithms accessible for all graduate students with sound background in calculus and probability theory. The book mainly describes the mathematical ideas and concepts together with the algorithms; it is not closely related to any programming language and is generally not discussing technical details of the algorithm that may depend on the implementation. The interested reader can use the source code of UNU.RAN to see possible solutions to the implementation details for most of the important algorithms. Of course the book can also be seen as the "documentation" explaining the deeper mathematics behind the UNU.RAN algorithms. Thus it should also be useful for all users of simulations who want to gain more insight into the way random variate generation works.

In general the chapters collect algorithms that can be applied to similar random variate generation problems. Only the generation of continuous one-dimensional random variates with known density contains so much material that is was partitioned into the Chaps. 4, 5, and 6. Several of the chapters are relatively self contained. Nevertheless, there are some dependencies that lead to the following suggestions: For readers not too familiar with random variate generation we recommend to read Chap. 2 and Chap. 3 before continuing with any of the other chapters. Chapters 4 to 7 may be read in an arbitrary order; some readers may want to consult Chap. 8 first to decide which of the methods is most useful for them. Chapters 9, 12, and 13 are self contained, whereas some sections of Chaps. 10 and 11 are are generalizations of ideas presented in Chaps. 4 and 6. Chapter 14 is self contained, but to appreciate its developments we recommend to have a look at Chap. 11 first. In Chap. 15 algorithms of Chaps. 2, 4, 11 and 12 are used for the different simulations.

Course Outlines

We have included exercises at the end of most of the chapters as we used parts of the book also for teaching courses in simulation and random variate generation. It is obvious that the random variate generation part of a simulation course will mainly use Chaps. 2 and 3. Therefore we tried to give a broad and simple development of the main ideas there. It is possible to add the main idea of universal random variate generation by including for example Sect. 4.1. Selected parts of Chap. 15 are useful for a simulation course as well.

A special course on random variate generation should start with most of the material of Chaps. 2 and 3. It can continue with selected sections of Chap. 4. Then the instructor may choose chapters freely, according to her or the students' preferences: For example, for a course with special emphasis on multivariate simulation she could continue with Chaps. 11 and 14; for a course concentrating on the fast generation of univariate continuous distributions with Chaps. 5, 7, and 8. Parts of Chap. 15 can be included to demonstrate possible applications of the presented algorithms.

What Is New?

The main new point in this book is that we use the paradigm of universal random variate generation throughout the book. Looking at the details we may say that the rigorous treatment of all conditions and all implications of transformed density rejection in Chap. 4 is probably the most important contribution. The second is our detailed discussion of the universal generation of random vectors. The discussion of Markov chain Monte Carlo methods for generating iid. random vectors and its comparison with standard rejection algorithms is new.

Throughout the book we present new improved versions of algorithms: Among them are variants of transformed density rejection, fast numerical inversion, universal methods for increasing hazard rate, indexed search for discrete distributions with unbounded domain, improved universal rejection algorithms for orthomonotone densities, and exact universal algorithms based on MCMC and perfect sampling.

Practical Assumptions

All of the algorithms of this book were designed for practical use in simulation. Therefore we need the simplifying assumption of above that we have a source of truly uniform and iid. (independent and identically distributed) random variates available. The second problem we have to consider is related to the fact that a computer cannot store and manipulate real numbers.

Devroye (1986a) assumes an idealized numerical model of arbitrary precision. Using this model inversion requires arbitrary precision as well and is therefore impossible in finite time unless we are given the inverse cumulative distribution function. On the other hand there are generation algorithms (like the series method) that require the evaluation of sums with large summands and alternating signs, which are exact in the idealized numerical model of arbitrary precision.

If we consider the inversion and the series method implemented in a modern computing environment (e.g. compliant with the IEEE floating point standard) then -- using bisection -- inversion may be slow but it is certainly possible to implement it with working precision close to machine precision. On the other hand, the series method may not work as -- due to extinction -- an alternating series might "numerically" converge to wrong values. We do not know the future development of floating point arithmetic but we do not want to use an idealized model that is so different from the behavior of today's standard computing environments. Therefore we decided to include numerical inversion algorithms in this book as long as we can easily reach error bounds that are close to machine precision (i.e. about 10-10 or 10-12). We also include the series method but in the algorithm descriptions we clearly state warnings about the possible numerical errors.

The speed of random variate generation procedures is still of some importance as faster algorithms allow for larger sample sizes and thus for shorter confidence intervals for the simulation results. In our timing experiments we have experienced a great variability of the results. They depended not only on the computing environment, the compiler, and the uniform generator but also on coding details like stack variables versus heap for storing the constants etc. To decrease this variability we decided to define the relative generation time of an algorithm as the generation time divided by the generation time for the exponential distribution using inversion which is done by -log(1-random()). Of course this time has to be taken in exactly the same programming environment, using the same type of function call etc. The relative generation time is still influenced by many factors and we should not consider differences of less than 25%. Nevertheless, it can give us a crude idea about the speed of a certain random variate generation method.

We conclude this introduction with a statement about the speed of our universal generators. The relative generation time for the fastest algorithms of the Chaps. 4, 5, and 7 are not depending on the desired density and are close to one; i.e. these methods are about as fast as the inversion method for exponential variates and they can sample at that speed from many different distributions.

[ARVAG]     Wolfgang Hörmann and Josef Leydold   (October 21st, 2003) Research supported by FWF