💾 Archived View for gemini.spam.works › mirrors › textfiles › programming › random.txt captured on 2022-04-29 at 00:33:12.
⬅️ Previous capture (2020-10-31)
-=-=-=-=-=-=-
Computer Generated Random Numbers
David W. Deley
deleyd@netcom.com
Copyright 1991
SYNOPSIS
This paper presents a number of techniques for analyzing a computer
generated sequence of random numbers. Some background theory in
probability theory and inferential statistics is presented, and a
number of empirical tests are presented along with example programs
which apply the tests to several popular random number generators
including the VAX Run-Time Library routine MTH$RANDOM used by the
FORTRAN intrinsic function RAN and by the VAX BASIC function RND, the
standard ANSI C rand() function used by VAX C, and the Data Encryption
Standard (DES) algorithm. For comparison the sample programs also
apply the tests to the obsolete RANDU random number generator.
CONTENTS:
1.0 INTRODUCTION
1.1 MTH$RANDOM USAGE
1.2 RANDU USAGE
2.0 PROBABILITY THEORY
2.1 ROLLING DIE
2.2 THE LANGUAGE OF PROBABILITY
3.0 INFERENTIAL STATISTICS
3.1 THE PROBABILITY OF ERROR
3.2 GENERAL TEST PROCEDURE
3.3 THE CHI-SQUARE TEST
3.4 THE KOLMOGOROV-SMIRNOV TEST
4.0 ANALYSIS OF A SEQUENCE OF RANDOM NUMBERS
4.1 ONE-DIMENSIONAL TESTS
4.2 TWO-DIMENSIONAL TESTS
4.3 THREE-DIMENSIONAL TESTS
4.4 HIGHER DIMENSIONAL TESTS
5.0 LINEAR CONGRUENTIAL RANDOM NUMBER GENERATORS
5.1 THE MTH$RANDOM ALGORITHM
5.2 THE RANDU ALGORITHM
5.3 STARTING THE MTH$RANDOM GENERATOR
5.4 STARTING THE RANDU GENERATOR
5.5 STARTING THE ANSI C GENERATOR
6.0 RATINGS FOR VARIOUS GENERATORS
6.1 REVIEW OF TESTS PERFORMED
6.2 TEST RESULTS
6.2.1 MTH$RANDOM
6.2.2 RANDU
6.2.3 C and ANSI C
6.2.4 Microsoft C
6.2.5 Turbo Pascal
7.0 THE SPECTRAL TEST
8.0 SHUFFLING
9.0 DES AS A RANDOM NUMBER GENERATOR
10.0 HOW GOOD A RANDOM GENERATOR DO YOU NEED?
11.0 WHAT IS A RANDOM NUMBER GENERATOR?
12.0 CONCLUSION
REFERENCES
---------------------------------------------------------------------------
1.0 INTRODUCTION
The VMS Run-Time Library provides a random number generator routine
called MTH$RANDOM. On the first call you supply an initial SEED value
to the routine, and it returns a random deviate uniformly distributed
between [0,1). Call the routine again and it returns another random
deviate uniformly distributed between [0,1). Calling the routine over
and over again produces a sequence of random deviates all supposedly
uniformly distributed between [0,1). Our question now is, is this
sequence of numbers really random?
Although there is nothing "random" about a completely deterministic
computer generated sequence of numbers, we can analyze the sequence of
numbers to see if there is any reason to doubt that the sequence could
have come from a truly random stochastic process.
1.1 MTH$RANDOM USAGE
The MTH$RANDOM routine is called by the FORTRAN intrinsic function RAN.
>From FORTRAN it is called as follows:
INTEGER SEED
R = RAN(SEED)
1.2 RANDU USAGE
The VMS FORTRAN Run-Time Library also contains a second random number
generator which implements the infamous RANDU random number generator
first introduced by IBM IN 1963. This turned out to be a poor random
number generator, but nonetheless it has been widely spread. If you
should ever come across it on some other machine, don't use it. VMS
FORTRAN abandoned the RANDU generator in 1978 in favor of the MTH$RANDOM
routine. The RANDU routine is considered obsolete and is now completely
undocumented, but it still exists in the VMS FORTRAN Run-Time Library
and will always exist to support any old programs which still use it.
>From FORTRAN the obsolete RANDU random number generator can be called in
any of three ways:
INTEGER*4 SEED
INTEGER*2 W(2)
EQUIVALENCE( SEED, W(1) )
R = FOR$IRAN( W(1), W(2) ) !Preferred way
CALL FOR$RANDU( W(1), W(2), R ) !Alternate way
CALL FOR$RANDU_W( W(1), W(2), R ) !Alternate way
R is the return value between [0,1). W(1) and W(2) together is the
seed value for the generator. This goes back to the PDP-11 days of 16
bit integers. SEED is really a 32 bit integer, but it was represented
as two 16 bit integers.
We will first cover some background theory in probability theory, and
then we will test the FORTRAN RAN and RANDU generators. Hopefully we
will find out why the infamous RANDU generator is considered bad.
---------------------------------------------------------------------------
2.0 ROLLING DIE
Consider the case of rolling a single die. A theoretician picks up the
die, examines it, and makes the following statement: "The die has six
sides, each side is equally likely to turn up, therefore the
probability of any one particular side turning up is 1 out of 6 or 1/6.
Furthermore, each throw of the die is completely independent of all
previous throws." With that the theoretician puts the die down and
leaves the room without ever having thrown it once.
An experimentalist, never content to accept what a cocky theoretician
states as fact, picks up the die and starts to roll it. He rolls the
die a number of times, and carefully records the results of each throw.
Let us first dispel the myth that there is such a thing as an ideal
sequence of random numbers. If the experimentalist threw the die 6
times and got 6 "one"s in a row, he would be very upset at the die. If
the same experimentalist threw the same die a billion times in a row and
never got 6 "one"s in a row, he would again be very upset at the die.
The fact is a truly random sequence will exhibit local nonrandomness.
Therefore, a subsequence of a sequence of random numbers is not
necessarily random. We are forced to conclude that no one sequence of
"random" numbers can be adequate for every application.
2.1 THE LANGUAGE OF PROBABILITY
Before we can analyze the resulting data the theoretician collected, we
must first review the language of probability theory and the underlying
mathematics. (Unfortunately, the limitations of my character cell
terminal prevent me from directly typing many of the mathematical
symbols needed in this discussion. So I'll make due as best I can and
hope the reader can follow me.)
A single throw of the die is called a "chance experiment" and is
designated by the capital letter E. An outcome of a chance experiment
E is designated by the Greek letter ZETA. If there is more than one
possible outcome of the chance experiment, (as there always will be,
else the analysis becomes trivial), the different possible outcomes
of the chance experiment are designated by ZETA_1, ZETA_2,...
The set of all possible outcomes of a chance experiment is called the
"Universal Set" or "Sample Space", and is denoted by the capital
letter S.
For the die throwing experiment E, we have:
Chance experiment:
E: one throw of the die
Possible outcomes:
ZETA_1 = . (one dot)
.
ZETA_2 = (two dots)
.
.
ZETA_3 = . (three dots)
.
. .
ZETA_4 = (four dots)
. .
. .
ZETA_5 = . (five dots)
. .
. . .
ZETA_6 = (six dots)
. . .
Sample Space (or Universal Set):
S = {ZETA_1, ZETA_2, ZETA_3, ZETA_4, ZETA_5, ZETA_6}
We now define the concept of a "random variable". A random variable is
a function which maps the elements of the sample space S to points on
the real number line RR. This is how we convert the outcomes of a
chance experiment E to numerical values. The random variable function
is denoted by a capital X. An actual value a random variable X takes
on is denoted by a lowercase x.
A natural random variable function X to define for the die rolling
experiment is:
X( ZETA_1 ) = 1.0
X( ZETA_2 ) = 2.0
X( ZETA_3 ) = 3.0
X( ZETA_4 ) = 4.0
X( ZETA_5 ) = 5.0
X( ZETA_6 ) = 6.0
Note that the faces of a standard die are marked with this random
variable function X in mind. The number of dots showing on the top
face of the die corresponds to the value the function X takes on for
that particular outcome of the experiment.
We now consider the generation of a random sequence of numbers
by repeating a chance experiment E a number of times. This can be
thought of as a single n-fold compound experiment, which can be
designated by:
E x E x . . . x E = En
(N times)
The sample space or universal set of all possible outcomes for our
compound experiment En is the Cartesian product of the sample spaces
for each individual experiment E:
Sn = S x S x . . . x S
(N times)
Any particular result of a compound experiment En will be an ordered
set of elements from set S. For example, if E is the die throwing
experiment, then the elements of set Sn for rolling the die N times are:
Sn = { ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_1 }
{ ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_2 }
{ ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_1, ZETA_3 }
. . . . . .
. . . . . .
{ ZETA_6, ZETA_6, ZETA_6, ZETA_6, ZETA_6, ZETA_6 }
The total number of elements in Sn is equal to the number of elements
in S raised to the Nth power. For the die throwing experiment, there
are 6 possible outcomes for the first throw, 6 possible outcomes for
the second throw, etc. so the total number of possible outcomes for
the compound experiment En is:
6 * 6 * . . . * 6 = 6**N = |Sn|
(N times)
The probability of any particular element in Sn is equal to the product
of the probabilities of the corresponding elements from S which
together comprise the element in Sn. If each element of S is equally
likely, then each element of Sn is equally likely. Which translated
means if each possible outcome of experiment E is equally likely, then
each possible outcome of experiment En is equally likely. For our die
throwing experiment, each of the 6 possible outcomes for throwing the
die once is equally likely, so each of the 6**N possible outcomes for
throwing the die N times is equally likely.
We can now take our random variable function X and apply it to each of
the N outcomes from our compound experiment En. This is how we convert
the outcome of experiment En into an ordered set of N random variables:
( X1, X2, ..., Xn )
In this way each specific element of our sample space Sn can be
transformed into a set of n ordered real values. We will use a
lowercase letter x to denote specific values that a random variable
function X takes on. So the results of a specific compound experiment
En when transformed into specific real values would look something
like this:
( x1, x2, ..., xn )
For example, if we rolled a die N times we might get the result:
( 3, 5, 2, 6, ..., 4 )
(N times)
Then again, if we rolled a die N times we might get the result
consisting of all ones:
( 1, 1, 1, 1, ..., 1 )
(N times)
The first result looks somewhat like what we would expect a sequence of
random numbers to be. The second result consisting of all ones however
tends to make us very unhappy. But, and here is the crux, EACH
PARTICULAR OUTCOME IS EQUALLY LIKELY TO OCCUR! Each outcome
corresponds to a particular element in Sn, and each element in Sn has
an equal probability of occurring. We are forced to conclude that ANY
sequence of random numbers is as equally likely to occur as any other.
But for some reason, certain sequences make us very unhappy, while
other sequences do not.
We are forced to conclude, therefore, that our idea of a good sequence
of random numbers is more than just anything you might get by
repeatedly performing a random experiment. We wish to exclude certain
elements from the sample space Sn as being poor choices for a sequence
of random numbers. We accept only a subspace of Sn as being good
candidates for a sequence of random numbers. But what is the
difference between a good sequence of random numbers and a bad sequence
of random numbers?
For our die rolling experiment, if we believe the die to be fair, and
roll the die N times, then we know all possible sequences of N outcomes
are equally likely to occur. However, only a subset of those possible
outcomes will make us believe the die to be fair. For example, the
sequence consisting of all ones {1,1,...,1} is as equally likely to
occur as any other sequence, but this sequence tends not to reinforce
our assumption that the die is fair. A good sequence of random numbers
is one which makes us believe the process which created them is random.
---------------------------------------------------------------------------
3.0 INFERENTIAL STATISTICS
Up until now we have been dealing with probability theory. We started
with a chance experiment in which the probability of each result is
given, and then worked our way forward to find out what the probability
of a particular compound experiment would be. We now turn to the field
of inferential statistics, where we attempt to do the reverse. Given a
sequence of numbers created by N trials of a chance experiment, we
attempt to work our way backwards and estimate what the defining
properties of the chance experiment are.
3.1 THE PROBABILITY OF ERROR
For our die throwing experiment, we would like to estimate the
probability of each outcome based on experimental data. If we throw
the die N times, then the best estimate we could make for the
probability of a particular outcome ZETA would be to count the number
of times ZETA occurred in N trials and form the ratio:
number of times ZETA occurred
------------------------------
total number of trials
If we define a frequency function f(ZETA) to be the number of times
ZETA occurred in N trials, we then would have:
f(ZETA)
Pr(ZETA) = -------
N
This is our estimate of the probability of ZETA based on N trials of the
chance experiment. Note that some people use this equation as the
definition for the probability of an event ZETA when N goes to infinity:
_ f(ZETA)
Pr(ZETA) = lim -------
N -> infinity N
For our die throwing experiment, assume we want to estimate the
probability of ZETA_1, (the side with one dot shows up). Say we throw
the die 6 times, and out of those 6 times ZETA_1 occurred twice. We
then calculate the probability of ZETA_1 as 2/6.
We now have the problem of deciding which value to accept as correct,
and how much faith we can place on our choice. Let us define two
possible hypotheses, which we will denote as H0 and H1.
H0 : The probability of ZETA_1 is 1/6
H1 : The probability of ZETA_1 is 2/6
In choosing which hypothesis to believe, there are two possible
errors we could make, known as a type I error and a type II error:
Type I error : Choose H1 when H0 is actually true.
Type II error : Choose H0 when H1 is actually true.
We can determine the probability of making each type of error by
examining the theoretical probability distributions for f(ZETA)
for each hypothesis being tested.
We rolled the die six times, so the possible number of times ZETA_1 can
occur is between 0 and 6 times. Let f(ZETA_1) be the number of times
ZETA_1 occurred in six rolls of the die. The following binomial
formula will calculate the probability of having ZETA_1 occur k times
in n trials:
/ n \
Pr( f(ZETA_1) = k ) = | | * p**k * q**(n - k)
\ k /
Where:
Pr( f(ZETA_1) = k ) = probability of exactly k occurrences
of ZETA_1 in n trials
f(ZETA_1) = number of occurrences of ZETA_1 in n trials
n = number of trials
p = Pr(ZETA_1) (probability of ZETA_1)
q = 1 - p
/ n \ n!
| | = ----------- (n things taken k at a time)
\ k / (n - k)! n!
With n = 6 trials, and assuming hypothesis H0 to be true (the
probability of throwing ZETA_1 is 1/6), we obtain the following table:
f(ZETA_1) Pr ( f(ZETA_1) | H0 ) (sideways bar graph)
----------- ----------------------
0 0.334897 |================
1 0.401877 |====================
2 0.200938 |==========
3 0.053583 |==
4 0.008037 |
5 0.000643 |
6 0.000021 |
As you can see from the numbers, the highest probability is for
throwing one ZETA_1, although there is still a high probability for
throwing zero or two ZETA_1's. If we repeat the binomial calculations
this time assuming hypothesis H1 is true (the probability of ZETA_1 is
2/6), we arrive at the following slightly different table:
f(ZETA_1) Pr ( f(ZETA_1) | H1 ) (sideways bar graph)
----------- ----------------------
0 0.087791 |====
1 0.263374 |=============
2 0.329218 |================
3 0.219478 |==========
4 0.082304 |====
5 0.016460 |=
6 0.001371 |
We can now calculate the probability of making a type I error and of
making a type II error.
Define:
M : f(ZETA_1) = 2
the event that ZETA_1 occurs twice in six rolls
H0 : The probability of ZETA_1 is 1/6
H1 : The probability of ZETA_1 is 2/6
Notation:
Pr(A|B) : The probability of A given B
From Bayes's Theorem we then have:
Pr(M|H0) Pr(H0)
Pr(H0|M) = ---------------
Pr(M)
applying the principle of total probability to Pr(M)
Pr(M|H0) Pr(H0)
Pr(H0|M) = ---------------------------------
Pr(M|H0) Pr(H0) + Pr(M|H1) Pr(H1)
and if we allow both hypothesis H0 and H1 to be equally likely
so that Pr(H0) = Pr(H1) = 1/2
Pr(M|H0)
Pr(H0|M) = -------------------
Pr(M|H0) + Pr(M|H1)
Similarly, we can repeat the above steps to obtain
Pr(M|H1)
Pr(H1|M) = -------------------
Pr(M|H0) + Pr(M|H1)
>From these two formulas and the tables above we can now calculate the
probability of making a type I error and a type II error.
First note that
Pr(M|H0) = probability of rolling ZETA_1 two out of 6 times
when probability of rolling ZETA_1 is 1/6
= 0.200938 (from table above)
Pr(M|H1) = probability of rolling ZETA_1 two out of 6 times
when probability of rolling ZETA_1 = 2/6
= 0.329218 (from table above)
Then from the above formulas and data we have:
Pr(M|H0)
Pr( type I error ) = Pr(H0|M) = -------------------
Pr(M|H0) + Pr(M|H1)
0.200938
= -------------------
0.200938 + 0.329218
= 0.379
Pr(M|H1)
Pr( type II error ) = Pr(H1|M) = -------------------
Pr(M|H0) + Pr(M|H1)
0.329218
= -------------------
0.200938 + 0.329218
= .621
In this case we have a 38% chance of being wrong if we choose to
believe Pr(ZETA_1) = 2/6 and a 62% chance of being wrong if we choose
to believe Pr(ZETA_1) = 1/6. Either way we're likely to be wrong. The
data does not allow us to differentiate between the two possibilities
with much confidence.
Let us repeat the above calculations this time assuming we rolled
ZETA_1 five out of six times. This would lead us to the estimate that
Pr(ZETA_1) = 5/6. Let us create a new hypothesis:
H1 : The probability of ZETA_1 is 5/6
and let us repeat the binomial calculations assuming hypothesis H1 is
true (the probability of ZETA_1 is 5/6), we arrive at the following
table:
f(ZETA_1) Pr ( f(ZETA_1) | H1 ) (sideways bar graph)
----------- ----------------------
0 0.000021 |
1 0.000430 |
2 0.008037 |
3 0.053583 |==
4 0.200938 |==========
5 0.401877 |====================
6 0.334897 |================
Following the same procedure as above:
Define:
M : f(ZETA_1) = 5
the event that ZETA_1 occurs five out of six rolls
H0 : The probability of ZETA_1 is 1/6
H1 : The probability of ZETA_1 is 5/6
Pr(M|H0)
Pr( type I error ) = Pr(H0|M) = -------------------
Pr(M|H0) + Pr(M|H1)
0.000643
= -------------------
0.000643 + 0.401877
= 0.001
Pr(M|H1)
Pr( type II error ) = Pr(H1|M) = -------------------
Pr(M|H0) + Pr(M|H1)
0.401877
= -------------------
0.000643 + 0.401877
= .998
In this case, given a choice between believing the theoretician who
claims Pr(ZETA_1) = 1/6 and believing our own tests which claim
Pr(ZETA_1) = 5/6, if we choose to believe the theoretician we have a
99.8% chance of being wrong, and if we choose to believe our own
experimental results we have a 0.1% chance of being wrong. Clearly the
odds are in favor of our being right. And we only threw the die six
times!
This then is why a sequence consisting of all ones { 1, 1, ..., 1 } is
unlikely to be from a compound experiment in which each outcome is
equally likely. If we hypothesize that each outcome is equally likely
(hypothesis H0), then all elements in our sample space Sn do have an
equal probability. But if we conjure up an alternative hypothesis
about the probability of outcomes (hypothesis H1), then all the
elements in our sample space Sn no longer have an equal probability.
The probability will be higher for some elements and lower for other
elements. Our hypothesis H1 is specifically formulated to give
whatever sequence we actually got the highest probability possible.
If the resulting probability of our sequence occurring under our
original H0 hypothesis is very low, then it's almost for certain that
it will loose out under an alternate hypothesis H1. Thus if all we
want to do is test whether or not to reject our original H0 hypothesis,
all we need to do is see what the probability of our sequence occurring
under our original H0 hypothesis is. If it is very low, we can reject
our H0 hypothesis under the implied knowledge that a better hypothesis
H1 probably exists. We don't actually have to formulate a better
hypothesis H1, we just have to know that we could make one if we cared
to.
Thus from here on out we won't go the whole distance of formulating an
alternative hypothesis H1 and then calculating the probabilities of
type I errors and type II errors. We'll just say if an actual outcome
of a compound experiment has a very low probability of occurring under
our H0 hypothesis (or null hypothesis), then that hypothesis is
probably not true.
3.2 GENERAL TEST PROCEDURE
Our general procedure then for testing the results of a compound
experiment is:
1. Formulate a null hypothesis H0 about the single chance
experiment which was repeated N times to create a sequence of
N values. To test a sequence of supposedly random numbers,
our null hypothesis H0 is that each outcome of the chance
experiment is equally likely, and that each trial of the
chance experiment is independent of all previous trials.
2. Formulate a real valued function g which somehow tests the
null hypothesis H0. To test a sequence of supposedly random
numbers, our function g might be one which counts the number
of occurrences of a particular outcome.
3. Under our null hypothesis H0, mathematically define a sequence
of N random variables:
( X1, X2, ..., Xn )
and apply our function g to the sequence of N random variables
creating a new random variable Z:
Z = g( X1, X2, ..., Xn )
Then determine the probability density function of Z either by
mathematical calculations or by finding a table of the
particular probability density function we're interested in.
4. Take the specific sequence of values supposedly obtained by
performing a chance experiment N times:
( x1, x2, ..., xn )
and apply the function g to obtain a specific value z
z = g( x1, x2, ..., xn )
5. Determine from the probability density function of Z how
likely or unlikely we are to obtain our value of z assuming
our null hypothesis H0 is true. If the probability is very
low, we may reject our hypothesis H0 as being most likely
incorrect.
3.3 THE CHI-SQUARE TEST
We come now to our first real test for a sequence of random numbers.
The test is the Chi-Square (ki-square) test. It tests the assumption
that the probability distribution for a given chance experiment is as
specified. For our die throwing experiment, it tests the probability
that each possible outcome is equally likely with a probability of 1/6.
The formula for our Chi-Square g function is:
Define:
N = number of trials
k = number of possible outcomes of the chance experiment
f(ZETA_i) = number of occurrences of ZETA_i in N trials
E(ZETA_i) = The expected number of occurrences of ZETA_i in N
trials. E(ZETA_i) = N*Pr(ZETA_i).
i=k [ f(ZETA_i) - E(ZETA_i) ]**2
CHISQ = SUM ----------------------------
i=1 E(ZETA_i)
We now go through the steps as defined above:
1. Formulate the null hypothesis H0. Our null hypothesis H0 for
the die throwing experiment is that each of the six possible
outcomes is equally likely with a probability of 1/6.
2. Formulate a real valued function g which somehow tests the null
hypothesis. Our g function is the chi-square test. (This is
not exactly the chi-square function, but for N large enough it
approaches the chi-square distribution.) Note that if the
result CHISQ is zero, then each possible outcome occurred exactly
the expected number of times. On the other hand, if the result
CHISQ is very large, it indicates that one or more of the
possible outcomes of the chance experiment occurred far less or
far more than would be expected if our null hypothesis H0 was
true.
3. Determine the probability density function of Z = g(X1,X2,...Xn).
As mentioned in step 2, our chi-square test function g is not
exactly the chi-square function, but for N large enough it
approaches the chi-square distribution of which tables are
available for in any book on probability theory or any book of
mathematical tables such as the CRC Standard Mathematical Tables
book. The stated rule of thumb is N should be large enough so
that every possible outcome is expected to occur at least 5
times. Thus for our die throwing experiment we should repeat it
at least 30 times, since 30*(1/6) = 5.
Each row of the chi-square distribution table is for a
particular number of "degrees of freedom" which we will
abbreviate to df. The value to choose for df is one less than
the number of possible outcomes for our chance experiment. Thus
for our die throwing experiment where there are 6 possible
outcomes, we have df = 6 - 1 = 5.
Below is the chi-square distribution for 5 degrees of freedom:
p=1% p=5% p=25% p=50% p=75% p=95% p=99%
0.5543 1.1455 2.675 4.351 6.626 11.07 15.09
This is the cumulative chi-square distribution. If we recall
the distribution itself is typically bell shaped, then the
extremes of the bell shaped curve where the probability is very
low corresponds to the extremes of the cumulative distribution.
Thus if we are below the p=5% mark or above the p=95% mark, the
corresponding probability is quite low. On the other hand, if
we are between the p=25% mark and the p=75% mark, the
corresponding probability is high.
(Note that some tables list the probability for (1 - p) instead
of p, so that 0.5543 is listed as p=99% instead of p=1%, and
15.09 is listed as p=1% instead of p=99%. It should be easy to
tell by inspection which probability the table is listing.)
4. We can now apply our chi-square test to a specific outcome of
rolling the die 30 times. Assume we threw the die 30 times and
obtained the following results:
f(ZETA_1) = 7
f(ZETA_2) = 5
f(ZETA_3) = 4
f(ZETA_4) = 6
f(ZETA_5) = 6
f(ZETA_6) = 2
Then our CHISQ test becomes:
(7-5)**2 (5-5)**2 (4-5)**2 (6-5)**2 (6-5)**2 (2-5)**2
CHISQ = -------- + -------- + -------- + -------- + -------- + --------
30*(1/6) 30*(1/6) 30*(1/6) 30*(1/6) 30*(1/6) 30*(1/6)
= 3.2
5. We can now check our result of 3.2 against the chi-square table
above in step 3 and note that it falls between the p=25% mark
and the p=50% mark. This is a high probability region and so we
are not led to believe our null hypothesis H0 is false.
Further notes on chi-square:
Note that not only are overly large values of CHISQ considered
highly improbable, but so are overly small values of CHISQ. Whereas
overly large values of CHISQ lead us to believe the probabilities of
each possible outcome of our chance experiment are not what we
thought them to be, overly small values of CHISQ lead us to believe
the N trials of our chance experiment were not independent.
As an example, if we threw a single die 6 times, there would be 6**6
or 46656 possible outcomes. If we calculate the CHISQ value for
each of the possible values and tabulate them we get the following
table:
CHISQ f(CHISQ)
0 720
2 10800
4 16200
6 9000
8 7200
12 2100
14 450
20 180
30 6
Although a CHISQ value of 0 happens 720 times, a CHISQ value of 2
happens 10800 times, making it much more likely. We can compare our
CHISQ values here with the chi-square distribution table above for 5
degrees of freedom. Our values of CHISQ = 2 through 12 which appear
to be the most likely ones lie above the p=5% mark and below the
p=95% mark on the table, where we would expect them to lie. Note
that this comparison was for N = 6, although we have as a rule of
thumb that N should be at least 30 for this test.
Another note about the chi-square distribution, most tables of the
chi-square distribution go up to 30 degrees of freedom. Above that
the distribution approaches the normal distribution. If your chance
experiment has more than 30 degrees of freedom (number of possible
outcomes - 1), you can use the following formula to convert the CHISQ
distribution for df > 30 into the normal distribution with mean = 0
and variance = 1:
given: CHISQ (and) df
Z = [ SQRT(24*CHISQ - 6*df + 16) - 3*SQRT(2*df) ] / 4
Here Z is a normally distributed variable with mean = 0 and
variance = 1. As an example, let's assume we are simulating
the random choosing of a single card from a deck of 52. There
are 52 possible outcomes. We will need to repeat this
experiment 52*5=260 times so that each possible outcome occurs
about 5 times. Our number of degrees of freedom is 52 - 1 = 51.
We calculate CHISQ and get CHISQ = 30.4 . We now calculate
Z from the above equation given CHISQ = 30.4 and df = 51:
Z = [ SQRT(24*30.4 - 6*51 + 16) - 3*SQRT(2*51) ] / 4
= -2.33
We now look up in a table of the cumulative normal distribution
and find this value of Z corresponds to the first 1% of the
distribution. Thus our CHISQ value is much lower than
expected. We could conclude that something is probably amiss.
One last note about chi-square. If you are simulating the flipping
of a coin, then the number of degrees of freedom is 1 (two possible
outcomes). Our rule of thumb about performing enough experiments so
that each possible outcome occurs at least 5 times is not adequate
in this case. Because the number of degrees of freedom is so low,
we need to compensate for this by increasing our number of trials.
Enough trials so that each outcome is expected to occur at least 75
times is required in this case.
3.4 THE KOLMOGOROV-SMIRNOV TEST
If we repeat the die throwing experiment 30 times, we can perform one
chi-square test on the results. If we repeat the die throwing
experiment 3000 times, we can either perform one chi-square test to get
an idea of the long term probabilities of each possible outcome, or we
can perform 300 chi-square tests for each block of 30 trials to get an
idea of the relative short term fluctuations. In the latter case, we
would expect the distribution of our 300 chi-square tests to approach
the chi-square distribution.
This brings us to our next subject, which is testing a set of data
points to see how closely it matches a given continuous distribution.
The Kolmogorov-Smirnov (KS) test compares a set of nc data points Snc
(e.g. set of a number of chi-square data points) with a given
cumulative distribution function P(x) (e.g. the chi-square cumulative
distribution).
Given a set of nc data points Snc, we can create a cumulative
distribution function Snc(x) from the data points. The function will
be a step function, rising from 0 to 1, with each step size equal to
1/nc. Below is an example which compares a cumulative distribution
step function Snc(x) for 4 data values against a continuous cumulative
distribution P(x) for a function which is uniform between 0 and 1 (the
diagonal line):
|
1.0| +/--
| /|
| +----/--+
| | /
| +-+/
| |/
| +/+
| /|
0.0 |/--+---------------
0 1
The KS test is now very simple. Define D as the maximum vertical
difference between these two functions:
D = max |Snc(x) - P(x)|
And define DQ as:
DQ = SQRT(nc)*D
where nc is the number of data points.
We can then go directly to a table for the KS test (similar to the
chi-square table), and read off where in the cumulative distribution
our DQ value lies, remembering that values below p=5% only occur 5% of
the time, and values above p=95% only occur 5% of the time. Tables for
this KS test are not as easy to find as tables for the chi-square test.
Fortunately, as in the chi-square case, if our number of data values nc
in set Snc is greater than 30, the distribution approaches something
which can be easily calculated as follows:
Given: D = max |Snc(x) - P(x)|
nc > 30 (number of data points)
1
p = 1 - --------------------------------
2*D**2 2*D 1
exp( ------ + ---------- + ----- )
1 3*SQRT(nc) 18*nc
----------------------------------------------------------------------
4.0 ANALYSIS OF A SEQUENCE OF RANDOM NUMBERS
Enough theory for now. Let's turn to the analysis of a sequence of
random numbers. There's nothing "random" about any specific sequence of
numbers. But we can analyze the sequence of numbers to see if we can
find any sufficient reason to doubt the hypothesis that the numbers
could have come from multiple independent trials of a chance experiment.
4.1 ONE-DIMENSIONAL TESTS
Our first test then is to see if each possible outcome appears to be
equally likely. Imagine we have a number of bins, and we toss balls at
the bins. After throwing a certain number of balls at the bins, we
count up how many balls fell in each bin. We then apply the chi-square
test to determine if the balls appear to have an equal probability of
falling in each bin.
The MTH$RANDOM and RANDU generators can be set up to generate numbers
between 1 and NBINS as follows:
I = INT( NBINS * RAN(SEED) ) + 1 !MTH$RANDOM
I = INT( NBINS * FOR$IRAN(SEED1,SEED2) ) + 1 !RANDU
Arbitrarily choosing NBINS = 30 (number of bins), SEED = 1 (initial SEED
value), and BALLS = 300, the chi-square test was performed 10 times on
the MTH$RANDOM generator with the following results:
30 BINS, 300 BALLS, 10 CHI-SQUARE TESTS, MTH$RANDOM
BALLS= 300 CHISQ= 35.1999969 PROB= 0.8019531
BALLS= 300 CHISQ= 22.8000031 PROB= 0.2143845
BALLS= 300 CHISQ= 36.7999992 PROB= 0.8485908
BALLS= 300 CHISQ= 19.8000011 PROB= 0.1009650
BALLS= 300 CHISQ= 48.7999992 PROB= 0.9878765
BALLS= 300 CHISQ= 29.3999996 PROB= 0.5556139
BALLS= 300 CHISQ= 22.8000011 PROB= 0.2143840
BALLS= 300 CHISQ= 36.5999985 PROB= 0.8432655
BALLS= 300 CHISQ= 29.3999996 PROB= 0.5556139
BALLS= 300 CHISQ= 18.6000004 PROB= 0.0688844
Here CHISQ is the chi-square value, and PROB is the percentage of random
sequences (with 30 - 1 = 29 degrees of freedom) which would have a lower
CHISQ value. An extremely high value of PROB, say above 0.990, would
indicate that some bins had a much higher than expected number of balls
in them while other bins have a much lower than expected number of balls
in them. In general, an extremely high value of PROB indicates that the
experimental results do not reinforce the hypothesis that each bin
had an equally likely chance of getting a ball. So far all of the PROB
values appear reasonable.
We can further analyze the results of the 10 chi-square tests above by
noting that the distribution of the 10 CHISQ values should approximate
the chi-square distribution (for 30 - 1 = 29 degrees of freedom). We
can apply the KS test to the 10 CHISQ values to see how well their
distribution approximates the chi-square distribution for d.f. = 29.
KS D= 0.2019531 PROB= 0.8093885
Here D is the KS D value which is the maximum distance between the
approximate and the real chi-square function, and PROB is the
percentage of times we would have a larger D value, meaning the
approximation could be worse. A very low value of PROB here means it
couldn't get much worse. Anything above a very low value means our
distribution bears at least some resemblance to the function we're
comparing it to. Again this result appears reasonable.
The following values were obtained when the above chi-square test was
applied 10 times to the RANDU generator:
30 BINS, 300 BALLS, 10 CHI-SQUARE TESTS, RANDU
BALLS= 300 CHISQ= 31.3999996 PROB= 0.6531839
BALLS= 300 CHISQ= 60.7999992 PROB= 0.9995093
BALLS= 300 CHISQ= 33.4000015 PROB= 0.7381005
BALLS= 300 CHISQ= 24.3999977 PROB= 0.2910264
BALLS= 300 CHISQ= 20.8000011 PROB= 0.1336691
BALLS= 300 CHISQ= 16.6000023 PROB= 0.0319657
BALLS= 300 CHISQ= 32.0000000 PROB= 0.6801265
BALLS= 300 CHISQ= 30.2000027 PROB= 0.5959312
BALLS= 300 CHISQ= 31.2000027 PROB= 0.6439425
BALLS= 300 CHISQ= 45.5999985 PROB= 0.9742958
And a KS test on the above 10 CHISQ values:
KS D= 0.2959312 PROB= 0.3452127
Again there is nothing as yet unusual about either the MTH$RANDOM or
RANDU generators.
We can also test the MTH$RANDOM and RANDU generators using the KS test
directly to see how well the output fits the hypothesized uniform
distribution.
KS TEST ON MTH$RANDOM (Start with SEED = 1)
KS NPOINTS= 10 D= 0.2142200 PROB= 0.7484154
KS NPOINTS= 100 D= 0.0944867 PROB= 0.3338285
KS NPOINTS= 1000 D= 0.0314864 PROB= 0.2746516
KS NPOINTS= 10000 D= 0.0079555 PROB= 0.5514056
KS NPOINTS= 100000 D= 0.0017763 PROB= 0.9105781
KS NPOINTS= 1000000 D= 0.0009270 PROB= 0.3565834
KS TEST ON RANDU (Start with SEED = 1)
KS NPOINTS= 10 D= 0.6555050 PROB= 0.0003705
KS NPOINTS= 100 D= 0.1326773 PROB= 0.0591586
KS NPOINTS= 1000 D= 0.0337385 PROB= 0.2050490
KS NPOINTS= 10000 D= 0.0063568 PROB= 0.8138239
KS NPOINTS= 100000 D= 0.0042999 PROB= 0.0495556
KS NPOINTS= 1000000 D= 0.0007990 PROB= 0.5457213
Note the only thing questionable here is the extremely low value of
PROB on RANDU for the first 10 points. This indicates that perhaps a
seed value of 1 for the RANDU generator is not a good place to start
it. Indeed, the first 10 numbers from RANDU starting with SEED = 1
are:
RANDU first 10 numbers starting with seed = 1
1 0.0001831
2 0.0032959
3 0.0444950
4 0.5339386
5 0.0068024
6 0.8734164
7 0.1705012
8 0.3222913
9 0.9906484
10 0.7260775
Notice that half of the first 10 values are either below 0.006 or above
0.99. This alone does not indicate a problem with the RANDU generator,
it does indicate that a SEED value of 1 is a poor choice to start it
with. Note that the MTH$RANDOM generator does not have this problem when
started with SEED = 1.
4.2 TWO-DIMENSIONAL TESTS
Our original hypothesis about our sequence of supposedly random numbers
was that it could have come from a repeated chance experiment in which
each possible outcome was equally likely and each test of the chance
experiment was independent of all previous tests. So far we have been
testing the sequence to see if each possible outcome appears to be
equally likely. The remainder of our tests now focus on checking the
hypothesis that each outcome is independent of all previous outcomes.
If we throw a die and get a 5, what is the probability that the next
throw of the die will be a 3? If the probability is still 1/6, we say
the next outcome is independent of the previous outcome. i.e. knowledge
of the previous outcome does not give us any help in predicting the next
outcome.
The mathematical definition of statistical independence between two
random variables can be written as:
Pr( X(i+1) | X(i) ) = Pr( X(i+1) )
The serial test in two dimensions checks to see if there is any
correlation between two consecutive outcomes of the random number
generator. This is a chi-square test performed on pairs of random
numbers. For the die throwing experiment, we would construct a total of
36 bins numbered BIN(1,1) - BIN(6,6). One bin for each possible pair of
numbers we may get. We then throw the die twice to get an ordered pair
of numbers, and drop a ball in the corresponding bin. We throw the die
enough times so that each bin is expected to get an average of at least
5 balls. Then we perform a chi-square test on the results to see if the
balls appear to be more or less evenly distributed among all the bins.
This serial test for pairs of outcomes was done on both the MTH$RANDOM
and RANDU generators with the following results:
SERIAL 2D BINS=30x30 BALLS = 9000 CHI-SQUARE TESTS = 10 MTH$RANDOM
BALLS= 9000 CHISQ= 895.7998657 PROB= 0.4761399
BALLS= 9000 CHISQ= 945.2001343 PROB= 0.8615244
BALLS= 9000 CHISQ= 883.6000366 PROB= 0.3633031
BALLS= 9000 CHISQ= 905.0000000 PROB= 0.5624363
BALLS= 9000 CHISQ= 902.3989868 PROB= 0.5382197
BALLS= 9000 CHISQ= 911.8001709 PROB= 0.6241364
BALLS= 9000 CHISQ= 932.4005737 PROB= 0.7863315
BALLS= 9000 CHISQ= 865.4000854 PROB= 0.2157318
BALLS= 9000 CHISQ= 909.5996704 PROB= 0.6043593
BALLS= 9000 CHISQ= 901.7994385 PROB= 0.5325246
KS test to see how well the 10 CHISQ values fit the chi-square
distribution:
KS D= 0.2667681 PROB= 0.4751010
SERIAL 2D BINS=30x30 BALLS = 9000 CHI-SQUARE TESTS = 10 RANDU
BALLS= 9000 CHISQ= 848.7999268 PROB= 0.1168594
BALLS= 9000 CHISQ= 891.7999878 PROB= 0.4386667
BALLS= 9000 CHISQ= 874.1998291 PROB= 0.2827876
BALLS= 9000 CHISQ= 884.1997681 PROB= 0.3687753
BALLS= 9000 CHISQ= 799.5999756 PROB= 0.0077433
BALLS= 9000 CHISQ= 841.8005981 PROB= 0.0864621
BALLS= 9000 CHISQ= 867.8007202 PROB= 0.2330545
BALLS= 9000 CHISQ= 892.6010132 PROB= 0.4461077
BALLS= 9000 CHISQ= 952.4000854 PROB= 0.8945085
BALLS= 9000 CHISQ= 921.5999146 PROB= 0.7069101
KS test to see how well the 10 CHISQ values fit the chi-square
distribution:
KS D= 0.3631590 PROB= 0.1430003
So far there is nothing blatantly wrong with either the MTH$RANDOM or
RANDU generators. All the chi-square values are within reason, and the
distribution of the chi-square values appears to fit the chi-square
distribution.
4.3 THREE-DIMENSIONAL TESTS
We can now expand the serial test to triples of numbers. Suppose you
are throwing a die, and you know that the results of the last two throws
were a 2 and then a 5. What is the probability that you will now throw
a 3? It should still be 1/6, because the outcome of this throw should
not be affected by the outcome of the previous two throws. We can test
for this by constructing 6**3 = 216 bins labeled BIN(1,1,1) -
BIN(6,6,6), one bin for each possible triple of outcomes. We then
throw the die three times to get a triple of numbers, and then drop
a ball in the corresponding bin. We throw the die enough times so that
each bin is expected to get an average of at least 5 balls. Then we
perform a chi-square test on the results to see if the balls appear to
be more or less evenly distributed among the bins.
This 3-dimensional test extends the test for statistical independence.
If we assume the first two variables are independent, then this test
checks to see if the third variable is also independent. The
mathematical definition of statistical independence between two random
variables can be written as:
Pr( X(i+2) | X(i+1),X(i) ) = Pr( X(i+2) )
This serial test for triples of outcomes was done on both the
MTH$RANDOM and RANDU generators. The generators were set up to simulate
the throwing of a 30 sided die. A total of 270000 balls were tossed so
that an average of 10 balls were expected to fall in each bin. A
chi-square test was then performed on the results. This was repeated a
total of 10 times, and a KS test was done on the resulting chi-square
values to see how well they fit the chi-square distribution. Below are
the results:
SERIAL 3D BINS=30x30x30 BALLS = 270000 MTH$RANDOM
BALLS= 270000 CHISQ= 27233.4375000 PROB= 0.8438070
BALLS= 270000 CHISQ= 26732.8027344 PROB= 0.1262939
BALLS= 270000 CHISQ= 26866.4550781 PROB= 0.2845250
BALLS= 270000 CHISQ= 26765.3710938 PROB= 0.1561499
BALLS= 270000 CHISQ= 26650.6250000 PROB= 0.0659529
BALLS= 270000 CHISQ= 26665.5117188 PROB= 0.0751096
BALLS= 270000 CHISQ= 27165.1523438 PROB= 0.7621238
BALLS= 270000 CHISQ= 26861.5625000 PROB= 0.2786521
BALLS= 270000 CHISQ= 27002.1171875 PROB= 0.5027421
BALLS= 270000 CHISQ= 27090.8613281 PROB= 0.6547577
KS test to see how well the 10 CHISQ values fit the chi-square
distribution:
KS D= 0.3162388 PROB= 0.2699622
SERIAL 3D BINS=30x30x30 BALLS = 270000 RANDU
BALLS= 270000 CHISQ= 454159.3750000 PROB= 1.0000000
BALLS= 270000 CHISQ= 452698.1250000 PROB= 1.0000000
BALLS= 270000 CHISQ= 455055.3437500 PROB= 1.0000000
BALLS= 270000 CHISQ= 455938.2812500 PROB= 1.0000000
BALLS= 270000 CHISQ= 454747.7187500 PROB= 1.0000000
BALLS= 270000 CHISQ= 453017.8125000 PROB= 1.0000000
BALLS= 270000 CHISQ= 453818.7812500 PROB= 1.0000000
BALLS= 270000 CHISQ= 454553.5937500 PROB= 1.0000000
BALLS= 270000 CHISQ= 453205.6875000 PROB= 1.0000000
BALLS= 270000 CHISQ= 454122.0625000 PROB= 1.0000000
KS test to see how well the 10 CHISQ values fit the chi-square
distribution:
KS D= 1.0000000 PROB= 0.0000000
It is immediately obvious that the RANDU generator has failed miserably
on the serial test for triples! The chi-square value is way too high
on every one of the 10 trials.
4.4 HIGHER DIMENSIONAL TESTS
We can easily expand the chi-square tests for higher dimensions. One
thing to note though, the number of bins increases exponentially as the
number of dimensions increases. Thus there is a practical limit as to
how high a dimension we can carry this test. In section 6 (Ratings for
Various Generators) we analyze 5 popular random number generators by
using the chi-square test for dimensions 1 through 8. Higher
dimensions are not as serious since such a large number of random
numbers are required before some nonrandomness can be detected.
Before we use these tests to analyze some random number generators we
should first take a look at what's inside these generators and how to
properly use them.
---------------------------------------------------------------------------
5.0 LINEAR CONGRUENTIAL RANDOM NUMBER GENERATORS
Both the MTH$RANDOM and the RANDU algorithms are specific instances of
linear congruential generators. The general form of a linear
congruential generator is:
SEED = (A * SEED + C) mod M
X = SEED/M
The user supplies an initial integer value for SEED. On each call to
the random number generator a new value of SEED is calculated by taking
the current value of SEED, multiplying by A, adding C, and taking the
remainder MOD M of the result. The new value of SEED is an integer
between 0 and M-1. This new value of SEED is then converted into a
floating point value between 0 inclusive and 1 exclusive by dividing
SEED by M. The result X is returned as a floating point random number
in the range [0,1).
First note some obvious properties of sequences generated by this method:
1. The modulus M is an upper bound on the number of different values
SEED can take on.
2. If the new value of SEED is the same as one we had before,
the sequence will begin to repeat itself in a cyclic
manner.
3. All sequences generated by a linear congruential formula will
eventually enter a cycle which repeats itself endlessly.
A good linear congruential formula will generate a long sequence of
numbers before repeating itself. The maximum length a sequence can
possibly be before repeating itself is M, the modulus. This is because
there are only M different possible values SEED can take on. A linear
congruential formula will generate a sequence of maximum length M if and
only if the following three conditions are met (See Knuth for a proof of
this theorem):
i) C is relatively prime to M;
ii) B = (A - 1) is a multiple of P, for every prime P dividing M;
iii) B = (A - 1) is a multiple of 4, if M is a multiple of 4.
5.1 THE MTH$RANDOM ALGORITHM:
We now take a look at the MTH$RANDOM function. The algorithm MTH$RANDOM
uses to generate successive random numbers is:
SEED = (69069*SEED + 1) MOD 2**32
X = SEED/2**32
This is a linear congruential generator with:
A = 69069
C = 1
M = 2**32
Note first that MTH$RANDOM satisfies the three necessary and sufficient
conditions above which insure that the sequence it generates is of maximum
length:
i) 1 is relatively prime to 2**32
since 1 is relatively prime to all numbers.
ii) 69068 is a multiple of 2, which is the only prime dividing 2**32.
iii) 69068 is a multiple of 4.
Thus the sequence generated by MTH$RANDOM is of maximal length,
generating all 2**32 possible values of SEED before repeating itself.
Note for the MTH$RANDOM function if SEED is initially an ODD value
then the new value of SEED will always be an even value. And if SEED
is an EVEN value, then the new value of SEED will be an ODD value.
Thus if the algorithm is repeatedly called, the value of SEED will
alternate between EVEN and ODD values.
5.2 THE RANDU ALGORITHM
Now let us look at the RANDU function. The algorithm RANDU uses to
generate successive random numbers is:
SEED = (65539*SEED) MOD 2**31
X = SEED/2**31
This is a linear congruential generator with:
A = 65539
C = 0
M = 2**31
Note first that the modulus M for RANDU is 2**31 whereas the modulus M
for MTH$RANDOM is 2**32. The maximum sequence length for RANDU is thus
at least 1/2 of what it is for MTH$RANDOM. Note also that if SEED is
initially an odd value, the new SEED generated will also be an odd
value. Similarly, if SEED is initially an even value, the new SEED
generated will also be an even value. Thus there are at least two
separate disjoint cycles for the RANDU generator, depending upon whether
you start with an even SEED value or and odd SEED value.
Actually, the situation is even worse than that. For odd SEED values
there are two separate disjoint cycles, one generated by the SEED value
1, and one generated by the SEED value 5. Although we can generate each
full cycle by initializing SEED to any value the cycle takes on, it is
convenient to refer to the smallest SEED value the cycle takes on as the
one which generates the cycle. This helps us differentiate between
different cycles. We will denote the sequence generated by the SEED
value 1 as <1>, and the sequence generated by the SEED value 5 as <5>.
The cycles <1> and <5> each contain 536,870,912 values. Together they
account for all of the (2**31)/2 possible odd SEED values. The
following table summarizes the cycles for odd values of SEED:
TABLE 5.2.1 RANDU WITH ODD VALUES OF SEED
CYCLE LENGTH OF CYCLE
<1> 536,870,912
<5> 536,870,912
Even values of SEED do not have it so nice. There are 30 different
disjoint cycles using even SEED values. Some of them are quite short:
TABLE 5.2.2 RANDU WITH EVEN VALUES OF SEED
CYCLE LENGTH OF CYCLE
<2> 268435456
<4> 134217728
<8> 67108864
<10> 268435456
<16> 33554432
<20> 134217728
<32> 16777216
<40> 67108864
<64> 8388608
<80> 33554432
<128> 4194304
<160> 16777216
<256> 2097152
<320> 8388608
<512> 1048576
<640> 4194304
<1024> 524288
<1280> 2097152
<2048> 262144
<2560> 1048576
<4096> 131072
<5120> 524288
<8192> 65536
<10240> 262144
<16384> 32768
<20480> 131072
<32768> 16384
<40960> 65536
<81920> 32768
<163840> 16384
---------
Total: 1073709056
There are a total of (2**31)/2 = 1,073,741,824 possible even SEED
values. So far we've accounted for 1,073,709,056 of them. The
remaining 32768 SEED values are ones for which the 31 bit binary
representation of them has the lower 16 bits set to 0. These SEED
values are treated by RANDU as if the SEED value were 1, and they
result in the cycle <1>.
We can see from this that even values of SEED should be avoided when
using RANDU since they can give us very short cycles.
-----------------------------------------------------------------------
5.3 STARTING THE MTH$RANDOM GENERATOR
The initial SEED value for the MTH$RANDOM generator can be any integer
value between 0 and (2**32)-1 = 4294967295. If we could "randomly"
choose a number in this range and insure that each number had an equal
probability of being chosen, we'd be set. However, human nature being
what it is, simply picking a number out of the air usually does not
satisfy the criterion for randomness, especially when there are over
four billion numbers to choose from.
For many applications it doesn't matter what the initial few numbers
from the random number generator are, and ANY initial seed value will
suffice. Note also some myths about choosing the initial seed value
for the MTH$RANDOM generator. The VAX FORTRAN manual recommends using
an odd value. This piece of advice is obviously left over from the
RANDU days. The obsolete RANDU generator does require an odd seed
value, but for the MTH$RANDOM generator the seed value alternates
between even and odd values at each iteration, so there is no advantage
to starting with an odd value.
Another incorrect comment about the MTH$RANDOM generator can be found
in the VMS source code comments. There it indicates that the
MTH$RANDOM generator might do poorly in a 3-D test for randomness, but
actually MTH$RANDOM does quite well in a 3-D test. It's the RANDU
generator which does poorly in a 3-D test.
More important than starting the MTH$RANDOM generator to get
one random sequence is the problem of restarting the generator to get
several different sequences. You may wish to run a simulation
several times and use a different random sequence each time.
One simple method would be to use consecutive seed values for each run.
For example, you could generate ten different random sequences using
initial seed values of 1 through 10. If we compare side by side these
ten random sequences the first thing we'll note is the first value of
each sequence is nearly identical. This is because there is a very
strong correlation between the initial SEED value and the first value
output from the generator. The simple relationship between SEED and
the first output value is approximately:
SEED
X ~= ------- mod 1
62183.7
where mod 1 here indicates when SEED/62183.7 becomes greater than 1 we
wrap around and start counting from 0 again.
Thus if SEED changes by 10 the first value changes by only (10/62183.7)
~= 0.00016 . Note that this correlation between the initial SEED value
and the first value output from the MTH$RANDOM generator will happen no
matter what the initial value of SEED is. Large initial values of SEED
are no better than small initial values of SEED.
Besides the initial expected correlation of the first value output from
the generator, there are also unexpected correlations between other
sequence elements. For the example above, not only are the first
values of all 10 sequences nearly identical, there is also a strong
correlation for elements 6, 13, 39, 41, 54, 60, 65, 72, 82, and 84.
You can play with program CORELA.FOR to see these correlations for
yourself. Program CORELA generates the first 100 values for ten
different sequences, and then compares them. Most of the elements are
quite random, but once in a while you come across an element which is
nearly identical in all ten sequences. File CORELA.DOC shows some
sample outputs with comments.
It seems to be not all that easy to get away from this problem of
occasional correlation between different sequences. If instead of
using the seed values {1,2,3,...,10} you used the seed values
{1,1001,2001,3001,...,9001}, you would still get occasional elements
which are strongly correlated, they just occur in slightly different
places. In this example they occur at elements 1, 11, 17, 32, and 80.
For many applications this occasional correlation may not be important.
Each sequence when taken as a whole is quite different from the other
ones.
One method for dealing with the correlation between the initial SEED
value and the first number output from the random number generator is
to simply throw away the first number. The second number output from
the generator has a much weaker correlation with the initial SEED
value. It is much more "random" in that very small changes in the
initial SEED value produce very large and somewhat unpredictable
changes in it's value. And since you're throwing away numbers you
might as well throw away the second number too. The third number
output from the MTH$RANDOM generator is quite random compared to the
initial SEED value.
There is another way to view this operation of throwing away the first
two values. We need to come up with a good "random" initial SEED value
to start the generator with. So we can use our random number generator
to help us generate a random SEED value to start our random number
generator with! We select an initial nonrandom SEED value, then run
our random number generator a few cycles to generate a random SEED
value. We then restart our random number generator with the new random
SEED value, and the output is then a properly initialized random
sequence. Note that this is the same as starting with a nonrandom SEED
value and throwing away the first few numbers output from the
generator.
Another method to deal with correlation is to use the shuffling
technique to be described in section 8 (Shuffling). This technique
essentially removes all problems of correlation and also greatly
improves the randomness of the numbers. If shuffling is used, then ANY
set of initial seed values may be used with confidence including
something as simple as {1,2,...10}.
5.4 STARTING THE RANDU GENERATOR
The RANDU generator should only be started with an odd SEED value, since
even initial SEED values can result in sequences with very short
periods. Note that for RANDU, if the initial SEED value is odd, then
all subsequent SEED values will be odd.
The rest of the problems mentioned for MTH$RANDOM also apply to RANDU,
except that instead of throwing away the first two random numbers
output you should instead throw away the first ten numbers output.
Actually you shouldn't be using the RANDU generator at all because of
it's very poor 3-D performance.
5.5 STARTING THE ANSI C GENERATOR
The ANSI C generator is very much like the MTH$RANDOM generator. All
comments for the MTH$RANDOM generator also apply to the ANSI C
generator. See the comments above for the MTH$RANDOM generator.
---------------------------------------------------------------------------
6.0 RATINGS FOR VARIOUS GENERATORS
Here are the results of testing 5 different commonly found generators.
The generators tested are:
1. MTH$RANDOM
used by VAX FORTRAN and VAX BASIC
2. RANDU
obsolete but still found on some systems
3. C and ANSI C
The standard rand() function for C and ANSI C as defined by
Brian W. Kernighan and Dennis M. Ritchie in their book "The C
programming language". Used by VAX C.
4. Microsoft C
The rand() function in the Microsoft C library v4.0.
5. Turbo Pascal
The random function in Turbo Pascal v6.0
6.1 REVIEW OF TESTS PERFORMED
We'll review once more the tests which were performed on each
generator.
6.1.1 THE 1-D TEST
The 1-D test is the frequency test. Imagine a number line stretching
from 0 to 1. We will use our random number generator to plot random
points on this line. First we divide the line into a number of
segments or "bins". Below is an example of dividing the unit line into
4 segments or bins:
1 2 3 4
|---------|---------|---------|---------|
0 .25 .50 .75 1.0
We then see how randomly the random number generator fills our bins.
If the bins are filled too unevenly, the Chi-Square test will give a
value that's too high indicating the points do not appear random.
For example, if all the points landed in bin 1 and no points landed in
any of the other bins, we would conclude that the generator did not
appear to be random. On the other hand, if the bins are filled too
evenly, the Chi-Square test will give a value that's too low, also
indicating the points do not appear random. For example, if after
plotting 1000 points each of the four bins had exactly the same number
of points in them, we would conclude that the generator was not random
enough. The values given for the 1-D test is the approximate number
of bins above which the random number generator fails the Chi-Square
test in one dimension.
6.1.2 THE 2-D TEST
The 2-D test is the serial test. We can imagine a unit square
stretching from (0,0) at one tip to (1,1) at the other tip. We will
use the random number generator to plot points on this unit square. To
plot a point we generate two random numbers and use them as the
coordinates. But first we divide the unit square into a number of
smaller squares or "bins", in a gridlike pattern. Below is an example
of dividing up the unit square into a total of 16 bins with 4 bins on a
side:
1.0|---------|---------|---------|---------|
| | | | |
| | | | |
.75|---------|---------|---------|---------|
| | | | |
| | | | |
.50|---------|---------|---------|---------|
| | | | |
| | | | |
.25|---------|---------|---------|---------|
| | | | |
| | | | |
|---------|---------|---------|---------|
0 .25 .50 .75 1.0
We then see how randomly the random number generator fills our bins,
the same as in the 1-D test. The value given for the 2-D test is the
approximate number of bins per side or bins per dimension (bpd) above
which the random number generator fails the Chi-Square test in two
dimensions.
If we were to look at the unit square after we had plotted a large
number of points, we would notice that rather than appearing random,
the points all lined up on parallel lines. All linear congruential
random number generators have this flaw. We will have more to say
about this in section 7 (The Spectral Test).
6.1.3 THE 3-D TEST
The 3-D test is the logical extension of the 2-D test and the 1-D test.
Imagine a unit cube stretching from (0,0,0) at one tip to (1,1,1) at
the other tip. We will use our random number generator to plot points
on this unit cube. To plot a point we generate three random numbers
and use them as the coordinates. We divide the unit cube into a number
of smaller cubes or "bins", and then see how randomly the random number
generator fills our bins, the same as we do for the 2-D test and 1-D
test. The value given for the 3-D test is the approximate number of
bins per side above which the random number generator fails the
Chi-Square test in three dimensions.
If we were to look at the unit cube after we had plotted a large
number of points, we would notice that rather than appearing random,
the points all lined up on parallel planes. All linear congruential
random number generators have this flaw. The infamous RANDU
generator has only 16 parallel planes covering the unit cube, hence it
performs extremely poorly on the 3-D test. We will have more to say
about this in section 7 (The Spectral Test).
6.1.4 THE 4-D TEST
The 4-D test and higher dimensional tests are again the logical
extension of the 3-D test, 2-D test, and 1-D test. For an n-D or
n-dimensional test, we imagine an n-dimensional hypercube, divided in a
gridlike manner into smaller hypercubes or "bins". We plot points in
this n-dimensional hypercube by generating n random numbers and using
them as the coordinates of a point. We then perform a chi-square test
to see how evenly or unevenly our bins get filled with points. The
value given for any n-D test is the approximate number of bins per
dimension (bpd) above which the random number generator failed the
chi-square test. (The total number of bins in an N-dimensional
hypercube would be the number of bins per dimension raised to the Nth
power.)
6.2 TEST RESULTS
Each generator was tested for dimensions 1 through 8. Here now are the
results:
6.2.1 MTH$RANDOM
DEFINITION:
SEED = (69069*SEED + 1) mod 2**32
X = SEED/2**32
returns real in range [0,1)
RATING:
Fails 1-D above 350,000 bpd (bins per dimension)
Fails 2-D above 600 bpd
Fails 3-D above 100 bpd
Fails 4-D above 27 bpd
Fails 5-D above 15 bpd
Fails 6-D above 9 bpd
Fails 7-D above 7 bpd
Fails 8-D above 4 bpd
Comments:
This generator is also used by the VAX FORTRAN intrinsic
function RAN, and by the VAX BASIC function RND.
6.2.2 RANDU
DEFINITION:
SEED = (65539*SEED) mod 2**31
X = SEED/2**31
returns real in range [0,1)
RATING:
Fails 1-D above 200,000 bpd
Fails 2-D above 400 bpd
Fails 3-D above 3 bpd
Fails 4-D above 6 bpd
Fails 5-D above 6 bpd
Fails 6-D above 4 bpd
Fails 7-D above 3 bpd
Fails 8-D above 3 bpd
Comments:
Note the extremely poor performance for dimensions 3 and
above. This generator is obsolete.
6.2.3 ANSI C ( rand() )
DEFINITION:
SEED = (1103515245*SEED + 12345) mod 2**31
X = SEED
returns integer in range [0, 2**31)
RATING:
Fails 1-D above 500,000 bpd
Fails 2-D above 600 bpd
Fails 3-D above 80 bpd
Fails 4-D above 21 bpd
Fails 5-D above 15 bpd
Fails 6-D above 8 bpd
Fails 7-D above 7 bpd
Fails 8-D above 5 bpd
6.2.4 Microsoft C v4.0 rand()
DEFINITION:
SEED = (214013*SEED + 2531011) mod 2**31
X = int( SEED/2**16 )
returns integer in range [0, 2**15)
Fails 1-D above 3000 bpd
Fails 2-D above 700 bpd
Fails 3-D above 84 bpd
Fails 4-D above 16 bpd
Fails 5-D above 15 bpd
Fails 6-D above 7 bpd
Fails 7-D above 6 bpd
Fails 8-D above 4 bpd
Comments:
The poor performance for 1-d is partly due to truncating
the lower 16 bits of SEED before returning X. X is
returned as a 15 bit positive signed integer.
6.2.5 Turbo Pascal v6.0 (random)
DEFINITION:
SEED = (134775813*SEED + 1) mod 2**32
X = int(SEED/2**16)
returns integer in range [0, 2**16)
Fails 1-D above 6000 bpd
Fails 2-D above 700 bpd
Fails 3-D above 87 bpd
Fails 4-D above 27 bpd
Fails 5-D above 13 bpd
Fails 6-D above 9 bpd
Fails 7-D above 7 bpd
Fails 8-D above 4 bpd
Comments:
The poor performance for 1-d is partly due to truncating
the lower 16 bits of SEED before returning X. X is
returned as a 16 bit unsigned integer.
----------------------------------------------------------------------------
7.0 THE SPECTRAL TEST
As stated earlier, for a linear congruential random number generator,
a 2-D plot will show that all the points lie on parallel lines, and a
3-D plot will show that all the points lie on parallel planes. In
general, an n-D plot would show that all the points lie on (n-1)
dimensional hyperplanes.
One test of a linear congruential random number generator is to see how
far apart these parallel planes are for each dimension. The spectral
test calculates the distances between these parallel planes for a
linear congruential random number generator given the modulus M and the
multiplier A of the generator. A good random number generator would
need to have the planes closely spaced for each dimension so that the
individual planes do not become apparent early on. RANDU, having only
16 widely spaced parallel planes in the 3-D unit cube, fails the
spectral test quite horribly.
We will not discuss further the spectral test here. There is a working
program SPECTRAL.FOR which you can play with. Further information on
the spectral test can be found in the book by Knuth.
---------------------------------------------------------------------------
8.0 SHUFFLING
A simple way to greatly improve any random number generator is to
shuffle the output. Start with an array of dimension around 100. (The
exact size of array is not important.) This is the array that will be
used to shuffle the output from the random number generator.
Initialize the array by filling it with random numbers from your random
number generator. Then when the program wants a random number,
randomly choose one from the array and output it to the program. Then
replace the number chosen in the array with a new random number from
the random number generator.
Note that this shuffling method uses two numbers from the random number
generator for each random number output to the calling program. The
first random number is used to pick an element from the array of random
numbers, and the second random number is used to replace the number
chosen.
The improvement in randomness is quite dramatic. We tried this
shuffling technique with the ANSI C random number generator. The
results are given below first without shuffling and then with shuffling
using an array size of 128 (a convenient power of 2):
ANSI C rand()
WITHOUT SHUFFLING: WITH SHUFFLING:
1-D Fails above 500,000 bpd Fails above 400,000 bpd
2-D Fails above 600 bpd Fails above 3100 bpd
3-D Fails above 80 bpd Fails above 210 bpd
4-D Fails above 21 bpd Fails above 55 bpd
5-D Fails above 15 bpd Fails above 24 bpd
6-D Fails above 8 bpd Fails above 14 bpd
7-D Fails above 7 bpd Fails above 9 bpd
8-D Fails above 5 bpd Fails above 7 bpd
Notice the dramatic improvement in performance when shuffling is
used. This shuffling technique works so well it may even fix up the
defective RANDU generator to the point where it is acceptable. (Only a
lack of time prevented this test from being performed.) This shuffling
technique is highly recommended for anyone who is even remotely
concerned about the quality of their random number generator.
Shuffling also takes care of the problem of points lying "mainly in the
planes" as was discussed earlier. Points plotted on a unit square
appear random rather than lying on parallel lines, and points plotted
on a unit cube appear random rather than lying on parallel planes. In
general points plotted on a unit n-dimensional hypercube no longer all
lie on parallel (n-1) dimensional hypercubes.
This shuffling technique is highly recommended.
---------------------------------------------------------------------------
9.0 DES AS A RANDOM NUMBER GENERATOR
The Data Encryption Standard (DES) is described by the reference "Data
Encryption Standard", 1977 January 15, Federal Information Processing
Standards Publication, number 46 (Washington: U.S. Department of
Commerce, National Bureau of Standards). This is the American National
Standard Data Encryption Algorithm X3.92-1981.
This algorithm was always intended to be implemented in hardware
because it is so complicated a software implementation would be too
slow for any practical use. Nevertheless software implementations
exist and a VAX MACRO implementation is provided in file DES.MAR
DES can be set up as a random number generator by feeding the output
back into the input the same as with a liner congruential generator.
The DES algorithm returns a 64 bit quadword. For the tests below the
low 32 bits was taken as one random unsigned longword and the high 32
bits was taken as another random unsigned longword. Thus we got two
random numbers for each iteration. Still it was incredibly slow.
Here are the results of testing DES as a 32 bit random number generator:
1-D PASS 1,000,000 bpd, fails above
2-D PASS 1100 bpd, fails above
3-D PASS 90+bpd, highest tested due to slowness of DES
4-D PASS 33+bpd, highest tested due to slowness of DES
5-D PASS 16+bpd, highest tested due to slowness of DES
6-D PASS 9+bpd, highest tested due to slowness of DES
7-D PASS 7+bpd, highest tested due to slowness of DES
8-D PASS 4 bpd, fails above
Note that DES is extremely slow. DES takes an hour to generate the
same number of random variables a shuffled linear congruential
generator can generate in about a minute an a half. And note a
shuffled linear congruential generator does better than DES.
---------------------------------------------------------------------------
10.0 HOW GOOD A RANDOM GENERATOR DO YOU NEED?
The exact determination of whether a particular random number generator
will be adequate for your needs is highly dependent upon your
application. There is a trade off between resolution and quantity. At
high resolution, it doesn't take very many random numbers before some
non-randomness can be detected. At low resolution, it takes quite a
large number of random numbers before some non-randomness can be
detected.
The chi-square tests above give us an idea of where the tradeoff lies.
Each test used an expected value of 10 balls per bin, and each ball
required DIM random numbers per ball, where DIM is the dimension of the
chi-square test being performed. The number of bins per dimension
determined the resolution required of the random numbers.
For example, the MTH$RANDOM number generator generates floating point
deviates between 0 and 1. The floating point number returned is an
F_floating datum type with 24 bits of accuracy. If we are doing a 1-D
chi-square test with 128 bins, we would convert the floating point
random number into an integer between 0 and 127 and then add 1. This
represents a loss of 17 bits of resolution, since we essentially
truncate the low 17 bits of the F_floating datum and keep the upper 7
bits which represents a number between 0 and 127. Because of this loss
of resolution, we can generate a lot more random numbers before any
non-randomness can be detected.
The chi-square tests above on the linear congruential generators
without shuffling were able to detect non-randomness in the generators
after about one million values were generated. This value is somewhat
approximate, but it gives you an idea of the order of magnitude above
which you cannot expect the sequence to appear random. If you only
need a small number of random numbers, say less than 100,000, if you do
not demand high resolution of the numbers, and if you're not concerned
about correlations between the initial SEED value and the first few
numbers from the generator, or about occasional correlations between
successive runs of the generator, or about the numbers all lying in
parallel lines and parallel planes, then the simple MTH$RANDOM or ANSI
C random number generators may be adequate for your needs. Despite
it's numerous drawbacks, this simple method of random number generation
is quite adequate for many applications.
If you are concerned about the quality of your random numbers, or if
you need more than 100,000 random numbers, or if you want to make sure
there is no easy relation between the initial SEED value and the first
few numbers output from the generator, or if you want to make sure
there will be no relation between successive runs of the generator,
then you can shuffle the numbers using the technique in section 8
(Shuffling). This produces excellent results, the amount of memory
required is quite small, and the small extra amount of time required is
usually quite trivial.
By shuffling the numbers you may be able to get a million or more
satisfactory random numbers from the MTH$RANDOM or ANSI C random number
generators. If significantly more than a million random numbers are
needed, then a closer look at your application is warranted to
determine if the random number generator will be adequate for the
purpose. But don't despair. Just because 50 million shuffled random
numbers from the ANSI C generator would fail all of the chi-square
tests doesn't mean they won't be quite adequate for you application.
---------------------------------------------------------------------------
11.0 WHAT IS A RANDOM NUMBER GENERATOR?
The field of Information Theory is where we find a mathematical
understanding of what a random number generator really is. Our ideal
random number generator is a zero-memory information source with
maximum entropy. We think of a source as emitting a sequence of
symbols from a fixed finite set S = {s1, s2, ... sn}. Successive
symbols are selected according to some fixed probability law. For a
zero-memory source, the source remembers nothing about what previous
symbols have been emitted; each successive symbol emitted from the
source is completely independent of all previous output symbols.
Maximum entropy indicates that each symbol in the set S is
equally likely to occur. Any other probability distribution would
result in a lower entropy, with some symbols being more likely than
others. You may recall entropy is also used in thermodynamics to
measure the amount or "randomness" in a particular system.
A linear congruential random number generator is an example of a
first-order Markov process. The next value depends entirely on
the previous value. Even the Data Encryption Standard (DES) is nothing
more than a first-order Markov process. A first-order Markov source
can emulate a zero-memory source up to a point.
A higher order Markov process has a better chance of emulating a
zero-memory process. By shuffling the output, we have converted the
first-order Markov process to something like an Mth-order Markov
process where M is the size of the shuffling deck we use. In this
case M pieces of information about previous outputs are kept instead of
just 1 piece of information. Thus, because of the additional memory
used, the output can resemble more closely that of a zero-memory
source. And indeed it does do a much better job.
For further understanding of zero-memory information sources and Markov
information sources consult a good book on Information Theory. The
book by Abramson in the references below is an excellent one which also
gives an interesting account of how languages such as English, Spanish,
French, and German, can be synthesized using a low order Markov
process.
---------------------------------------------------------------------------
12.0 CONCLUSIONS
In this paper we have examined some probability theory, the concept of
a random sequence coming from a random process, and ways of analyzing a
sequence of numbers to determine how likely it is they came from a
random process.
We then examined the inside of a linear congruential random number
generator which is the most prevalent type of random number generator
in use today. We dissected 5 commonly found linear congruential
generators, and then analyzed them using the chi-square and KS tests.
We found it was possible to detect non-randomness in each of these
generators after a few million iterations, except for the defective
RANDU generator, which failed the 3-D test after less than 800
iterations.
A technique for shuffling the random numbers was given and tested with
excellent results. Shuffling the output of any random number generator
using the technique given is therefore highly recommended.
Analysis of the Data Encryption Standard (DES) as a random number
generator was also performed. A software implementation of the Data
Encryption Standard DES was found to be too slow for any practical use
and the randomness of the results was found to be less than that
obtained by shuffling the output of a good linear congruential
generator.
Finally, a brief mathematical explanation of what an ideal random
number generator really is and how it differs from a computerized
random number generator was given. For further reading on this topic,
consult a book on Information theory and read about the different types
of information sources defined there.
===========================================================================
REFERENCES:
NUMERICAL RECIPES: The Art of Scientific Computing
William H. Press,
Harvard-Smithsonian Center for Astrophysics
Brian P. Flannery,
EXXON Research and Engineering Company
Saul A. Teukolsky,
Department of Physics, Cornell University
William T. Vetterling,
Polaroid Corporation
Cambridge University Press, 1986
This book contains a wealth of knowledge along with source code
examples on how to solve various kinds of programming problems.
There are several programming language versions of this book
available, including FORTRAN, Pascal, C, and BASIC. Machine
readable forms in various languages of all the sample routines
are also available.
THE ART OF COMPUTER PROGRAMMING: Vol 2. Seminumerical Algorithms
Knuth, Donald E.
Stanford University.
2nd edition. Reading, Mass.: Addison-Wesley 1981.
Knuth's 3 volume set on computer programming is the comprehensive
source which all others turn to.
INFORMATION THEORY AND CODING
Abramson, Norman
Associate Professor of Electrical Engineering, Stanford University
McGraw-Hill Book Company, Inc. 1963 McGraw-Hill electronic
Sciences Series
An excellent book on Information Theory. Describes a zero memory
source and an Nth order Markov process. Also shows how language
can be synthesized by an Nth order Markov process. For example,
you could synthesize something which sounded French, and if you
didn't know French, you'd say it was French. The same is true
for random numbers; you can synthesize them.
THE C PROGRAMMING LANGUAGE
Kernighan, Brian W.
Ritchie, Dennis M.
Englewood Cliffs, N.J. : Prentice Hall, c1978.
2nd edition. Englewood Cliffs, N.J. : Prentice Hall, c1988.
Considered the book which defines the C programming language.
"A Revised Algorithm for the Spectral Test" Algorithm AS 193.
Hopkins, T.R.
University of Kent at Canterbury, UK.
'Applied Statistics' Vol 32, pg 328-335 1983
This article gives a version of the spectral test coded in FORTRAN.
The program modified for VAX/VMS is presented in SPECTRAL.FOR
THE GENERATION OF RANDOM VARIATES
Newman, Thomas G.
associate professor of mathematics Texas Tech University
Odell, Patrick L.
professor of mathematics and statistics Texas Tech University
New York, Hafnew Publishing Company 1971
Series Title: Griffin's statistical monographs & courses ; no. 29.
DISTRIBUTION SAMPLING FOR COMPUTER SIMULATION
Lewis, Theodore Gyle
University of Southwestern Louisiana
Lexington, Massachusetts, D.C. Heath and Company 1975
ADVANCED ENGINEERING MATHEMATICS
Erwin Kreyszig,
Professor of Mathematics, Ohio State University, Columbus Ohio.
5th edition. John Wiley & Sons, 1983.
Chapter 23: Probability and Statistics
PROBABILITY AND STOCHASTIC PROCESSES FOR ENGINEERS
Carl W. Helstrom
Department of Electrical Engineering and Computer Sciences,
University of California, San Diego.
New York, Macmillan Publishing Company, 1984
Advanced probability theory. College or graduate level.
PROBABILITY AND STATISTICS
Edwards, Allen L.
Professor of Psychology, University of Washington
New York. Holt, Rinehart and Winston, Inc. 1971
A first course in probability and statistics, applied statistics
and experimental design.
DETECTION OF SIGNALS IN NOISE
Whalen, Anthony D.
Bell Telephone Laboratories
New York, Academic Press, Inc. 1971
A graduate level book. Chapter 1 is a review of probability and
chapter 2 describes random processes.
NOISE AND ITS EFFECT ON COMMUNICATION
Blachman, Nelson M.
Sylvania Electronic Systems, Mountain View, California
McGraw-Hill Book Company 1966
McGraw-Hill Electronic Sciences Series
ENCYCLOPEDIA OF COMPUTER SCIENCE
edited by Anthony Ralston
Petrocelli (New York, 1976)
See:
Random Number Generation (pp. 1192-1197)
by G. Marsaglia