crng

Random-number generators (RNGs) implemented as Python extension types coded in C.

version 1.2 (released 23 Oct 2002)

Per Kraulis


Copyright notice

Copyright © 2000-2002 Per Kraulis

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License (file gpl.txt) for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

The MT19937 code has been copyrighted by Matsumoto & Nishimura, (see the copyright notice in that section), and is governed by the GNU Lesser General Public License (LGPL, file lgpl.txt, previously called 'GNU Library General Public License').


Official website


Introduction

The Python module crng implements random-number generators (RNGs) based on several different algorithms producing uniform deviates in the open interval (0,1), i.e. exclusive of the end-point values 0 and 1. A few continuous and integer-valued non-uniform deviates are also available. Each RNG algorithm is implemented as a separate Python extension type. The RNG types are independent of each other, but have very similar interfaces. The entire module is implemented as one single C source code file.

There are already several RNGs available for Python, notably the standard module random (implementing various random-number distributions) based on whrandom (the Wichmann-Hill algorithm). These are coded in Python, and hence slow. The purpose of the crng module is to provide efficient implementations coded in C of several modern RNG algorithms. The design allows easy comparison and interchange of the RNGs in Python code.

The crng module contains a C interface that allows other Python extension modules to directly call the C function for each RNG implementation. One may thereby attain execution speeds close to those possible for a direct implementation in C, while retaining the flexibility of the Python interface.

Several recent and reasonably well-tested RNGs are available in this module. This author is by no means an expert on the subject, but having read some of the recent literature, I would recommend using MT19937; it appears to combine speed of execution with good results in several standard tests of RNG quality. It has an astronomical period of approximately 4.3*106001.


Installation

The source code for the crng module consists of one and only one C code file crng.c, which is made available to the Python system by compiling it into a shared-object file or DLL, and placing it in an appropriate directory in the Python search path. The script file crng_pickle.py should be installed in the same directory, if it is needed.

The module has been tested with Python 1.5.2, and reportedly works with Python 2.0 also. Modifications have been made to the C source code to allow direct compilation using common Windows compilers.

Shared-object files for a few specific systems are available for direct download from the official web site. Otherwise, you must compile it yourself from the source code distribution.

Unpacking the distribution

The distribution is a gzipped tar file. To unpack it, do:

    gunzip crng-1.2.tar.gz
    tar xvf crng-1.2.tar
    cd crng-1.2

The files can also be downloaded one by one from the official web site.

Using setup.py with Distutils

If you have Distutils installed, you can use the setup.py script to perform the installation. Refer to the Distutils documentation.

Using old-style installation on Unix

If you do not have Distutils installed, you must use the old-style installation procedure. The following information is for Unix (including Linux) systems. The author has no experience with software development on Windows systems, so there is currently no information for Windows.

Fetch a copy of the file Makefile.pre.in to the current directory. Common locations for this file are the directories /usr/lib/python2.1/config/ or /usr/local/lib/python2.1/config/. Check the value of the sys.path variable (see below) for hints, or ask your local Python administrator. Do the following, using the appropriate directory:

    cp /usr/lib/python2.1/config/Makefile.pre.in .

Next, create the Makefile from the Setup.in file, and compile the module.

    cp Setup.in Setup
    make -f Makefile.pre.in boot
    make

This will create a shared-object file crng.so (at least on the Unix systems it has been tested on). This file must be copied or moved to a directory in the Python search path. To display your search path, do the following:

    python -c "import sys; print sys.path"

A common location is the /usr/lib/python2.1/site-packages directory, but this may require administrator privileges. Do the following, using the appropriate directory:

    cp crng.so /usr/lib/python2.1/site-packages

Also copy the crng_pickle.py script to the same directory.

You should now be able to run the example scripts that are part of the distribution:

    python validate.py
    python profile.py
    python display.py

Usage

All RNGs in the crng module are used in basically the same way:

  1. create an instance of the chosen RNG type, giving the appropriate seeds and initialization parameters,
  2. call that instance repeatedly to obtain new random numbers.

    >>> import crng
    >>> r=crng.ParkMiller(seed=123)
    >>> r()
    0.000962643418909
    >>> r(2)
    (0.179147941609, 0.939454629989)

