As simulation arid Monte Carlo continue to play an increasing role in statistical research, careful attention must be given to problems which arise in implementing and documenting collect ions of random number generators. This paper examines the value of theoretical as well as empirical evidence in establishing the quality of generators, the selection of generators to comprise a good basic set, the techniques and efficiency of implementation, and the extent of documentation. Illustrative examples are drawn from various current sources.

David C. Hoaglin

English

Research Papers in Economics

NBER WORKING PAPER SERIES
IMPLEMENTING AND DOCUMENTING
RANDOM NUMBER GENERATORS
David C. Hoaglin*
Working Paper No. 75
COMRJTERRESEARCH CEN'ITR FOR ECONOMICS AND MANEL tENT SCIENCE
National Bureau of Economic Research,Inc.
575 TechnologySquare
Cambridge,Massachusetts 02139
March 1975
Preliminary: not for quotation
NBERworking papers are distributed informally and in limited
numbers for comments only. They should not be quoted without
written permission.
This report has not undergone the review accorded official NBER
publications; in particular, it has not yet been submitted for
approval by the Board of Directors.
NBER Computer Research Center and Harvard University, Department
of Statistics. Research supported in part by National Science
Foundation Grant GJ-1l54X3 to the National Bureau of Economic'
Research, Inc.Abstract
As simulation arid Monte Carlo continue to play an iricrasing role in statis-
tical research, careful attention must be given to problems which arise in
implementing and documenting collect ions of random number generators. This
paper examines the value of theoretical as well as empirical evidence in
establishing the quality of generators, the selection of generators to
comçxr'ise agood basic set, the techniques andefficiency of implementation,
and the extent of documentation. Illustrative examples are drawn from van-
out current sources.
.
.Contents
1. Introduction •1
2. Uniform Generators •1
3. A Basic Library •2
4. Implementation 2
5. Documentation 3
6. Sunniry 3
References 3BIPLUIEWrIN3 AND DOCUMDING
RANDOMNUMBER GENERA'IORS
1. INTRODUCTION
Across the fields of statistics arid coiñputer
science, from theoretical to applied, simulation
and !bnte Carlo continue to play a significant
role. The variety of clever applications is great,
but often it seems that the technical foundations
are shaky. The random-number generators, on which
this whole experimentation structure rests, are
still all too often incautiously selected, haphaz-
arly implemented, and inadequately documented.
Taking examples from arrong the available genera—
tors, algorithms, routines, and libraries, the
balance of this paper examines prevailing prac-
tices in selecting, implementing, and documenting
random—number generators and offers some reccinen-
dations.
First, however, we shouldlookfurtheratthe
questionof prevalence: How widely are simulation
and Tbnte Carlo used in statistical research?In
preparing a position paperonpublication of
computation-based results [8] David Andrews andI
hadoccasionto go carefully throughthe1973
volumesofBiometrikaand Journalof the American
Statistical Association, counting pap of vari-
ous kinds. We found that 20% of all papers in—
volved simulation results, and the individual per-
centages in the three hodies of papers (JASA
Applications, JASA Theory and Methods, and
Biometrika) derted surprisingly little from the
overall figure. Even without comparable data for
the computer science literature, the overall con-
clusion is clear: simulation is an important com-
ponent in quite a lot of research.
2. UNIFORMGENERATORS
A source of uniform ranaom numbers is at the heart
of aljrost all algorithms for generating non-
uniform distributions, so it deserves a lot of
attention and is the natural place to start. This
observationis hardly new, butinview of the
typicalqualityof. available generators it still needs hasis. Forexample, wherever there is a
computerwitha word size of 32bits, one is likely
tofindthepoorgeneratorRANDU [10] --itsper-
sistence has beenrernrkable.Ofcourse, this
generator(along with rrost others in comiron use) is
nultiplicative-corigruentjal,arid all such genera-
tors are well-3-iowntoproduceoutputsequences
whichhave regularstructure. Specifically, the
set of all n-tuples (Xi,Xi1,..., Xj+nl) fO alattice in Euclidean n-space [13]. For some
high-accuracymultidimensional calculations this
type of defect may render all congr'ueritial genera-
tors undesirable. Various schemes for permuting,
raridomly sFniffling, or otherwise modifying the out-
put of such basic generators offer reasonable
improvement, but there is still much to be learned.
Returning now to the basic congruentialgenerators,
weshould be aware of one important argument in
theirfavor: they are the easiest to analyze
theoretically. As a result we can determine (over
the full period of the generator) several indica-
tive properties of these generators and thus have
a much clearer picture of what we can expect of
them and what we cannot. The lattice structure of
the generators provides the basis for the two rrst
widely used theoretical tests: the spectral test—2—
12]arid the lattice test [2, 14]. Briefly,
spectral test in n dimensions looksat(in a
ndardized reciprocal scale) the distance
tween adjacent hyperplanes in the mostwidely
ratedfamily of hyperplanes inthe lattice of
uplesproduced by the generator. The lattice
St looks at the length ratio of longest arid
)rtest sides in a reduced basis for thelattice
n-tuples.It is not surprising, then, that
se two tests are rather closely related. Still,
may help our understanding to apply both tests
study the results in dimensions 2 through 6.
it actualtest criteria should we use?in
ith'snotation [12] for the spectral test we
Lculate C for n2,3,4,5,6 andrequirethatall sevalues be at least1 .Thisis the more
ingent of the two criteria Knuth suggested,but
entempirical experience [9] suggests that it
not unreasonably difficult to find multipliers
Lch meet this requirement. For the lattice test
us use I to denote the length ratio of long-
side to shortest side. Marsaglia suggested
;] requiring L 2, and this seems sensible,
un for n:2,3,4,5,6. (Foniulated in such simple
ns as "Cs >1"and ??I <2",the lattice test
iears to be the more stringent of the two). to
Lnts surrunarize this discussion of theoretical
)perties: these theoretical tests are much
:ter for screening congruential generators than
}cuown empirical tests, and no congruential
ierator should be put into use without passing
!Tfl.
3. A BASIC LIBRARY
:'s tum now to what we can do with a carefully
)sen uniform generator. What other generators
)uld we put with it to form a serviceable set
'most simulation purposes? We should not have
h difficulty agreeing on a basic library, and
'easonable list should look raich like this:
Continuous distributions
Uniform (0,1)
Gaussian (0,1)
Exponential
Gamma(and x2)
Beta (and F)
Student's t
Discrete distributions
Uniform
Binomial
Poisson
iumber of other distributions may suggest them-
.ves for specific applications, but any reason-
.e library should support those in thebasic
;t.
eachnon-uniform distribution we of course tto use exact and efficient algorithms and
id such methods as inverse c .d .f.approximations
and the Central Limit Theorem. Clever exact
methods for the Gaussian and the exponential have
been available for about ten years, but until
recently the situation for the general case in
most of the other distribution families was not so
encouraging. Fortunately a number of new algo-
rithmshaveappeared during the last year or so to
improve matters. Thework of Dieter and Ahrens
(for example,[1,6]) is particularlynoteworthy
here;their acceptance-rejection methods for the
gamnadistributionand the beta distribution [6]
havethe attractive property ofrequiring nearly
constanttime regardless of the parameter(s) of
thedistribution.
4 I1LEMENTATION
Nowhow shouldwe go about implementingour chosen
random-numbergenerators? Therangeof issues
here is quite broad --fromchoosing the level of
programming language to being careful, in an
assembly-language urd.fonrt generator, not to throw
away significant bits when converting the result
to floating-point. Let's look at some of the
questions from the top down.
1. What should the generator return? For most
applications the convenient form of output is a
vector of random numbers, and this means we will
be producing subroutines instead of functions. In
some cases a function might be better, but we
would want to balance this against the overhead
(both in progranining effort and storage space) of
adding the function foth to the library.
2. How should calling sequences (parameter lists)
bestructured? To integrate the routines as a
library, we would put the comion parameters first,
as in
RANDOM (X, N, [other parameters]).
3. How should we organize the way in which non—
uniform generators use a basic uniform generator
or generators? Here it is likely to be cleaner if
each routine which requires a source of uniform
randoii numbers actually incorporates its own.
This would consume little space and eliminate con-
siderable subroutine linkage especially in the
nore complicated rejection algorithms. By remov-
ing "side-effect" interactions among different non—
unifonn generatcrs it should make complex simula-
tion programs easier to debug.
4. How should we handle starting values? In
order to reproduce results the user must be able
to set the starting value(s) (usually for the
basic uniform genenator( s)) and recover the
current value at any point in the sequence. For
the user who wants a "random" start we can provide
a routine which uses the system clock or some
other such source.
5. In what language should we proguam the
tors? This question may receive more varied
ers than the previous ones. Many generators,
especially uniform ones, have in the past been
coded in assembly language because the result runs
.
.faster and because most higher-level languages
,
don'thave the primitives for the operations in-
volved. Now, however, it seems preferable to use
higher-level languages (such as FORTRAN or FL/I)
as much as possible. For one thing, this is the
only sensible way to approach portability from one
line of computers to another, and having machine-
independent generators will facilitate replication
of simulation studies, something we have largely
neglected. Another important consideration is the
correctnessof the implementation: assembly-
language generators arelikelyto have more bugs,
and those bugs will be hardertoisolate. One
IBM/360 assembler implementation of Marsaglia's
rectangle-wedge-tail algorithm for the Gaussian
distribution[L]providesa good example. Because
theprogrammer misused one of the machine instruc-
tions, the generator produced an excessive number
of deviates with large magnitude (like 5 and 6!).
It's reasonable to admit that one can gain a good
deal of speed in most random-number generators by
codingthemin assembly language, but the conclu-
sion has to be that we should never start at that
level. Program the generator in a higher-level
languageanddebug it thoroughly so that there
will be a well-understood version to compare the
assenbly-language one against.
6. Whatshouldwe do about testing? The simple
answer,ofcourse, is "Be thorough". Thisisold
advice, but manygeneratorsseem not to get a very
extensivewor1ut. Forexample, apparently the
only testapplied to the Gaussiangenerator [4]
mentionedpreviously wasachi-squared test based
ondividing the reallineinto 20 intervals of
equalprobability content. Since each taillies
entirely within one of these intervals, there was
no check on the tail part of the algorithm. A
simple probability plot would have exposed the
problem almost immediately. This example suggests
a natural strategy: the testing should be
designed to cover each segment of a complicated
algorithm(inaddition to the performance of the
whole). This is valuable whenthe implementation
isin a higher-level language, andit isvital
whenassemblylanguage is involved.Testingalso
should reveal something about the comparative
speedofthe algorithm because this is often a
nre complicated question thantheoreticalcalcu-
lations (of such things as the average numberof
uniform deviates used inarejection algorithm)
can answer.Forexample, W.M.Gentleman told me
recently thatona Honeywell 6000-series computer
the Gaussian algorithm of Brent[3]runsabout35%
slower thanthe1964 algorithm of Marsaglia and
Bray[15]. InfonBation like this is machine-
dependent butstillquite useful.
5.DOCUMENTATION
Finally we cometo documentation--themost
important stepin making a generator orlibrary
accessibleto users. Here the procedure is
straightforward,but lapses are frequentenough
to demand a briefdiscussion. There are two basic
aspects: use of the generator andits "pedigree".
—3—
Documentationdescribinguse is whateveryprogram-
mer will readiimnediately,anditshould start with
a precise statement of what the generator is and
what it produces. (It may be that this goes with-
out saying, but an earlier (1 July 1973) edition of
the IMSL Library1 Manual[11] did nottellwhat
congruent ial generator was implementedin the sub-
routineGGU1;it was necessary to read the assembly
code [7].Fortunatelythis is nolonger truein
the latest edition.) Otheressential details for
usearethecalling sequence or parameter list,
restrictions onparameters(forexample, thestart-
ing value), what other generators are used, and
thedefault initialization.
Toestablish a generator's "pedigree", supporting.
documentation should report the specific algorithm
(with information onitsefficiency), relevant
theoretical properties (especially for uniform
generators, including any embeddedina non-uniform
one), the sourcesof anyprevious implementations
onwhich thepresent one is based, arid the results
of testing.Together,these should give the user
anadequately detailedpicture of the generator.
6. S1J?ARY
This paper has briefly endeavoredtogive an up-to-
dateconsumer 'sview of random-number generators.
Specific recormnendationscoveruniform generators,
thecompositionof a basic library,and principles
of implementationand documentation. While a
number of actual exarnpls indicate that currently
available generators and libraries often fall short
of the best that we 1.a-iow how to do, it is reassur-
ing tonotethat most of the tools needed for sub-
stantial improvement are ready to hand. We should
now expect (and perhaps demand) the gap between
possibility and prtice to close rapidly.
REFERENCES
[1] J. H. Ahrens and U. Dieter, "Computer Methods
for Sampling from Gariuna, Beta, Poisson, and
Binomial Distributions," Computing 12 (1974),
2232146.
[2] W. A. Beyer, R. B. Roof, and Dorothy
Williamson, "The lattice Structure of Multi-
plicativeCongruential Pseudo-Random Vectors,"
MathematicsofComputation 25 (April 1971),
345—363.
[3] Richard P. Brent,"Algorithm 488:A Gaussian
Pseudo-Random Number Generator,"
tionsof the ACM 17,12 (December 1974), 704-
706.——
['4] lovickEdward Cannon III, "Pseudo Random
Number Generators for Statistical Applica-
tions," Technical Report 69, Department of
Statistics and Computer Sciences, University
of Georgia, August 1971.
[5] R. R. Coveyou andR.D. MacPherson, "Fourier
Analysisof UniformRandom NumberGenerators,"Journal of the Association for Computing
Machiner1Tl967), 100—119.
[6] U. Dieter andJ. H.Ahrens,"Acceptance-
rejectionTechniques for Sampling from the
Ganvnaand Beta Distributions," Technical
Report 83, Department of Statistics, Stanford
University,May29, 1974.
[7] David C.Hoaglin, Remarks on the IMSL
Randan Number Generator GGtJ1," unpublished.
[8] David C. Hoaglin and David F. Andrews, "The
Reporting of Computation-based Result in
Statistics."In revision for The Mierican
Statistician.
[9] David C. Hoaglin andCordon Sande,"A Study
of MultipliersforPseudo-Random Number
Generatorswith Modulus 231_i."Presentedat
Joint Statistical Meetings, St. Louis,
Missouri,August 1974.
[10]IBM Corporation, System/360 Scientific Sub-
routine Package (360A—CM-03X) Version III,
Programmer'sManual. H20—0205—3,1968.
[11] InternationalMathematical and Statistical
Libraries,Inc.,The IMSL Librax1Reference
Manual, Edition (TTRANIV,S/370-
360).
[12] D. E. Knuth,TheArt of ComputerProgramming,
Volume 2: Seumerical Algorithms. Addison—
Wesley,1969.
[13] C. Marsaglia, "Regularities in Congruential
RandomNumberGenerators,"Numerische
matik 16 (1970), 8—10.
[14] G. Marsaglia, "TheStructure ofLinear Congru-
entialSequences,"ApplicationsofNumber
Theory toNumerical Analysis (S. K. Zaremba,
editor)249—285.Academic Press, 1973.
[15] G. Marsaglia andT.A. Bray,"AConvenient
Method for Generating Normal Variables,"
SIANReview6, 3 (July 19614), 260—264.
S

Bray,&quot;AConvenient Method for Generating Normal Variables,&quot;

Brent,&quot;Algorithm 488:A Gaussian Pseudo-Random Number Generator,&quot;

Computer Methods for Sampling from

Fourier Analysisof UniformRandom NumberGenerators,&quot;Journal of the Association for Computing

H.Ahrens,&quot;AcceptancerejectionTechniques for Sampling from the Ganvnaand Beta Distributions,&quot;

Hoaglin andCordon Sande,&quot;A

Knuth,TheArt of ComputerProgramming, Volume 2: Seumerical Algorithms.

Pseudo Random Number Generators for Statistical Applications,&quot;

Regularities

Remarks on the IMSL Randan Number Generator GGtJ1,&quot; unpublished.

The lattice Structure of MultiplicativeCongruential Pseudo-Random Vectors,&quot;

The Reporting of Computation-based Result in Statistics.&quot;In revision for The Mierican Statistician.

TheStructure ofLinear CongruentialSequences,&quot;ApplicationsofNumber Theory toNumerical Analysis

Implementing and Documenting Random Number Generators

Implementing and Documenting Random Number Generators

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics