
*NE582 Monte Carlo* *Return to Course Outline*

Lesson 3: Pseudo-Random Numbers

This lesson is the one and only lesson in this course that has to do with the INPUT to the "black box" introduced in Lesson 2:

which is a sequences of uniformly distributed random numbers between 0 and 1 that we denote with the Greek letter xsi (x).

As implied by the figure, these form the material from which every Monte Carlo method creates its answer. And, as with any human enterprise, the quality of the inputs has a strong influence on the quality of the outputs. In the early days of Monte Carlo (the 40s and 50s), it was common for researchers to suggest physical devices--based on truly random physical phenomena such as nuclide decay--to generate truly random numbers. As computers have gotten faster and faster, this has gotten less and less feasible. (For nuclide decay for example, we need to DETECT about 50 particles for each random number generated. If we were to try to keep up with a 1 gigaHz clock speed at, say, a 10% counting efficiency, we would need a (50 particles detected/random number)*(1.e9 numbers/sec needed)/(10% counting efficiency) ~= 10 curie source!

The most common sources of "random" numbers nowadays are the pseudo-random numbers (which can be designed to be "random enough").

Pseudo-random number generation

Linear congruential generators

For the most part, the reading in the Numerical Recipes book (which I recommend you get a copy of if you plan to do much computational work in your career) covers the material well. With the homework problem P-1, I think you should have covered all of this subject that we want to.

Here are the principal ideas that I want you to remember:

Any random number generator from a computer is not truly random (since it is repeatable).
This repeatability is actually very important to us for validation reasons. Also, this repeatability is the basis of some Monte Carlo perturbation methods.
The most common form is that of "linear congruential generators" which generate a series of bounded integers, with each one using the one before it with the formula:

which makes 1 impossible, but 0 possible.

There are several problems that can show up in random number generators:

Periodicity: The string of numbers is finite in length and so will repeat itself eventually.

Bad choices of a,b, and m can lead to correlation among successive numbers (e.g., a very low number always followed by another very low number). It can be shown that a "full" period of m can be assured by choosing:

a(mod 4)=1

m as a power of 2

b odd.

The finiteness of the string also means that certain domains of the problem may be unreachable.

We have to be careful not to depend on the later digits of a uniform deviate: they may be less random than the early digits.

Example: Find the pseudo-random series corresponding to m=8, a=5, b=3, i0=1.
Answer:

# i ai+b Mod m

0 1 8 0 0.000

1 0 3 3 0.375

2 3 18 2 0.25

3 2 13 5 0.625

4 5 28 4 0.5

5 4 23 7 0.875

6 7 38 6 0.75

7 6 33 1 0.125

Extending the period using more than one generator

One consequence of use of the linear congruential generators is that no more than m random numbers can be generated before the series begins to repeat; a second consequence is that m has to less than or equal to the largest integer that the computer can handle (since you have to do arithmetic on the computer using it). Together, these tell us that we can have periods of no more than the largest integer that the computer can handle, which is

, where N is the number of bits used to denote integers.

Fortunately, you can get around this limitation by combining the results of several generators. If you use N different generators, add the N random numbers, and keep the fractional part, the resulting random number is properly distributed and the series has a period equal to the PRODUCT of the periods of the N generators. For example, the FORTRAN subroutine used to get a (double precision) random number for the KENO and MORSE computer codes at ORNL is the following:

      DOUBLE PRECISION FUNCTION FLTRN()
C#######################################################################
C
C                AUDIT TRAIL INFORMATION
C
C     DATE THE MODULE WAS LAST PERMANENTLY UPDATED:     95/05/04
C     TIME THE MODULE WAS LAST PERMANENTLY UPDATED:     16:20:35
C     PROGRAMMER NAME:                                  L.M.PETRIE
C     MODULE NAME:                                      ULFLTRN
C     CURRENT ARCHIVING LEVEL NUMBER:                   00002
C     CURRENT NUMBER OF PERMANENT UPDATES:              00002
C     DATE OF LAST ACCESS BY LIBRARIAN:                 95/05/04
C     DATASET NAME: X4S.SCALE4.MASTER
C
C#######################################################################
      INTEGER X, Y, Z, X0, Y0, Z0, XM, YM, ZM
      DOUBLE PRECISION V, X1, Y1, Z1
      EXTERNAL RANDNUM
      COMMON /QRNDOM/ X, Y, Z
      SAVE /QRNDOM/
      PARAMETER ( MX = 157, MY = 146, MZ = 142 )
      PARAMETER ( X0 = 32363, Y0 = 31727, Z0 = 31657 )
      PARAMETER ( X1 = X0, Y1 = Y0, Z1 = Z0 )
      X      = MOD(MX*X,X0)
      Y      = MOD(MY*Y,Y0)
      Z      = MOD(MZ*Z,Z0)
      V      = X/X1 + Y/Y1 + Z/Z1
      FLTRN = V - INT(V)
      RETURN
      END

Close examination of this coding reveals that it uses 3 linear congruential generators (with b=0, m=prime number, and period=m). The individual generators have periods of about 32000 (which makes them usable for 16-bit integer computer systems), but combined they have a limiting period of over 3xE13, which should be large enough for any (present-day) application.

Quasi-random number generation

In addition to the pseudo-random sequences, another family of "random numbers" that can be used (with appropriate care) to feed our Monte Carlo process are the quasi-random numbers. I left them for the end of the lecture because:

1. They are not really random, so they have to be used very carefully.
2. They are not really used very much for Monte Carlo methods in radiation transport.

But, because they are of such importance to the general field of Monte Carlo (i.e., beyond transport methods), you should have a taste of what they are, how to get them, and what they are useful for.

"Quasi"-random numbers are not random at all; they are a deterministic, equal-division, sequence that are "ordered" in such a way that they can be used by Monte Carlo methods.

The easiest to generate are Halton sequences, which are found by "reflecting" the digits of prime base counting integers about their radix point. Clear as mud, right? A simple example makes it perfectly clear how to do it.

Okay. You pick a prime number to be the base (e.g., 2). You then simply count in that base, like the first two columns below:

Base 10	Base 2	"Reflected" Base 2	x=Reversed Base 10
1	1	0.1	0.5
2	10	0.01	0.25
3	11	0.11	0.75
4	100	0.001	0.125
5	101	0.101	0.625
6	110	0.011	0.375
7	111	0.111	0.875
8	1000	0.0001	0.0625
...	...	...	...

The third column shows what we mean by "reflecting" the number. When translated back to base 10, it becomes the next random number in the sequence.

If you follow the formula, you can see that after M=2^N-1 numbers in the sequence, the sequence consists of the equally spaced fractions (i/(M+1), i=1,2,3,...M). Therefore, the sequence is actually deterministic, not stochastic. The beauty of the sequence is that the re-ordering results in the an even coverage of the domain (0,1) even if the number of random numbers chosen does not correspond to 2^N-1.

Because you need ALL of the numbers in the sequence to cover the (0,1) range (which is not true of the pseudo-random sequence), it is important that all of the numbers of the sequence be used on the SAME decision. For a Monte Carlo process that has more than one decision, a different sequence must be used for each decision. The most common way this situation is handled is to use each of the low prime numbers in order--i.e., use the sequence based on 2 for the first decision, the one based on 3 for the second, then 5, 7, 11, etc. for the subsequent decisions.

For reference, here is a subroutine that delivers the Halton sequence for base 2:

      function rand2(r)                                                 rand2 1
      data idig/0/                                                      rand2 2
      real digit(50)                                                    rand2 3
      integer idigit(50)                                                rand2 4
      double precision rn,prev                                          rand2 5
      if(r.ne.0.)then                                                   rand2 6
        idig=0                                                          rand2 7
      endif                                                             rand2 8
      n=2                                                               rand2 9
      rn=1.d0/2.d0                                                      rand2 10
      prev=1.d0                                                         rand2 11
      if(idig.eq.0)then                                                 rand2 12
        do i=1,50                                                       rand2 13
          idigit(i)=0                                                   rand2 14
          prev=prev*rn                                                  rand2 15
          digit(i)=prev                                                 rand2 16
        enddo                                                           rand2 17
        idig=1                                                          rand2 18
      endif                                                             rand2 19
      idigit(1)=idigit(1)+1                                             rand2 20
      i=1                                                               rand2 21
      dowhile(idigit(i).eq.n)                                           rand2 22
        idigit(i)=0                                                     rand2 23
        i=i+1                                                           rand2 24
        idigit(i)=idigit(i)+1                                           rand2 25
        if(i.gt.idig)idig=i                                             rand2 26
      enddo                                                             rand2 27
      rand2=0.                                                          rand2 28
      do i=1,idig                                                       rand2 29
        rand2=rand2+digit(i)*idigit(i)                                  rand2 30
      enddo                                                             rand2 31
      return                                                            rand2 32
      end                                                               rand2 33

You can adapt it to other bases by changing the name and the "2" to the new base in lines 9 and 10. Normally, one would reset the sequence by passing any argument OTHER than 0 to the function; thereafter, passing 0 gets the next number in the sequence.

Pseudo or Quasi?

Well? Which should we use: pseudo-random or quasi-random? Most Monte Carlo particle transport methods that you will encounter have already made the decision for you, and most of them use pseudo-random. The reason for this is that particle transport problems are, to the Monte Carlo theorist, INFINITE dimensional problems. Before we can understand what this means and why this gives pseudo-random numbers the edge for our transport problems, we need to define a few terms.

Dimension, Hypercubes, and Discrepancy

Here are the definitions and a brief discussion about each of these important terms AS USED in Monte Carlo.

Dimension: The number of random numbers used for each estimate obtained.

Since this is a different meaning of the word "dimension" than we usually see, let's examine it a little closer.

We began the course by saying that a Monte Carlo method consists of two steps. The first step was to design a numerical experiment and the second was to run it. No matter how revolutionary, innovative, etc., we are in the first step, the second looks the same: We "draw out" a certain number of random numbers and, from them, produce a single number as the resulting estimate of whatever we are looking for.

The process is mathematically indistinguishable from estimating the average value of a function f(x1,x2,x3,...,xN) of N dimensions, each dimension having a domain (0,1). This is best seen in a example.

Example: Find the expected value of the sum of two fair dice (i.e., 6 sides with values 1, 2, 3, 4, 5, and 6 -- each of which is equally likely to come up).
Probability 101 point of view: The expected value of a single die is the average of its 6 values, or 3.5. Therefore, two dice have an expected value of 7.
Analog Monte Carlo point of view: Devise a Monte Carlo method (numerical experiment) that mimics the action:
1. Pick an integer value among the first 6 integers. (This could be accomplished by "pulling" a random number, multiplying it by 6, adding 1 to it, and discarding any fractional part.)
2. Pick another integer value among the first 6 integers. (Repeat step #1 with a different random number.)
3. Sum the values from step 1 and step 2 to form the first estimate, s1.
Repeat steps 1-3 N times. The resulting ANSWER is the sum of the scores divided by N.
Functional point of view: Looking at the Monte Carlo process we just ran, we can see that the process was one that turned two random numbers into an estimate. Although our THOUGHT processes involved chance and probability, the result does not--the same two random numbers would produce the same estimate every time. In other words, the Monte Carlo process actually defines a DETERMINISTIC function of the random numbers that are used. In this case, the function is a two dimensional histogram (which actually looks like the UT main library, for those of you who have seen it!):

The two dimensions, x and y, correspond to the two random numbers that are used in the Monte Carlo process and have domains limited to (0,1). The answer itself is the average of the function, which corresponds to the integral of the function:

The example is applicable to every Monte Carlo process. All them turn a set of random numbers into estimates, so all of them are equivalent to multi-dimensional functions. And the ANSWERS to the Monte Carlo processes are averages of the estimates, so correspond to integrations of the functions. Of course, the number of dimensions corresponds to the number of random numbers used to get the answer.

We can actually say a little more than this, because each of the dimensions are integrated from 0 to 1, so the functions are defined on a hypercube, which is a multidimensional domain where each dimension is bounded by 0 and 1.

Our final new term to introduce is that of the discrepancy of a sequence of random points in a hypercube (i.e., the random numbers used to get all of the estimates that go into making our Monte Carlo answer). Discrepancy is just a quantitative measue of how evenly the random points COVER the hypercube.

[Other concepts covered in slides but not yet introduced into lecture: offsetting, Koksma-Hwalka theorem, error bounds on Halton.]


*Return to Course Outline* © 2005 by Ronald E. Pevey. All rights reserved.