The initialization call at object creation is specific to each RNG type, and is described in the sections for each RNG below.

The deviates of distributions other than the uniform (0,1) distribution are generated by an appropriate mathematical transformation of the values from a uniform (0,1) distribution. Therefore, the non-standard deviate objects must all be initialized with a basic RNG as argument. By default, ParkMiller is used, but any basic RNG from the crng module may be specified.

    >>> n=crng.NormalDeviate(rng=crng.Ranlux(), stdev=4.0)
    >>> n(4)
    (0.960612428619, 6.30989834512, -0.972826635646, -2.09557141754)

All RNGs implement the same methods, except for the initialization call. The members of the basic RNGs (if any) are usually read-only. The members of the deviates other than the uniform (0,1) distribution are usually read/write.

A few utility functions for sampling and shuffling are provided, which take as arguments a basic RNG and a sequence, i.e. a string, tuple or list. The result is a new sequence of the same type as the input.


Methods available for all RNG objects

Method call Returned value Comment
r.next(n=1) if n==1: single float
if n>1: tuple of floats
where float 0.0<f<1.0
The returned floating point value(s) are the next random number(s) in the series.
r(n=1) if n==1: single float
if n>1: tuple of floats
where float 0.0<f<1.0
Object instance call, which is equivalent to calling the method next.
r.compute(n=1) None Compute the next n random number in the series, where n is a positive integer, but do not return any values. This may be useful for profiling the C implementation of the RNG.
r.density(x) single float The density of the distribution function at x.
str(r) string
e.g. '<crng.ParkMiller object at 0x80a6860>'
The built-in str() function calls an internal function defined for the object. It returns an informal string representation of the RNG instance, which contains the type of the RNG and its memory address.
repr(r) string
e.g. 'crng.ParkMiller(seed=314159265)'
The built-in repr() function calls an internal function defined for the object. It returns a string representation of the RNG instance that is a valid Python code string, which can be used in an eval expression.


Saving and restoring an RNG object

The current state of any instance of a crng type can be retrieved through its members. This allows creation of an identical copy of an instance. One may also save an instance and retrieve it at a later stage without loss of information.

Basically, there are two independent ways to save and retrieve an instance of an crng type.

repr and eval()

All RNG types (basic as well as non-standard deviates) define an internal function called by the built-in function repr, which produces a string representation of the instance that is valid Python code. The string contains the module name prepended to the type name. This can be used to create copies of an RNG instance:

    >>> import crng
    >>> r1 = crng.ParkMiller(seed=123)
    >>> s = repr(r1)
    >>> print s
    crng.ParkMiller(seed=123)
    >>> r2 = eval(s)
    >>> print r1(), r2()
    0.000962643418461 0.000962643418461

Note that the string returned by the repr function is quite long for the RNGs that contain large state vectors (mainly Ranlux and MT19937). The non-standard deviates produce strings that also contain a representation of the basic RNG they use.

In Python 1.5.2, the string created by repr for a float value is sometimes not an exact representation of the value. This problem carries over to the representation of the RNG types. This has been fixed in Python 1.6 and later.

pickle/cPickle

The standard Python module pickle (and its optimized C implementation cPickle) can be used to store an RNG type in a safe manner. The script crng_pickle.py must be imported by a script that does a save and/or retrieve operation; this script performs the required setup as it is imported.

Here are two example scripts showing how to use pickle/cPickle with the crng types.

The first script saves an crng instance to a file using pickle. Note that the script crng_pickle must be imported by this script.

    import crng
    import crng_pickle
    import pickle

    rng = crng.Ranlux(luxury=4, seed=54987)
    rng.compute(1000)

    file = open('rng.save', 'w')
    pickler = pickle.Pickler(file)
    pickler.dump(rng)
    rng.compute(10)
    print rng(2)                    # (0.77230155468, 0.246806561947)

The second script retrieves the stored crng instance from the file using pickle. Note that this script does not need to import neither the crng_pickle nor the crng modules.

    import pickle

    file=open('rng.save')
    unpickler = pickle.Unpickler(file)
    rng = unpickler.load()
    rng.compute(10)
    print rng(2)                    # (0.77230155468, 0.246806561947)


