Random Number Generation for Discrete Variables

The authors present a surprisingly efficient random number generator that applies in principle to any discrete univariate probability mass function without alteration.

by Colin Rose and Murray D. Smith

In this article we present a generator for random sampling from discrete distributions that is both general and efficient: general because it can in principle be implemented for any discrete distribution, and large-sample efficient relative even to distribution-specific solutions such as those implemented in the Statistics Packages.

Introducing RNG [ ]

Let univariate X have discrete domain W = {(Subscript[x, 0]), (Subscript[x, 1]), ... }, with known probability mass function [pmf]  f(x) = P(X = x), and cumulative distribution function [cdf]  F(x) = P(X ≤ x), such that (∑_(x ∈ Ω))f(x)  =  1. Consider the problem of generating pseudo-random copies of X. One well known approach is the inverse method:  if u is a pseudo-random drawing from U(0,1), the uniform distribution defined on the unit interval, then x = (F^(-1))(u) is a pseudo-random copy of X. Of course, this method is only desirable if the inverse function of F is computationally tractable. Unfortunately, many discrete distributions do not satisfy this condition, so alternative approaches are often required; see [Peterson and Kronmal 1983] and [Stuart and Ord 1994, chapter 9]. In this paper we present a discrete pseudo-random number generator entitled RNG[ ] that is virtuous in two respects. First, it is universal Ð it applies in principle to any discrete univariate pmf without alteration. This is achieved by constructing (F^( -1))(u) as a lookup table, instead of trying to do so symbolically. Second, this approach is surprisingly efficient. Given that pluralism and efficiency are usually mutually incompatible, the attainment of both goals is particularly pleasing.

The package function RNG[n, f] generates n pseudo-random copies of a discrete variable X, with pmf f. It allows f to take either a Function Form or a List Form.  Function Form inputs f(x) as a function of x. It requires that the elements of Ω are consecutive integers which increment by unity. Most common (and even not so common) discrete distributions satisfy this condition  Ð  see the Poisson, Binomial and Yule examples below. By contrast, List Form inputs f(x) as a list {f((Subscript[x, 0])), f((Subscript[x, 1])),..., f((Subscript[x, m]))}, for finite m  Ð  see the Frequency Table example below. List Form is especially useful if the elements of Ω are not consecutive integers.

Poisson

Let us suppose X ~ Poisson(6),  X ∈ W = {0, 1, 2, 3, ...} :

\(f = \(E\^\(-λ\)\ λ\^x\)\/\(x!\) /. λ -> 6; 
    
domain[f] = {x, 0, ∞}; \)

Note that the term domain[ ] must always be entered along with f, as it passes important information onto RNG[ ]. Here then are 10 000 pseudo-random copies of X:

aa = RNG[10000, f]; // Timing

We can contrast the empirical distribution of data set aa with the true distributon of X:

empfreq = Table[Count[aa,i],{i,Min[aa],Max[aa]}]/Length[aa]//N
truefreq = Table[f, {x, Min[aa], Max[aa]}] // N

To represent the above as a histogram, first load the package

<<Graphics`Graphics`

and then BarChart the data. Here, Tik[] denotes a package-defined function that modifies Tick marks on the horizontal axis:

BarChart[empfreq, truefreq, Tik[aa], PlotRange -> All];

The left column (red) denotes the empirical pmf, whilst the right column (blue) denotes the true density. One obtains a superb fit because RNG is an exact solution. Since the Poisson distribution is included in Mathematica's Statistics Package, we can use it to establish a performance benchmark.  After loading the master package,

<<Statistics`Master`

generate 10 000 copies from Poisson(6) using the customised routine contained in this package:

pdist = PoissonDistribution[6];
Table[Random[pdist], {10000}]; // Timing

For this example, RNG is over 250 times faster than the Statistics Package, even though RNG is a general solution that has not been specially optimised for the Poisson case. Of course, results vary by both distribution and parameter value, as illustrated by the tables below. An X Factor of 100 means that RNG is 100 times faster than the Statistics Package (i.e. takes 1/100th of the time). Results are based on a sample size of n = 10000, are rounded for presentation, and computed on a PowerMac 8500/120. Parameters are ordered as per the Statistics`DiscreteDistributions` package.

		       ({{Poisson (λ)}, {({{ λ, .5, 2, 20, 50 }, {X Factor, 100, 169, 314, 259 
                      }})}})	        ({{LogSeries(θ)}, {({{ θ, .01, .2, .5, .75 , .99}, {X Factor, 85, 95, 121, 160 , 153}})}})

   ({{Hypergeometric (m, r, 2 m)}, {({{ , r=m/5, r=m/3, r=m}, {m=15, 240X, 249X, 279X}, {m=30, 318X, 365X, 542X}, {m=60, 460X, 658X, 1441X  }})}})            ({{Binomial(m, p)}, {({{ , p=.1, p=.5, p=.99}, {m=1, 5X, 5X, 5X}, {m=10, 7X, 6X, 8X}, {m=50, 17X, 12X, 15X}})}})  	     ({{Negative Binomial(m, p)}, {({{ , p=.05, p=.5, p=.95}, {m=1, 2X, 4X, 4X}, {m=10, 4X, 5X, 5X}, {m=50, 7X, 9X, 13X}})}})
   
In every case examined here, the X Factor favours RNG.

Binomial

The Binomial distribution is often encountered in practice.  We generate from it using the Function Form approach. For example,

\(p = .3; m = 10; 

g = Binomial[m, x]\ p\^x\ \((1 - p)\)\^\(m - x\)\)domain[g] = {x, 0, m}

Here are 20 pseudo-random copies from this distribution:

RNG[20, g]

Yule

Because RNG[ ] is a general solution, it can generate random numbers from, in principle, any discrete distribution, and not just those distributions which appear in the Statistics Package.  Consider for instance the following Yule distribution:

h = 
      ((ϕ     ϕ!)/Pochhammer)[x, ϕ + 1] /. 
        ϕ -> 6domain[h] = {x, 1, ∞}

where, for this example, parameter ϕ ∈ {1, 2, 3, ...} has been set to 6.  It is easy to generate copies of X with RNG:

aa = RNG[10000, h]; // Timing

We again contrast the true and empirical distributions:

empfreq = Table[Count[aa,i],{i,Min[aa],Max[aa]}]/Length[aa]//N
truefreq = Table[h, {x, Min[aa], Max[aa]}] // N

BarChart[empfreq, truefreq, Tik[aa], PlotRange -> All];

A frequency table

The three previous examples have been conveniently tackled using what we term the Function Form of input.  In each case, the domain[ ] was defined over a subset of ℤ (the set of integers) with unit length between consecutive elements. Discrete variables may not always be defined in this way.  For example, suppose variable X is distributed as follows:     

             ({{f(x):, .4, .1, .2, .3}, {   x   :, -7,  3.4,  8.9,  13}})

We represent such examples using a "List Form":

f  = {0.4, 0.1, 0.2, 0.3};
domain[f] = { -7, 3.4, 8.9, 13 };

Here are 10 pseudo-random numbers from f:

RNG[10,  f]

And here are a few more ...

aa = RNG[10000, f]; // Timing

The empirical pmf and histogram are now given by:

