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.

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 

Uniform uniform distribution
Constant return a constant
Exponential negative exponential distribution
Cauchy Cauchy distribution
Normal normal distribution
ChiSq non-central chi-squared distribution
Gamma gamma distribution
Pareto Pareto distribution
Poisson Poisson distribution
Binomial binomial distribution
NegativeBinomial negative binomial distribution

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

PosGenX Positive random numbers with a decreasing density
SymGenX Random numbers from a symmetric unimodal density
AsymGenX Random numbers from an asymmetric unimodal density
PosGen Positive random numbers with a decreasing density
SymGen Random numbers from a symmetric unimodal density
AsymGen Random numbers from an asymmetric unimodal density
DiscreteGen Random numbers from a discrete distribution
SumRandom Sum and/or product of random numbers
MixedRandom Mixture of random numbers

Each of these classes has the following member functions  

Real Next() Get a new random number
char* Name() Name of the distribution
ExtReal Mean() Mean of the distribution
ExtReal Variance() Variance of the distribution

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.
 

RandomPermutation Draw numbers without replacement
RandomCombination Draw numbers without replacement and sort

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.
 

This version of newran does not do memory clean-up with the simulated exceptions.

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)

Borland 5.01A OK
Borland 3.1 See notes on 16 bit under testing
Microsoft 5.0 OK
Watcom 10A OK
Sun CC 4.2 OK
Gnu G++ 2.7.2 (Linux), 2.8.0 (Sun) OK

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).
 

If you are compiling on a PC with a 16 bit compiler you will need to set n_large to be 8000 rather than 1000000 and n to 8000 rather than 200000 in tryrand.cpp. This will change the output from the test program from that given in tryrand.txt although the general appearance should be the same.

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";
Note that the probability density pdf must drop to exactly 0 for the argument large enough. For example, include a statement in the program for pdf that, if the value is less than 1.0E-15, then return 0.

SymGenX:

This corresponds to PosGenX for symmetric distributions.
 

Note that the probability density pdf must drop to exactly 0 for the argument large enough. For example, include a statement in the program for pdf that, if the value is less than 1.0E-15, then return 0.

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";
Note that the probability density pdf must drop to exactly 0 for the argument large (large positive and large negative) enough. For example, include a statement in the program for pdf that, if the value is less than 1.0E-15, then return 0.

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.
 

Note that the probability density p(x) must drop to exactly 0 for the argument, x, large enough. For example, include a statement in the program for p(x) that, if the value is less than 1.0E-15, then return 0.

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:

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:

Distribution type Method Example
Continuous finite unimodal density (no parameters, can calculate density) Use PosGenX, SymGenX or AsymGenX.
Continuous finite unimodal density (with parameters, can calculate density) Derive a new class from PosGen, SymGen or AsymGen, over-ride Density. Gamma2
Can calculate inverse of distribution Transform uniform random number. Pareto
Transformation of supported random number Derive a new class from the existing class ChiSq1
Transformation of several random numbers Derive new class from Random; generate the new random number from the existing generators. ChiSq
Density with infinite singularity Transform a random variable generated by PosGen, SymGen or AsymGen. Gamma1
Distribution with several modes Breakdown into a mixture of unimodal distributions.
Linear or quadratic combination of supported random numbers Use SumRandom.
Mixture of supported random numbers Use MixedRandom.
Discrete distribution (< 100 possible values) Use DiscreteGen. Poisson2
Discrete distribution (many possible values) Use PosGen, SymGen or AsymGen. Poisson1

 

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:

readme.txt readme file
newran.htm this file
newran.h header file for newran
newran.cpp main code file
extreal.h header file for extended reals
extreal.cpp code file for extended reals
boolean.h definition of bool
include.h option file
tryrand.cpp test file
tryrand1.cpp called by tryrand - histograms of simple examples
tryrand2.cpp called by tryrand - histograms of advanced examples
tryrand3.cpp called by tryrand - statistical tests
tryrand4.cpp called by tryrand - test permutations
hist.cpp called by tryrand - draw histogram
tryrand.txt output from tryrand
nrcc.mak make file for CC compiler
nrgnu.mak make file for gnu G++ compiler
nrwat.mak makefile for Watcom 32 bit 10a compiler

 

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:

 

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