C interface

The crng module provides an interface to allow other Python extension modules written in C to access the C implementation of the RNGs directly. Notice that the interface is not designed to allow any C code to access the crng functions, but only to provide access via Python objects to the crng implementation for Python extension modules written in C.

The interface consists of a CObject containing a pointer to the C function for producing the next value for the RNG. The CObject is a private member of the crng instance. It is not visible in the __members__ attribute of the instance. In order to use it, the application Python code must provide the crng instance as an argument to the client Python function (a Python extension coded in C). The C implementation of the client function fetches the CObject from the instance, and gets the C pointer to the RNG function, which can then be used in the C code.

An example shows how the Python application code creates a crng instance, which is passed to the client Python extension, in this case a Monte Carlo integrator:

  import crng
  import MonteCarloIntegrator
  rng = crng.MT19973(seed=123)
  mci = MonteCarloIntegrator(func=function_to_integrate, rng=rng)
  integral = mci.integrate(region=(0,1))

And here are the relevant parts of the C code for the object initialization function in the Python extension module MonteCarloIntegrator:

  static PyObject *
  MonteCarloIntegrator_new (PyObject *self, PyObject *args)
  {
    PyObject *func, *rng, *cobj;
    double (*rng_next_value) (PyObject *rng);
    double rng_value;

    if (! PyArg_ParseTuple(args, "OO", &(func), &(rng)))
      goto error;

    /* ... code removed for clarity */

    cobj = PyObject_GetAttrString(rng, "_crng_basic_next_value");
    if (cobj && PyCObject_Check(cobj)) {
      rng_next_value = (double (*)(PyObject *)) PyCObject_AsVoidPtr(cobj);
    } else {
      PyErr_SetString(PyExc_TypeError,
                      "rng must be a basic RNG from module crng");
      goto error;
    }

    /* ... code removed for clarity */

    rng_value = rng_next_value(rng);

    /* ... code removed for clarity */

  }

The C function prototypes and return values are the same for the basic RNGs, but may be different for the non-basic deviate RNGs, depending on the type of the deviate. The CObject member name and the declaration for the C function pointer it contains are given for each RNG below. Note that the same CObject member name may return pointers to a C function with different prototypes depending on the RNG.


Scripts

A few Python scripts are included in the distribution.

validate.py

The script validate.py tests whether the crng implementation of a few specific RNGs produces values that agree with values given in the original reference for the RNG. Thus, the tests validate both the portability and the current implementation of the RNGs.

Currently, the following RNGs are validated by the script: ParkMiller, WichmannHill, Ranlux and MT19937.

profile.py

The script profile.py contains code to perform profiling (CPU timing) measurements for the RNGs. The results for two different systems are given in the profiling table.

display.py

The script display.py uses the Tkinter module to display some non-standard distributions. Using the radiobuttons in the menu, it is possible to change the distribution for one or both axes of the coordinate system. For each change, a new set of random values are generated and the corresponding points are displayed.

crng_pickle.py

The script crng_pickle.py performs the setup required to use the standard Python pickle/cPickle modules with crng instances. It uses the standard Python copy_reg module to do this.

This script should be installed in the same location as the crng module compiled shared-object file crng.so. This script needs only to be imported; no other actions are required to perform the setup. See the section Saving and restoring an RNG object for information on how to use it.


Profiling (CPU timings)

The CPU execution times were measured by running the compute function of each RNG for at least 5 million iterations. This measures the speed of the C implementation, with minimal overhead for the Python system. When using the next function in Python code, the execution speed will be significantly slower due to the overhead of the Python system.

The script profile.py contains code to perform profiling measurements.

