Alfred A. Brooks   100 Wiltshire Drive   Oak Ridge, TN 37830-4505
   Phone (& fax/call): 423 482-1559   E-mail: brooks50@comcast.net
File: test.doc					Date: 03/28/97 6:49:09 PM


Random Number Generators from Ranked Order Distribution Functions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The probabilistic approach is often limited or abandoned due to the lack 
of distribution functions. The following is an approach [1] which will 
give sequences which approximate a required cumulative distribution 
function and which improve in accuracy in the limit as the number of 
points increases. The desired cumulative probability function (CPF) for 
a set of n+1 data points is approximated by a cumulative probability 
sequence (CPS) comprising n intervals of uniform probability. The CPS 
may be used as a random number generator or in other ways.

Consider the set of n+1 data points drawn from the same varying 
population: x0, x1, ..., xn. Place them in rank order to form the 
sequence: x(0), x(1), ...x(k),... , x(n). Let the probability finding 
exactly one measurement in the range x(k) to (x(k+1) be 1/n. Then, the 
cumulative probability sequence of an observation being less than x(k) 
is: (k/n), i.e., the CPS is: 0, 1/n, 2/n, ..., k/n, ..., n/n=1. 

The CPF can be approximated by linear interpolation in the CPS and the 
equation is:

           j-1      (x-x(j-1))
F(x'<=x) = ---  + ---------------,  x(j-1)<x<=x(j),  j=1, 2, 3,...,n
            n     n*(x(j)-x(j-1))

This cumulative probability sequence has the following properties:

1) In the limit it approaches the corresponding cumulative distribution 
function.

2) The CPS and CPF are bounded by x(0) and x(n) which tend to the 
limiting bounds.

3) The cumulative probability sequence can be regarded as an equal 
probability histogram for which there is a simple and fast random number 
generator algorithm [2] based on the inverse distribution technique.

Inverse Distribution Theorem: If F(x) is a cumulative distribution 
function such that the probability of a measurement m <= F(m), and f(x) 
is its corresponding probability density function such that a 
measurement, m, lies between m-dx and m+dx, then if F(x) is distributed 
uniformly between 0 and 1, x is distributed randomly according to f(x).

For Monte Carlo simulation, there is little need for the probability 
density forms and they are included only for completeness.

The cumulative distribution function can be converted into its 
probability density function (PDF) by differentiation of the CPF.


The derivative in the i-th interval and the PDF are:


                1
f(x+-dx) = -------------,    x(j-I) < x < x(j),   j = 1, 2, 3,   ., n
           (x(j)-x(j-1))

This form may be converted to the more familiar uniform interval of 
measurement form by combining the appropriate panels into uniform 
intervals partitioning as necessary. This form has the following 
characteristics.

1) The probability density sequence can be regarded as a probability 
histogram and can be used as a random number generator2 but it is not as 
fast as the CPS form.

2) As the number of points becomes large it approaches the corresponding 
probability density function as a limit.

3) The sequence and distribution function (in the limit) are bounded by 
x(0) and x(n).

4) For n=0, the distribution function degenerates to a single point 
value, x(0).

5) For n=1, the distribution function is a uniform distribution of range 
x(0) to x(n).

6)For n>1, the distribution function is an n-panel histogram that may 
approximate but does not necessarily correspond precisely to any simple 
distribution function.

This method has the advantage that no assumption is made concerning the 
shape of the distribution function other than it is linear over short 
distances. The method also has an advantage for use in Monte Carlo 
methods in that the random number generator is faster than most 
selection techniques and can be programmed to obtain n as the first 
member of the CPS array. The table can then be expanded at will simply 
by editing its source file as new data becomes available.

[1]  I am indebted to Kenneth Redus of Macteck for bringing this method to 
my attention.
[2]  J. C. McGrath and D.C. Irving; Techniques for Efficient Monte Carlo 
Simulation, Vol. II - Random Number Generation; ORNL-ORSIC-38, p. 88; 
Apr 1975; Oak Ridge National Laboratory, Oak Ridge TN, 37830



The following is a random number generator based on the CPF described 
above:

/* ------------------------------------------------- *
 * RANKRN.C - Ranked Data Inverse Distribution       *
 * Given a CPF, F(x), and Uniform Random Number it   *
 * solves for x, a random number from the PDF, f(x). *
 * F(x) is the n+1 data values of a CPF having       *
 * constant propability intervals = 1/n              *
 * Data Values in xa: (float)n, x[0], x[1],..., x[n] *
 * Validated by variance & moments 1-4               *
 * ------------------------------------------------- */
#include <stdlib.h>
float urand(void); /* Uniform Ramdom Number Generator (0 to 1) */
float rankrn (float *xa); /* Returns a random number from f(x) */

float rankrn(float *xa)
{int i; float r;
 /* In F(x) calculate the interval index & partial interval  reset 
index*/
 r=xa[0]*urand(); i=(int)r; r-=(float)i; i++;
 /* Return the corresponding interpolated data value from f(x) */
 return(xa[i]+r*(xa[i+1]-xa[i]));}
/* ------------------------------------------------- */