empfreq = Map[Count[N[aa], #]&, domain[f]//N] / Length[aa] //N

BarChart[empfreq, f,  i=1; Ticks -> {Map[{i++, #}&, domain[f]], Automatic}];

How RNG works

To illustrate how RNG works, consider a Poisson density X ~ Poisson(1) defined on W = {0, 1, 2, ...}:

f = 
      (((E^(-λ)) (λ^x))/(x!)) /. λ -> 1domain[f] = {x, 0, ∞}

We proceed in 2 steps.  First, construct F(x) for every x ∈ {0, 1, 2, 3, ..., c}, where c is the smallest integer such that
1-F(c) < ε, for small ε > 0.  Here we set the tolerance level to ε = (10^(-6)):

i = 0; F = 
      FixedPointList[N[# + f /. x -> i++]&, 0, 
        SameTest -> (1 - #2 < (1/(10^6))&)]

Second, given F, create the inverse function (F^(-1))(u) in the form of a lookup table.  For instance,  if 0 < u < 0.367, let x = 0, and so on ...

j = 0;
Ergo[u_] = Which@@Flatten[Append[MapThread[{#1<u<#2, j++}&, 
              {Drop[F,-1], Rest[F]}], {True, j-1}]]

Then a Poisson(1) pseudo-random number is given by:

Ergo[Random[]]

The RNG function is actually a little more complex than this because it may also impose 'left censoring', it is then sorted by bin size, and is finally compiled. The above nevertheless conveys the basic concept.

On Quality

RNG adopts an exact methodology. Nevertheless, its quality can be effected by (i) censoring, and (ii) numerical precision. This section discusses these points in turn.

Censoring: A lookup table like Ergo is naturally limited to a finite number of elements. Thus, if variable X has infinite support (e.g., when X is a Poisson variable), one has to censor the tail(s) of the distribution.  For example, in step 1 above, the pmf was censored at a point c in such a way that less than ε = (10^(-6)) of the density mass was disturbed.  In some cases it may also be desirable to censor the left tail of the distribution.  Doing so yields smaller lookup tables and a more efficient generator. The RNG[n, f ] function automatically chooses upper and lower censoring points (say (OverBar[c]), (UnderBar[c]) ) by imposing default tolerance settings of (OverBar[ε]) = (UnderBar[ε]) = (10^(-6)).  The related package function RNGBounds[ f ] identifies the points (OverBar[c]) and (UnderBar[c]) , and calculates the probability mass that is theoretically effected by censoring.  For example, for our Poisson(1) example:

RNGBounds[f]

The printed output tells us that when RNG is used at its default settings, it can generate copies of X from the set \(TraditionalForm\`\(Ω\^*\)\) = {0, 1, ..., 9}.  By censoring at \(TraditionalForm\`c\&_\) = 9, outcomes \(TraditionalForm\`\(Ω\_*\)\)  = {10, 11, 12, ...} are reported as 9.  Thus, the censored mass is not lost; it is merely shifted to the censoring point. The density mass shifted in this way corresponds to f(10) + f(11) + f(12) + ... =

(∑_(x = 10)^∞) f  // N

as reported by RNGBounds[ ] above.  Since this is smaller than \(TraditionalForm\`ε\&_\), it satisfies the tolerance constraint.  Other methods of accounting for the censored mass include:  allocating it equally across all elements in \(TraditionalForm\`\(Ω\^*\)\) by scaling; truncation (proportionate scaling); or assigning it entirely to the mode of X.

At the default tolerance setting, the effect of censoring will usually be unobservable, unless one works with very large data sets (larger than (10^6)values).  For instance, if aa denotes a simulated set of random numbers, and Max[aa] < (OverBar[c]), then (OverBar[c]) has clearly not effected the random number generation process.  Let us try this out:

aa = RNG[10000, f];  Max[aa]

Here, Max[aa] < 9, so censoring could not have effected the data set.  In fact, if the simulation sample size was as large as n = (10^7), we would expect on average to observe only one copy of X lying in the censored region \(TraditionalForm\`\(Ω\_*\)\).

Setting your own e values
If censoring at X = 9 is not desirable (for instance if n = (10^8)), one might want to impose tighter tolerance levels. The function RNGBounds[ f, (UnderBar[ε]), (OverBar[ε]) ] can be used to inspect the effect of arbitrary tolerance settings, while RNG[n, f, (UnderBar[ε]), (OverBar[ε]) ] imposes those settings on the generator.  For example:

RNGBounds[f, (10^(-6)), (10^(-11))]

Thus, RNG[n, f, (10^(-6)), (10^(-11))] will generate n copies of X whose pmf is  f(x), with outcomes restricted to the integers in \(TraditionalForm\`\(Ω\^*\)\) = {0, 1, ..., 13};  censoring occurs at (OverBar[c]) = 13 which results in just over 4.5×(10^(-12)) of the density mass being shifted to that point.  Here are 20 pseudo-random numbers generated in this fashion:

RNG[20, f, (10^(-6)), (10^(-11))]

RNG and RNGBounds are defined for ε ⩾ (10^(-15)). Setting ε outside this interval is not meaningful (see below) and may cause problems. [ It is also assumed that ε < .25, though this constraint should never be binding.]  In List Form examples, the distribution is not censored, so RNGBounds does not apply and by design will not evaluate. For Function Form examples, we recommend that whenever RNG is applied, the printed output from RNGBounds should also be inspected.

Numerical precision: RNG uses machine-precision numbers in its calculations (typically 16-digit precision). Moreover, it makes use of Random[ ] which also has 16-digit precision. For both these reasons, it is not meaningful to set tolerance settings smaller than roughly (10^(-15)).  If precision greater than 16 digits is required, it can be obtained by manually re-configuring the RNG function to use arbitrary-precision numbers. To do so, you will need (i) to change N[expr] to N[expr, p] where p denotes the desired precision;  (ii) to remove the Compiling function as the latter uses machine-precision mathematics, and (iii) change Random[ ] to Random[Real, {0, 1}, \(TraditionalForm\`\(p\^+\)\)] where \(TraditionalForm\`\(p\^+\)\) > p  [see The Mathematica Journal, 6(3), p.17]. Alternatively, [Berdnikov, Turtia and Compagner 1996]  provide a high-quality (0,1) random-number generator that can be used as an input stream instead of Random[ ] in a manual setup. The attainment of extra numerical precision is of course not without cost: moving from machine-precision to arbitrary-precision numbers is computationally expensive.

Other caveats:  By constructing a lookup table, RNG trades off a small fixed cost in return for a lower marginal cost. This tradeoff will be particularly beneficial if a large number of copies are required. If only a few copies are needed, it may not be worthwhile. The fixed cost itself is proportional to the size of the lookup table. For instance, a discrete uniform distribution such as f = (10^(-6)) defined on Ω = {1, 2, 3, ..., (10^6)} will require constructing a huge lookup table with a million elements. Here, a technique such as  Random[Integer, {1,(10^6)}] is clearly more appropriate.

Acknowledgements
We thank Robert Bartels for helpful suggestions.

References
Berdnikov, A.S., Turtia, S.B., and A. Compagner (1996), A MathLink Program for High-Quality Random Numbers, The Mathematica Journal, 6(3), Summer 1996, 65-69. 

Peterson, A. V. and Kronmal, R. K.,  (1983), Analytic Comparison of Three General Purpose Methods for the Computer Generation of Discrete Random Variables, Applied Statistics, 12(3), 276-286.   

Stuart, Alan and Ord, K. (1994), Kendall's Advanced Theory of Statistics (volume 1), sixth edition, Edward Arnold.

About the Authors

Theoretical Research Instituten
  Bellevue Hill
  Sydney,  NSW 2023
  colinr@extro.ucc.su.oz.au)

 Dr Murray D. Smith   
 Department of Econometrics
 University of Sydneyn
Sydney,  Australia    2006
murrays@sue.econ.su.oz.au