RNG Pentium III 500 MHz, Linux Red Hat 6.1, gcc (2.91.66) -g -O2 SGI O2 (MIPS R5200), IRIX 6.5, cc MIPSpro 7.30 -O -n32
106 iterations per second microseconds per iteration 106 iterations per second microseconds per iteration
ParkMiller() 5.5 0.18 4.8 0.21
WichmannHill() 1.7 0.59 1.0 0.96
LEcuyer() 4.0 0.25 3.3 0.30
Ranlux(luxury=0) 7.4 0.14 6.0 0.17
Ranlux(luxury=1) 5.1 0.20 3.7 0.27
Ranlux(luxury=2) 3.1 0.32 2.2 0.45
Ranlux(luxury=3) 1.5 0.65 1.1 0.94
Ranlux(luxury=4) 0.93 1.1 0.63 1.6
Taus88() 9.8 0.10 6.4 0.16
MRG32k3a() 1.7 0.59 2.0 0.51
MT19937() 6.0 0.17 5.0 0.20
rand
(C stdlib; not in crng)
6.1 0.16 4.4 0.23
whrandom.py
(Python standard library)
0.065 15.43 - -


Changes in version 1.2

The changes in version 1.2 are only fixes for bugs. No new features have been implemented.


Changes in version 1.1

The changes in version 1.1 are fairly minor. Nothing has changed in the basic functionality and interface, nor are there any new RNG implementations.


To do


References

M. Galassi, J. Davies, J. Theiler, B. Gough, R. Priedhorsky, G. Jungman, M. Booth, "GNU Scientific Library 0.5" (1999). http://sources.redhat.com/gsl/.

J.E. Gentle (1998), "Random Number Generation and Monte Carlo Methods", Springer, ISBN 0-387-98522-0.

D.E. Knuth (1998), "The Art of Computer Programming", Vol. 2, "Seminumerical Algorithms", 3rd edition, Addison-Wesley, ISBN 0-201-89684-2. Chapter 3: Random numbers.

W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery (1992) "Numerical Recipes in C", 2nd edition, Cambridge University Press, ISBN 0-521-43108-5. Chapter 7: Random numbers.


Acknowledgments

Many thanks for miscellaneous contributions to this module since its first release: Alex Martelli, Paul Moore, Tim Peters, Anton Vredegoor.


Basic RNGs: uniform (0,1) distribution

All basic RNGs implemented in the crng module produce values uniformly distributed in the open interval (0,1), i.e. exclusive of the end-point values 0 and 1.

ParkMiller
Initialization r = ParkMiller(seed=314159265)

seed: initial seed, positive integer

Members Same as in the initialization call; read-only.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period 2.1*109
Reference Stephen K. Park and Keith W. Miller, Communication of the ACM (1988) 31, 1192-1201. "Random Number Generators: Good Ones Are Hard to Find."
Comment This is the 'minimal standard' RNG proposed by the authors in the reference. The implementation is a translation of the Pascal code given as "Integer Version 2" in the reference.

WichmannHill
Initialization

r = WichmannHill(seed1=314, seed2=159, seed3=365)

seed1: initial first seed, integer such that 0<seed<=30268

seed2: initial second seed, integer such that 0<seed<=30306

seed3: initial third seeds integer such that 0<seed<=30322

Members Same as in the initialization call; read-only.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period 6.95*1012
References B.A. Wichmann and I.D. Hill, Applied Statistics (1982) 31, 188-190. "Algorithm AS 183. An Efficient and Portable Pseudo-random Number Generator."

B.A. Wichmann and I.D. Hill, Applied Statistics (1984) 33, 123. "Correction: Algorithm AS 183"

See also A. Ian McLeod, Applied Statistics (1985) 34, 198-200. "A Remark on Algorithm AS 183."

The standard Python module whrandom.py.

Comment The McLeod paper points out that the algorithm may produce some zero values due to rounding errors under certain circumstances when 32 bit precision is used. The proposed amendment to fix this problem has not been implemented, since the computations in this implementation are done in double precision.

LEcuyer
Initialization r = LEcuyer(seed1=314159265, seed2=314159263)

seed1, seed2: initial seeds, positive integers

Members Same as in the initialization call; read-only.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period 2.3*1018
Reference Pierre L'Ecuyer, Communications of the ACM (1988) 31, 742-749+774. "Efficient and Portable Combined Random Number Generators."
Web site http://www.iro.umontreal.ca/~lecuyer/
Comment This is a combination of two Multiplicative Congruential Generators (MLCGs). The implementation is a translation of the code from fig. 3 of the reference.

