April 4, 2017---Class 23--- Pseudorandom Numbers Generators, Serial Correlation Tests
Activities:
Properties of Pseudo Random Number Generators
---------------------------------------------
We looked at some properties of Linear Congruential PRNGs.
x_{n+1} = (a x_n + c) mod m
a is the multiplier
c is the increment
m is the modulus.
If c=0, this is a multiplicative generator, otherwise, it is mixed.
We start the sequence with a seed x_0.
m determines the maximum possible length of the sequence. It is
necessary to pick a and c carefully to get the maximum length
sequence.
Here is a theorem: The length of the period is m if the following
three conditions all hold.
a) c is relatively prime to m
b) (a-1) is a multiple of p for every prime factor of m
c) if m is a multiple of 4, (a-1) is a multiple of 4.
We looked at the case, m=32, a=3, c=4. This does not obey the
theorem and there is no sequence longer than 8. There is even a
sequence of length one. So much for randomness! To explore this
generator use the Mathematica definition
myrand[x_] := Mod[ 3x+4, 32]
The NestList command can be used to study a sequence:
seed = 3
NestList[myrand, seed, 32]
Adjusting a and c, we can easily find sequences of length 32;
however, that is not enough to ensure that the sequence looks random.
Another example to try is the following:
myrand[x_] := Mod[ 13x+9, 32]
A Mathematica notebook April1_2013.nb may be found in my home
directory with some of the tests we did, including how to make a
set of random numbers with a triangular probability density.
Serial Correlation Tests
------------------------
In addition to having a uniform probability density, a pseudo random
number generator should have zero autocorrelation. There are tests
based on looking at pairs or triplets of numbers. For linear
congruential random number generators, there is a well known
problem that like the rain in Spain they tend to fall mainly in
the planes. This can be demonstrated by showing a two dimensional
scatter plot of (r_n, r_{n+1}). If you would like to see this, there
is an executable in ~sg/bin called random_lc:
random_lc | axis m 0 | xplot
will show you the correlations. The code is in ~sg/src and this
makes use of a poor linear congruential random number generator. You
would not see the same pattern with the fsr random number generator.
We looked in detail at another test of correlations based on
triplets of numbers. If we divide the interval (0,1) into L bins,
and look at triplets in the resulting L^3 cells, we derived that
the probability of having a cell remain unoccupied is exponentially
falling. Specifially, if we examine the occupancy after tL^3
triplets have been generated, the number of unoccupied cells should
be L^3 Exp[-t].
The box_rand.c code in ~sg/src implements this code for L=10.
(You can change the define statement near the top of the code to
use other values of L.) To compile the code use the command:
cc -o box_rand box_rand.c
Running box_rand
----------------
We ran a few tests with the fsr random number generator and saw that
the number of unoccupied cells was decreasing. Then students wrote
a shell script to ran the code 20 or 50 times, putting the reults
in a file box_rand.apr3. Using awk and aver, it was easy to get
the average of the 50 trials and plot the result with a semilog
scale using axis. The result clearly shows exponential decay of
the number of unoccupied boxes.
Some notes about shell scripting
--------------------------------
Here is an example of a bash script to run box_rand 50 times
#!/bin/bash
for i in {1..50}; do
box_rand 10
done
The key new point is that the {1..50} notation saves us from typing
a long list of integers. In csh or tcsh we can use a while loop:
#!/bin/csh
set count = 1
while ($count <= 50)
echo $count
@ count++
end
This just uses echo to print the numbers at the terminal. Key issues
are initialization of count with the set command before the while
loop, and the way to increment count at the end of the loop in the
line before the end statement. @ at the beginning of a line is use
to tell the shell that some aritmetic calculation is coming. The
result does not need a dollar sign, but if you use an equal sign, you
need the dollar sign for shell variables on the RHS. ++ is used to
increment a variable, just like in C or C++. I always use a tab
after the @ symbol, but I think only a space is required. If you
don't want to use the convienience of ++, you can also do this:
@ count = $count + 1
The following was not actually done in class, however, you can find the
result for the fsr random number generator in ~sg/src/misc/fsr.autocorrs.ax
Please plot the results in that file.
Autocorrelations of Pseudo Random Number Generators, or
comparing drand48 and fsr
-------------------------------------------------------
Before class I prepared codes to test the drand48 random number
generator that comes with Unix. The codes can be found in my src/misc
directory. random_drand48.c puts out two numbers per line and
can be used for a scatter plot or to calculate autocorrelations.
box_drand48.c can be used to test triplets of random numbers as
we did before with box_rand.c for the fsr random number generator.
We compared the autocorrelation of the two random number generators
out to a lag of 20. Doing about 50 trials in each case, we found
that drand48 seemed to have a smaller dispersion than fsr, but that
it seemed to have more negative values of autocorrelation than
positive values. (I also passed around a similar graph from a test
of a linear congruential random number generator that had some much
larger values of the autorcorrelations for some lags.)
We created lists of random numbers as long as 250,000 and looked
at autocorrelations with lags as large as 250. We also histogrammed
the autocorrelations and found a distribution that appears Gaussian.
Here is a list of files used in an earlier class on this topic.
-rw-r--r-- 1 sg staff 920 Nov 13 2014 temp.nov13
-rw-r--r-- 1 sg staff 57 Nov 13 2014 sum_errors.a
-rw-r--r-- 1 sg staff 25750 Nov 13 2014 fsr.autocorrs.ax
-rw-r--r-- 1 sg staff 482962 Nov 13 2014 fsr.autocorrs
-rw-r--r-- 1 sg staff 25746 Nov 13 2014 drand48.autocorrs.ax
-rwxr-xr-x 1 sg staff 13856 Nov 13 2014 aver
-rw-r--r-- 1 sg staff 9607 Nov 13 2014 aver.c
-rw-r--r-- 1 sg staff 482944 Nov 13 2014 drand48.autocorrs
-rwxr-xr-x 1 sg staff 8640 Nov 13 2014 random_drand48
-rw-r--r-- 1 sg staff 1550 Nov 13 2014 box_rand16.ax
-rw-r--r-- 1 sg staff 20425 Nov 13 2014 box_drand16.dat
-rwxr-xr-x 1 sg staff 9008 Nov 13 2014 box_drand48
-rw------- 1 sg staff 1191 Nov 13 2014 box_drand48.c
-rw-r--r-- 1 sg staff 25532 Nov 13 2014 box_rand16.dat
-rwxr-xr-x 1 sg staff 13112 Nov 13 2014 box_rand16
-rwxr-xr-x 1 sg staff 9016 Nov 13 2014 box_rand