Newran02A - a random number generator library

30 August, 1998

Copyright (C) 1989, 1995, 1998: R B Davies

Permission is granted to use or distribute but not to sell.

• Overview
• Getting started
• Descriptions of the classes to be accessed by the user
• Descriptions of the supporting classes
• Generating numbers from other distributions
• Other people's code
• Files included in this package
• Class structure
• To do
• History
• To online documentation page

This is a C++ library for generating sequences of random numbers from a wide variety of distributions. It is particularly appropriate for the situation where one requires sequences of identically distributed random numbers since the set up time for each type of distribution is relatively long but it is efficient when generating each new random number. The library includes classes for generating random numbers from a number of distributions and is easily extended to be able to generate random numbers from almost any of the standard distributions.

Newran02, fixes some bugs, includes more extensive testing, has been updated to be compatible with newmat09 and includes some additional distributions. Newran02A includes an updated version of the exception library.

Comments and bug reports to mailto:robertd@netlink.co.nz.

For updates and notes see http://webnz.com/robert/.

Overview

The following are the classes for generating random numbers from particular distributions 

The following classes are available to the user for generating numbers from other distributions  

Each of these classes has the following member functions  

These 4 functions are declared virtual so it is easy to write simulation programs that can be run with different distributions.

Real is typedefed to be either float or double. See customising. Note that Next() always returns a Real even for discrete distributions.

ExtReal is a class which is, in effect either a Real or one of the following: PlusInfinity, MinusInfinity, Indefinite or Missing. I use ExtReal so that I can return infinite or indefinite values for the mean or variance of a distribution.

There are two static functions in the class Random.

   void Random::Set(double)

must be called with an argument between 0 and 1 to set up the base random number generator before Next() is called in any class.

   double Random::Get()

returns the current value of the seed.

There are two classes for doing combinations and permutations.
 

Further details of all these classes including the constructors are given below.

Getting started

Customising

The file include.h sets a variety of options including several compiler dependent options. You may need to edit include.h to get the options you require. If you are using a compiler different from one I have worked with you may have to set up a new section in include.h appropriate for your compiler.

Borland, Turbo, Gnu, Microsoft and Watcom are recognised automatically. If none of these are recognised a default set of options is used. These are fine for AT&T, HPUX and Sun C++. If you using a compiler I don't know about, you may have to write a new set of options.

There is an option in include.h for selecting whether you use compiler supported exceptions, simulated exceptions, or disable exceptions. Use the option for compiler supported exceptions if and only if you have set the option on your compiler to recognise exceptions. Disabling exceptions sometimes helps with compilers that are incompatible with my exception simulation scheme.
 

If your compiler recognises bool as required by the standard activate the statement #define bool_LIB. This will deactivate my Boolean class. (I have included a command in the make file to do this automatically with Gnu G++).

Activate the appropriate statement to make the element type Real to mean float or double.

Activate the namespace option if your want to use namespaces and have a compiler that really does support them.

Compiling

You will need to compile newran.cpp, myexcept.cpp and extreal.cpp and link the resulting object files to your programs. Your source files which access newran will need to have newran.h as an include file.

Compilers

I have tested newran02 with the following compilers (all PC ones in 32 bit console mode)

I have included make files for Watcom 10A, CC and Gnu G++. See files section.

Testing

The files tryrand.cpp, tryrand1.cpp, tryrand2.cpp, tryrand3.cpp, tryrand4.cpp, hist.cpp run the generators in the library and print histograms of the resulting distributions. Sample means and variances are also calculated and can be compared with the population values. The correct results are in tryrand.txt (although there might be slight differences in the format of the output from different compilers).
 

The test program tryrand.cpp includes a simple test for memory leaks. This is valid for only some compilers. It seems to work for Borland C++ in console mode but not for Gnu G++ or Microsoft C++, where it almost always (presumably incorrectly) suggests an error.

Descriptions of the classes to be accessed by the user:

Random:

This is the basic uniform random number generator, used to drive all the others. The Lewis-Goodman-Miller algorithm is used with Marsaglia mixing. While not perfect, and now superseded, the LGM generator has given me acceptable results in a wide variety of simulations. See Numerical Recipes in C by Press, Flannery, Teukolsky, Vetterling published by the Cambridge University Press for details. The LGM generator does pass the Marsaglia diehard tests when you include the mixing. (It doesn't pass without the mixing). Nevertheless it should be upgraded. Ideally the basic generator should be recoded in assembly language to give the maximum speed to all the generators in this package. You can access the numbers directly using Next() but I suggest you use class Uniform for uniform random numbers and reserve Random for setting the starting seed and as the base class for the random number generators.

Uniform:

Return a uniform random number from the range (0, 1). The constructor has no parameters. For example

   Uniform U;
   for (int i=0; i<100; i++) cout << U.Next() << "\n";

prints a column of 100 numbers drawn from a uniform distribution.

Constant:

This returns a constant. The constructor takes one Real parameter; the value of the constant to be returned. So

   Constant C(5.5);
   cout << C.Next() << "\n";

prints 5.5.

Exponential:

This generates random numbers with density exp(-x) for x>=0. The constructor takes no arguments.

   Exponential E;
   for (int i=0; i<100; i++) cout << E.Next() << "\n";

Cauchy:

Generates random numbers from a standard Cauchy distribution. The constructor takes no parameters.

   Cauchy C;
   for (int i=0; i<100; i++) cout << C.Next() << "\n";

Normal:

Generates standard normal random numbers. The constructor has no arguments. This class has been augmented to ensure only one copy of the arrays generated by the constructor exist at any given time. That is, if the constructor is called twice (before the destructor is called) only one copy of the arrays is generated.

   Normal Z;
   for (int i=0; i<100; i++) cout << Z.Next() << "\n";

ChiSq:

Non-Central chi-squared distribution. The method uses ChiSq1 to generate the non-central part and Gamma2 or Exponential to generate the central part. The constructor takes as arguments the number of degrees of freedom (>=1) and the non-centrality parameter (omit if zero).

   int df = 10; Real noncen = 2.0;
   ChiSq CS(df, noncen);
   for (int i=0; i<100; i++) cout << CS.Next() << "\n";

Gamma:

Gamma distribution. The constructor takes the shape parameter as argument. Uses Gamma1, Gamma2 or Exponential.

   Real shape = 0.75;
   Gamma G(shape);
   for (int i=0; i<100; i++) cout << G.Next() << "\n";

Pareto:

Pareto distribution. The constructor takes the shape parameter as argument. I follow the definition of Kotz and Johnson's Continuous univariate distributions 1, chapter 19, page 234, with k = 1. The generator uses a power transform of a uniform random number.

   Real shape = 0.75;
   Pareto P(shape);
   for (int i=0; i<100; i++) cout << P.Next() << "\n";

Poisson:

Poisson distribution: uses Poisson1 or Poisson2. Constructor takes the mean as its argument.

   Real mean = 5.0;
   Poisson P(mean);
   for (int i=0; i<100; i++) cout << (int)P.Next() << "\n";

Binomial:

Binomial distribution: uses Binomial1 or Binomial2. Constructor takes n and p as its arguments.

   int n = 50; Real p = 0.25;
   Binomial B(n, p);
   for (int i=0; i<100; i++) cout << (int)B.Next() << "\n";

NegativeBinomial:

Negative binomial distribution. Constructor takes N and P as its arguments. I use the notation of Kotz and Johnson's Discrete distributions. Some people use p = 1/(P+1) in place of the second parameter.

   Real N = 12.5; Real P = 3.0;
   NegativeBinomial NB(N, P);
   for (int i=0; i<100; i++) cout << (int)NB.Next() << "\n";

PosGenX:

This uses an arbitrary density satisfying the previous conditions to generate random numbers from that density. Suppose Real pdf(Real) is the density. Then use pdf as the argument of the constructor. For example

   PosGenX P(pdf);
   for (int i=0; i<100; i++) cout << P.Next() << "\n";

SymGenX:

This corresponds to PosGenX for symmetric distributions.
 

AsymGenX:

Corresponds to PosGenX. The arguments of the constructor are the name of the density function and the location of the mode.

   Real pdf(Real);
   Real mode;
   .....
   AsymGenX X(pdf, mode);
   for (int i=0; i<100; i++) cout << X.Next() << "\n";

PosGen:

PosGen is not used directly. It is used as a base class for generating a random number from an arbitrary probability density p(x). p(x) must be non-zero only for x>=0, be monotonically decreasing for x>=0, and be finite. For example, p(x) could be exp(-x) for x>=0.

The method is to cover the density in a set of rectangles of equal area as in the diagram (indicated by ---).

   |
   x
   |xx------
   |  xx    |
   |    xxx |
   |.......xxx---------
   |        | xxxx     |
   |        |     xxxx |
   |        |.........xxxxx------------
   |        |          |   xxxxx       |
   |        |          |        xxxxxx |
   |        |          |..............xxxxxx----------------------
   |        |          |               |    xxxxxxx               |
   |        |          |               |           xxxxxxx        |
   |        |          |               |                  xxxxxxxx|
   +===========================================================================

The numbers are generated by generating a pair of numbers uniformly distributed over these rectangles and then accepting the X coordinate as the next random number if the pair corresponds to a point below the density function. The acceptance can be done in two stages, the first being whether the number is below the dotted line. This means that the density function need be checked only very occasionally and on the average only just over 3 uniform random numbers are required for each of the random numbers produced by this generator.

See PosGenX or Exponential for the method of deriving a class to generate random numbers from a given distribution.
 

SymGen:

SymGen is a modification of PosGen for unimodal distributions symmetric about the origin, such as the standard normal.

AsymGen:

A general random number generator for unimodal distributions following the method used by PosGen. The constructor takes one argument: the location of the mode of the distribution.

DiscreteGen:

This is for generating random numbers taking just a finite number of values. There are two alternative forms of the constructor:

   DiscreteGen D(n,prob);
   DiscreteGen D(n,prob,val);

where n is an integer giving the number of values, prob a Real array of length n giving the probabilities and val a Real array of length n giving the set of values that are generated. If val is omitted the values are 0,1,...,n-1.

The method requires two uniform random numbers for each number it produces. This method is described by Kronmal and Peterson, American Statistician, 1979, Vol 33, No 4, pp214-218.

SumRandom:

This is for building a random number generator as a linear or multiplicative combination of existing random number generators. Suppose RV1, RV2, RV3, RV4 are random number generators defined with constructors given above and r1, r2, r0 are Reals and i1, i3 are integers.

Then the generator S defined by something like

   SumRandom S = RV1(i1)*r1 - RV2*r2 + RV3(i3)*RV4 + r0;

has the obvious meaning. RV1(i1) means that the sum of i1 independent values from RV1 should be used. Note that RV1*RV1 means the product of two independent numbers generated from RV1. Remember that SumRandom is slow if the number of terms or copies is large. I support the four arithmetic operators +, -, * and / but cannot calculate the means and variances if you divide by a random variable.

Use SumRandom to quickly set up simple combinations of the existing generators. But if the combination is going to be used extensively, then it is probably better to write a new class to do this.

Example: normal with mean = 10, standard deviation = 5:

   Normal N;
   SumRandom Z = 10 + 5 * N;
   for (int i=0; i<100; i++) cout << Z.Next() << "\n";

Example: F distribution with m and n degrees of freedom:

   int m, n;
   ... put values in m and n
   ChiSq Num(m); ChiSq Den(n);
   SumRandom F = (double)n/(double)m * Num / Den;
   for (int i=0; i<100; i++) cout << F.Next() << "\n";

MixedRandom:

This is for mixtures of distributions. Suppose rv1, rv2, rv3 are random number generators and p1, p2, p3 are Reals summing to 1. Then the generator M defined by

   MixedRandom M = rv1(p1) + rv2(p2) + rv3(p3);

produces a random number generator with selects its next random number from rv1 with probability p1, rv2 with probability p2, rv3 with probability p3.

Alternatively one can use the constructor

   MixedRandom M(n, prob, rv);

where n is the number of distributions in the mixture, prob the Real array of probabilities, rv an array of pointers to random variables.

Normal with outliers:

   Normal N; Cauchy C;
   MixedRandom Z = N(0.9) + C(0.1);
   for (int i=0; i<100; i++) cout << Z.Next() << "\n";

or:

   Normal N;
   MixedRandom Z = N(0.9) + (10*N)(0.1);
   for (int i=0; i<100; i++) cout << Z.Next() << "\n";

RandomPermutation:

To draw M numbers without replacement from start, start+1, ..., start+N-1 use

   RandomPermutation RP;
   RP.Next(N, M, p, start);

where p is an int* pointing to an array of length M or longer. Results are returned to that array.

   RP.Next(N, p, start);

assumes M = N. The parameter, start, has a default value of 0.

The method is rather inefficient if N is very large and M is much smaller.

RandomCombination:

To draw M numbers without replacement from start, start+1, ..., start+N-1 and then sort use

   RandomCombination RC;
   RC.Next(N, M, p, start);

where p is an int* pointing to an array of length M or longer. Results are returned to that array.

   RC.Next(N, p, start);

assumes M = N. The parameter, start, has a default value of 0.

The method is rather inefficient if N is large. A better approach for large N would be to generate the sorted combination directly. This would also provide a better way of doing permutations with large N, small M.

ExtReal

A class consisting of a Real and an enumeration, EXT_REAL_CODE, taking the following values:

• Finite
• PlusInfinity
• MinusInfinity
• Indefinite
• Missing

The arithmetic functions +, -, *, / are defined in the obvious ways, as is << for printing. The constructor can take either a Real or a value of EXT_REAL_CODE as an argument. If there is no argument the object is given the value Missing. Member function IsReal() returns true if the enumeration value is Finite and in this case value of the Real can be found with Value(). The enumeration value can be found with member function Code().

ExtReal is used at the type for values returned from the Mean and Variance member functions since these values may be infinite, indefinite or missing.

Descriptions of the supporting classes:

ChiSq1:

Non-central chi-squared with one degree of freedom. Used as part of ChiSq.

Gamma1:

This generates random numbers from a gamma distribution with shape parameter alpha < 1. Because the density is infinite at x = 0 a power transform is required. The constructor takes alpha as an argument.

Gamma2:

Gamma distribution for the shape parameter, alpha, greater than 1. The constructor takes alpha as the argument.

Poisson1:

Poisson distribution; derived from AsymGen. The constructor takes the mean as the argument. Used by Poisson for values of the mean greater than 10.

Poisson2:

Poisson distribution with mean less than or equal to 10. Uses DiscreteGen. Constructor takes the mean as its argument.

Binomial1:

Binomial distribution; derived from AsymGen. Used by Binomial for n >= 40. Constructor takes n and p as arguments.

Binomial2:

Binomial distribution with n < 40. Uses DiscreteGen. Constructor takes n and p as arguments.

AddedRandom, SubtractedRandom, MultipliedRandom, ShiftedRandom, ReverseShiftedRandom, ScaledRandom, RepeatedRandom, SelectedRandom, AddedSelectedRandom:

These are used by SumRandom and MixedRandom.

Generating numbers from other distributions:

Other people's code:

The gamma function and Shell sort routines are adapted from Numerical Recipes in C by Press, Flannery, Teukolsky, Vetterling published by the Cambridge University Press.

Files included in this package:

Class structure:

The following diagram gives the class hierarchy of the package.

ExtReal.......................... Extended real numbers

Random........................... Uniform random number generator
 |
 +---Constant.................... Return a constant
 |
 +---PosGen...................... Used by PosGenX etc
 |    |
 |    +---PosGenX................ Positive random #s from decreasing density
 |    |
 |    +---Exponential............ Negative exponential rng
 |    |
 |    +---Gamma1................. Used by Gamma (shape parameter < 1)
 |    |
 |    +---SymGen................. Used by SymGenX etc
 |         |
 |         +---SymGenX........... Random numbers from symmetric density
 |         |
 |         +---Cauchy............ Cauchy random number generator
 |         |
 |         +---Normal............ Standard normal random number generator
 |              |
 |              +---ChiSq1....... Used by ChiSq (one df)
 | 
 +---AsymGen..................... Used by AsymGenX etc
 |    |
 |    +---AsymGenX............... Random numbers from asymmetric density
 |    |
 |    +---Poisson1............... Used by Poisson (mean > 8)
 |    |
 |    +---Binomial1.............. Used by Binomial (n >= 40)
 |    |
 |    +---NegativeBinomial....... Negative binomial random number generator
 |    |
 |    +---Gamma2................. Used by Gamma (shape parameter > 1)
 |
 +---ChiSq....................... Non-central chi-squared rng
 |
 +---Gamma....................... Gamma random number generator
 |
 +---Pareto...................... Pareto random number generator
 |
 +---DiscreteGen................. Discrete random number generator
 |
 +---Poisson2.................... Used by Poisson (mean <= 8)
 |
 +---Binomial2................... Used by Binomial (n < 40)
 |
 +---Poisson..................... Poisson random number generator
 |
 +---Binomial.................... Binomial random number generator
 |
 +---SumRandom................... Sum of random numbers
 |
 +---MixedRandom................. Mixture of random numbers
 |
 +---MultipliedRandom............ Used by SumRandom
 |    |
 |    +---AddedRandom............ Used by SumRandom
 |    |
 |    +---SubtractedRandom....... Used by SumRandom
 |
 +---ShiftedRandom............... Used by SumRandom
 |    |
 |    +---ReverseShiftedRandom... Used by SumRandom
 |    |
 |    +---ScaledRandom........... Used by SumRandom
 |
 +---NegatedRandom.......... .... Used by SumRandom
 |
 +---RepeatedRandom.............. Used by SumRandom
 |
 +---AddedSelectedRandom......... Used by MixedRandom
 |
 +---SelectedRandom.............. Used by MixedRandom

RandomPermutation................ Random permutation
 |
 +---RandomCombination........... Sorted random permutation

To do:

• More modern alternatives to the LGM generator;
• Additional generator classes;
• Better methods for combinations and permutations with large N and small M;
• An alternative (slower) class for Poisson with short setup time.

History:

August, 1998 - update exception package; work around problem with MS VC++ 5

January, 1998 - version compatible with newmat09

1995 - newran version, additional distributions

1989 - initial version

Go to top

To online documentation page