Ranlux
Initialization

r = Ranlux(luxury=3, seed=314159265, state=None)

luxury: the quality level, integer such that:
0: Equivalent to original RCARRY by Marsaglia and Zaman. Has a very long period, but fails many tests.
1: Considerable improvement in quality over level 0. Passes the gap test, but still fails the spectral test.
2: Passes all known tests, but theoretically still defective.
3: Any theoretically possible correlations have a very small chance of being observed.
4: Highest possible luxury, all 24 bits chaotic.

seed: initial seed, positive integer

state: set the initial state from a state tuple (available as a member in another Ranlux instance). This argument, if given, overrides the luxury and seed arguments. Warning: only values obtained from the state member in another Ranlux instance should be used; attempts to use a state tuple built from scratch may result in a dangerously invalid instance.

Members

luxury: the quality level; read-only

state: the tuple reflecting the current state; read-only

Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period approximately 5.2*10171
References Martin Lüscher, Computer Physics Communications (1994) 79, 100-110. "A portable high-quality random number generator for lattice field theory simulations."

F. James, Computer Physics Communications (1994) 79, 111-114. "RANLUX: A Fortran implementation of the high-quality pseudorandom number generator of Lüscher."

F. James, Computer Physics Communications (1996) 97, 357. "Erratum. RANLUX: A Fortran implementation of the high-quality pseudorandom number generator of Lüscher."

See also Kenneth G. Hamilton and F. James, Computer Physics Comunications (1997) 101, 241-248. "Acceleration of RANLUX."
Comment The implementation is based on the Fortran 77 version written by Fred James, as posted to a Usenet News group on 16 March 1995 by Byron Bodo.

Taus88
Initialization r = Taus88(seed1=314159265, seed2=314159263, seed3=314159261)

seed1: first initial seed, integer >=2
seed2: second initial seed, integer >=8
seed3: third initial seed, integer >=16

Members Same as in the initialization call; read-only.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period Approximately 288 = 3.1*1026
Reference Pierre L'Ecuyer, Mathematics of Computation (1996) 65, 203-213. "Maximally equidistributed combined Tausworthe generators."
Web site http://www.iro.umontreal.ca/~lecuyer/
Comment The implementation is from figure 1 of the reference.

MRG32k3a
Initialization r = MRG32k3a(s10=314159265, s11=314159263, s12=314159261, s20=314159259, s21=314159257, s22=314159255)

s10, s11, s12: initial seed vector s1 elements, positive integers
s20, s21, s22: initial seed vector s2 elements, positive integers

Members Same as in the initialization call; read-only.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period Approximately 2191 = 3.1*1057
Reference Pierre L'Ecuyer, Operations Research (1999) 47, 159-164. "Good parameters and implementations for combined multiple recursive random number generators."
Web site http://www.iro.umontreal.ca/~lecuyer/
Comment The implementation is from figure 1 of the reference.

MT19937
Initialization

r = MT19937(seed=314159265, state=None)

seed: initial seed, positive integer

state: set the initial state from a state tuple (available as a member in another MT19937 instance). This argument, if given, overrides the seed argument. Warning: only values obtained from the state member in another MT19937 instance should be used; attempts to use a state tuple built from scratch may result in a dangerously invalid instance.

Members

state: the tuple reflecting the current state; read-only

Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_basic_next_value or _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Period 219937-1 = approximately 4.3*106001
Reference Makoto Matsumoto and Takuji Nishimura, ACM Transactions on Modeling and Computer Simulation (1998) 8, 3-30. "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator."
Web site http://www.math.keio.ac.jp/matumoto/emt.html
Comment The implementation uses the C code available at the web site (above), which is copyrighted by Matsumoto and Nishimura, and has been released by them under the GNU Lesser General Public License (LGPL).
Copyright Copyright © 1997, 1999 Makoto Matsumoto and Takuji Nishimura.

When you use this, send an email to: matumoto@math.keio.ac.jp with an appropriate reference to your work.

