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.

Comments and bug reports to robert at statsresearch.co.nz [replace
at by you-know-what].

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.

There are three classes for generating random numbers where we want to
vary the parameters at each call. It would be inefficient to use the previous
classes since you would need to set up a random number object at each call.

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, 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 (very old) compiler does not recognises bool deactivate the
statement #define bool_LIB. This will activate my Boolean class.

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.0, 5.5, 5.6

OK

Microsoft 5.0, 6.0, 7.0

OK

Watcom 10A

OK

Sun CC 6

OK

Gnu G++ 2.96 (Linux), 2.95 (Sun)

OK

Intel 5 for Windows or Linux

OK

I have included make files for Borland 5.5, CC and Gnu G++. See files
section. You can use my genmake program to generate
make files for other compilers.

Testing

The files tryrand.cpp, tryrand1.cpp, tryrand2.cpp, tryrand3.cpp,
tryrand4.cpp, tryrand5.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 results from Borland 5.5 are in tryrand.txt. Other compilers may
give different but still correct results. In particular, Microsoft C++ compilers
give different results when global optimisation is set. This is most likely due
to different round-off error but may also be due to different order of
evaluation of expressions. The appearance of the histograms in the output should
be similar to that in tryrand.txt and the statistical tests should still be
passed. There are some notes in tryrand.txt for assessing the output.

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. It is very dubious if you are going to call it more than 100 million
times in a simulation. 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.

Remember that you need to call Random::Set(double seed) at the beginning of
your program. See the overview and tryrand.cpp.

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

Binomial:

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.

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.

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

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.

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.

VariPoisson

Use this class if you want to generate a Poisson random variable but you
want to change the parameter frequently, so using the
Poisson class would be inefficient. There is one member function

int VariPoisson::iNext(Real mu);

which returns a new Poisson random number with mean mu. To generate
100 Poisson random numbers with means 1,2,...,100 use the following program

The constructor is slow so put it outside any loop. The individual
calls to iNext should be quite fast. The method is approximate for mu >= 300. The
constructor is not in any class hierarchy and iNext is not virtual. This class is somewhat beta-ish and may change in a future release of
newran.

VariBinomial

Use this class if you want to generate a Binomial random variable but you
want to change the parameters of the frequently, so using the
Binomial class would be inefficient. There is one member function

int VariBinomial::iNext(int n, Real p);

which returns a new Binomial random number with number of trials n and
probability of success p. To generate
100 Binomial random numbers with n = 1,2,...,100 and p = 0.5 use the following program

The constructor is slow so put it outside any loop. The individual
calls to iNext should be quite fast. The method is approximate if both
n*p > 200 and n*(1-p) > 200. The
constructor is not in any class hierarchy and iNext is not virtual. This class is somewhat beta-ish and may change in a future release of
newran.

VariLogNormal

Use this class if you want to generate a log normal random variable and you
want to change the parameters of the frequently. There is one member function

Real VariLogNormal::Next(Real mean, Real sd);

which returns a new log normal random number with mean mean and
standard deviation sd. Note that mean and
sd are the mean and standard deviation of the log normal distribution and not of
the underlying normal distribution. To generate
100 log normal random numbers with mean = 1,2,...,100 and sd = 1.0 use the following program

The
constructor is not in any class hierarchy and Next is not virtual. This class is somewhat beta-ish and may change in a future release of
newran.

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.

Other people's code:

The gamma function is adapted from Numerical Recipes in C
by Press, Flannery, Teukolsky, Vetterling published by the Cambridge University Press.
The Shell sort and quick sort are adapted from Algorithms in C++ by
Sedgewick published by Addison Wesley.