The MT19937 code by Matsumoto and Nishimura is distributed under the GNU Lesser General Public License (formerly known as Library General Public License) as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with this library (file lgpl.txt); if not, write to the Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.


Deviates of continuous distributions

The deviates of the continuous distributions other than the uniform (0,1) distribution are computed by an appropriate mathematical transformation of values taken from the uniform (0,1) distribution. The extension types for the deviates have been designed so that any basic RNG from the crng module can be used. The choice is made at initialization, but may be changed after creation. Note that only basic RNGs are valid in this context; e.g. a GammaDeviate object cannot be used as a RNG for another GammaDeviate.

UniformDeviate
Initialization r = UniformDeviate(rng=ParkMiller(), a=0.0, b=1.0)

rng: a basic RNG from the crng module
a, b: the limits of the open interval, such that a<b

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)

ExponentialDeviate
Initialization r = ExponentialDeviate(rng=ParkMiller(), mean=1.0)

rng: a basic RNG from the crng module
mean: the deviate mean value, must be positive

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Description For a radioactive substance with a long half-life, the time between successive decay events is drawn from an exponential distribution. Similarly, the waiting times between the arrival of customers to an office can be described by an exponential distribution, provided the average influx is constant. Or, for a set of radioactive nuclei at time 0, the decay events occur at the times {X}, where {X} are exponential deviates. The half-life of the exponential is related to the mean by t0.5=mean*ln(2).

NormalDeviate
Initialization r = NormalDeviate(rng=ParkMiller(), mean=0.0, stdev=1.0, stored=None)

rng: a basic RNG from the crng module
mean: the deviate mean value
stdev: the deviate standard deviation, must be positive
stored: next raw random value (if any); saved from the previous calculation

Members

Same as in the initialization call; read/write, except:

stored: the next raw (mean and stdev have not been applied) random value (if any), saved from the previous calculation (every other call); read-only

Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Description Also called the Gaussian deviate; the density function is a Gaussian.
References G.E.P. Box and M.E. Muller, Annals Math. Stat. (1958) 29, 610-611. "A note on the generation of random normal deviates."

GammaDeviate
Initialization r = GammaDeviate(rng=ParkMiller(), order=2.0, scale=1.0, direct=(order<=12))

rng: a basic RNG from the crng module
order: the order of the gamma deviate, must be positive
scale: the scale factor; must be positive
direct: use the direct method (log of the sum of exponential random values), else the rejection method; integer interpreted as boolean

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Description If an event is the product of two or more sub-events, each having an exponential distribution, then that event has a gamma distribution. The number of sub-events is the order of the gamma distribution.
Reference J.H. Ahrens, Ann. Inst. Stat. Math. (1962) 13, 231-237.
Comment The optimal choice between the direct method and the rejection method for calculating the deviate depends to some extent on the properties of the compiler and the hardware. The default choice made in the current implementation is based on profiling (CPU timings) on a Pentium III 500 MHz system running Linux Red Hat 6.1, and on a SGI O2 (MIPS R5200) system running IRIX 6.5. Note that using a fractional order value significantly decreases the computational efficiency of the direct method.

BetaDeviate
Initialization r = BetaDeviate(rng=ParkMiller(), a=1.0, b=1.0)

rng: a basic RNG from the crng module
a, b: the parameters, must be positive

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects
C interface Member name: _crng_next_value
C function pointer: double (*next_value) (PyObject *rng)
Comment The calculation uses two gamma deviates which in turn call the specified basic RNG. The choice of using the direct method or the rejection method for calculating the gamma deviates is made automatically based on the values of the parameters a and b.


Deviates of integer-valued distributions

The deviates of the integer-valued distributions are computed by an appropriate mathematical transformation of values taken from the uniform (0,1) distribution. The extension types for the deviates have been designed so that any basic RNG from the crng module can be used. The choice is made at initialization, but may be changed after creation. Note that only basic RNGs are valid in this context; e.g. a PoissonDeviate object cannot be used as a RNG for another PoissonDeviate.

PoissonDeviate
Initialization r = PoissonDeviate(rng=ParkMiller(), mean=1.0, direct=(mean<12.0))

rng: a basic RNG from the crng module
mean: the deviate mean value, must be >=0.0
direct: use the direct method, else the rejection method; integer interpreted as boolean

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects, except that the functions that normally produce floats return integers.
C interface Member name: _crng_next_value
C function pointer: long (*next_value) (PyObject *rng)
Description For a radioactive substance with a long half-life, the number of decay events in a given time span is drawn from a Poisson distribution. Similarly, the number of customers arriving to an office in a given time span can be described by a Poisson distribution, provided the average influx is constant.
Comment The optimal choice between the direct method and the rejection method for calculating the deviate depends to some extent on the properties of the compiler and the hardware. The default choice made in the current implementation is based on profiling (CPU timings) on a Pentium III 500 MHz system running Linux Red Hat 6.1, and on a SGI O2 (MIPS R5200) system running IRIX 6.5.

BinomialDeviate
Initialization r = BinomialDeviate(rng=ParkMiller(), p=0.5, n=2)

rng: a basic RNG from the crng module
p: the probability of the event, in the interval [0,1]
n: the number of trials, positive integer

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects, except that the functions that normally produce floats return integers.
C interface Member name: _crng_next_value
C function pointer: long (*next_value) (PyObject *rng)
Description If an event occurs with a given probability, and a given number of trials are made, then the binomial distribution describes the number of successful trials.
Comment The algorithm implemented contains a tunable parameter, which affects the computational efficiency. Currently, this parameter is hardwired, and its value was chosen by empirical tests.

GeometricDeviate
Initialization r = GeometricDeviate(rng=ParkMiller(), p=0.5)

rng: a basic RNG from the crng module
p: the probability of the event, in the interval (0,1]

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects, except that the functions that normally produce floats return integers.
C interface Member name: _crng_next_value
C function pointer: long (*next_value) (PyObject *rng)
Description The number of independent trials to yield a success, for a given probability of the success of each trial.

BernoulliDeviate
Initialization r = BernoulliDeviate(rng=ParkMiller(), p=0.5)

rng: a basic RNG from the crng module
p: the probability, in the interval [0,1]

Members Same as in the initialization call; read/write.
Methods (), next(), compute(), density(), str(), repr()
See Methods available for all RNG objects, except that the functions that normally produce floats return integers.
C interface Member name: _crng_next_value
C function pointer: long (*next_value) (PyObject *rng)
Description The repeated flipping of a coin produces a Bernoulli series for a probability 0.5 if heads and tails is represented by 1 and 0, respectively.


Sampling and shuffling

The following are utility functions for sampling and shuffling of a sequence object (string, tuple or list). For choose, sample and shuffle, the result is a sequence of the same type as the input sequence.

choose
Usage newseq = choose(seq, n, rng)

newseq: a new sequence of the same type as seq
seq: the sequence (string, tuple or list) from which to choose
n: the number of objects to choose, such that 0<=n<=len(seq)
rng: a basic RNG from the crng module

Description The objects in the output sequence are chosen without replacement from the input sequence. The relative order of the objects in the output sequence is preserved.

sample
Usage newseq = sample(seq, n, rng)

newseq: a new sequence of the same type as seq
seq: the sequence (string, tuple or list) from which to sample with replacement
n: the number of objects to sample, such that 0<=n<=len(seq)
rng: a basic RNG from the crng module

Description The objects in the output sequence are chosen with replacement from the input sequence. The output sequence may contain multiple references to the same object. The relative order of the objects in the output sequence is random.

shuffle
Usage newseq = shuffle(seq, rng)

newseq: a new sequence of the same type as seq
seq: the sequence (string, tuple or list) to shuffle
rng: a basic RNG from the crng module

Description The objects in the output sequence are the same as in the input sequence, but in random order compared to the input sequence.

stir
Usage stir(list, rng)

list: the list to be shuffled
rng: a basic RNG from the crng module

Description The list is shuffled in-place.

pick
Usage obj = pick(seq, rng)

obj: the object picked from seq
seq: the sequence (string, tuple or list) from which to pick one object
rng: a basic RNG from the crng module

Description One object is picked at random from the given sequence.


$Date: 2002/10/23 12:11:20 $