Pseudorandom Number Generators and the Metropolis Algorithm

1
Pseudorandom Number Generators and the
Metropolis Algorithm

• Jane Ren

2
Outline

• Introduction to random number generators
• The Marsaglia random number generator
• SPRNG (A Scalable Library for Pseudorandom Number
  Generation)
• Textbook assignment 3.2.3 with Marsaglia random
  numbers
• Assignment 3.2.3 with random numbers from another
  type of random number generator.
• Test results
• Conclusions

3
Properties of Good Random Number Generators

• Randomness
• Uniformity
• Computational Efficiency
• Long period length
• Unpredictability
• Repeatability
• Portability
• Homogeneity
• Good random numbers are not easy to find

4
The Marsaglia Random Number Generator

• The programs in our textbook use the Marsaglia
  generator.
• A combinations of two generators
• Lagged Fibonacci series
• In In-r In-s mod 224, r 97, s 33
• Arithmetic Series
• I k, I 2k, I 3k, mod (224 - 3), k
  7654321
• The final step
• Take xn from a lagged Fibonacci series
• Take yn from a n arithmetic series
• If x gt y, the final number is x y
• Else the final number is x y 1

5
Sample Output from the Marsaglia Generator

• Assignment a0102_02 with seed1 1 and seed2
  6.
• Run 1
• RANMAR INITIALIZED.
• idat 1, xr 0.894103587
• idat 2, xr 0.034489155
• idat 3, xr 0.330096781
• idat 4, xr 0.613509536
• extra xr 0.710489452
• Run 2
• MARSAGLIA CONTINUATION.
• idat 1, xr 0.710489452
• idat 2, xr 0.429685593
• idat 3, xr 0.142630935
• idat 4, xr 0.369295061
• extra xr 0.894588530

6
Some Advantages of he Marsaglia Generator

• The lagged Fibonacci series has a period of (224
  1) 294 2110. This is very good for a
  pseudorandom number generator.
• The Marsaglia random number generator is portable
  because it uses the IEEE standard representation
  of 32 bit floating point numbers.

7
Some Drawbacks of the Marsaglia Generator

• Low discretization step size 2-24
• We have already seen this problem from assignment
  a0204_03.
• Not as fast as some other generators, such as the
  congruential generator.
• It fails the birthday spacing test.
• The birthday spacing test chooses m birthdays in
  a year of n days. If a certain day appeared more
  than once, then that value should be
  asymptotically Poisson distributed with mean
  m3/(4n).

8
3.2.3 Assignment for section 3.2 Page 152

• a0302_01
• Write a Metropolis program to generate the
  normal distribution through the Markov process
• x ? x x 2a(xr 0.5)
• where xr is a uniformly distributed random
  number in the range 0,1) and a is a real
  constant. Use a 3.0 and the initial value x
  0.0 to test your program. (i) Generate a chain
  of 2000 events. Plot the empirical peaked
  distribution function for these event in
  comparison with the exact peaked distribution
  function. Monitor also the acceptance rate. (ii)
  Repeat the above for a chain of 20000 events.

9
Assignment a0302_01

• Output for a chain of 20 000 events
• RANMAR INITIALIZED.
• Metropolis Generation of Gaussian random numbers
• Ndat20000, a3.0000
• Acceptance rate found 0.491700
• CPU time for running the Metropolis Program with
  Marsaglia 0.126u (average of 10 runs)
• Comparison between the empirical peaked
  distribution and the exact peaked normal
  distribution

10
Scalable Library for Pseudorandom Number
Generation (SPRNG)

• Software package released under the GNU General
  Public License.
• There are 6 kinds of random number generators in
  SPRNG
• Combined Multiple Recursive Generator
• 48 Bit Linear Congruential Generator with Prime
  Addend
• 64 Bit Linear Congruential Generator with Prime
  Addend
• Modified Lagged Fibonacci Generator
• Multiplicative Lagged Fibonacci Generator
• Prime Modulus Linear Congruential Generator
• Mersenne Twister Generator (will be added)
• Link http//sprng.cs.fsu.edu/

11
Why SPRNG?

• The random number generators in SPRNG are
  considered the top generators.
• The generators have never failed the standard
  tests.
• Including some of the largest random number tests
  with up to 1013
• The tests used were collisions, coupon,
  equidistribution, gap, maximum of t,
  permutations, poker, runs up, serial, sum of
  distributions, Ising Metropolis, Ising Wolff,
  and Random Walk.
• What if one of the SPRNG generators was used for
  assignment a0302_01?
• All of the random number generators passed the
  empirical tests, so which one do we want to use?

12
Descriptions of Some of the Tests SPRNG
Generators passed

• Permutation Test
• Divide the random sequence into subsequences of a
  certain length.
• Rank each value in the subsequence by magnitude.
• Count the frequency each permutation occurs.
• Conduct a Chi-Square test with these counts.
• Compare the results.
• Collision Test
• Measure the number of times the number falls at a
  position that already had a number in it.
• Coupon Collectors Test
• Compare actual period with expected period

13
The Generators

• Combined Multiple Recursive Generator
• Period 2219
• Fast because it is a linear congruential
  generator
• 48 Bit Linear Congruential Generator with Prime
  Addend
• Period 248
• Number of distinct streams is of the order of 219
  
• Fast because it is a linear congruential
  generator
• 64 Bit Linear Congruential Generator with Prime
  Addend
• Period 264
• Number of distinct streams is of the order of 108
  
• Fast because it is a linear congruential
  generator

14
The Generators

• 4. Modified Lagged Fibonacci Generator
• Period 231(2l 1), where l is the lag.
• The default period is 21310
• Number of distinct streams is of the order of
  231(l 1)-1 , default 21310
• 5. Multiplicative Lagged Fibonacci Generator
• Period 261(2l 1), default 281
• Number of distinct streams is of the order of
  263(l 1)-1 , default 281
• 6. Prime Modulus Linear Congruential Generator
• Period 264 - 2
• Number of distinct streams 258
• Fast because it is a linear congruential
  generator

15
Is there a better generator?

• The current 6 random number generators in SPRNG
  work sufficiently well.
• The Mersenne Twister generator MT19937 is not a
  part of SPRNG yet, but it is being added to
  SPRNG.
• I chose to use Mersenne Twister as a substitute
  for the Marsaglia generator for a0302_01.

16
Why Mersenne Twister?

• Astronomical period of 219937 1
• No other generator has achieved this!
• Marsaglia has 2110
• 623 dimensional equidistribution with up to 32
  bit accuracy in a working area of only 624 words.
• No other generator has achieved this!
• Can be employed in practical important
  simulations because it is based on a good
  definition of randomness.
• Many other generators are only random in a
  particular area.
• Performance exceeds that of integer-large-modulus
  generators.
• Efficient random number generation.
• Passed tests many stringent tests, including the
  diehard tests (invented by Marsaglia) and tests
  on the higher dimensional uniformity, including
  the spectral test and the k-distribution test.
• Most suitable for Monte Carlo because it
  emphasizes on the most significant bits.

17
Mersenne Twister The Algorithm

• Based on the recurrence (mod 2)
• l, s, t, and u are integers, b and c are bit
  masks of word size, x0n-1 is an array of n
  unsigned integers of word size, and u, ll, and a
  are unsigned constant integers of word size.

18
Tempering Matrix

• The Mersenne Twister uses a tempering matrix, so
  that multiplication by the tempering matrix is
  very efficient and fast.
• Example of a tempering Matrix
• Let x (xw-1, xw-2, xw-3,, x0)
• a (aw-1, aw-2, aw-3,, a0)
• xA can be calculated with only binary operations.
• xA shiftright(x) if x0 0
• xA shiftright(x) XOR a if x0 1

19
Mersenne Twister How does It Work?

• Step 0
• Create a mask for the upper bits and another one
  for the lower bits.
• Step 1
• The x array is initialized with a seed
• Step 2
• y is the upper bits of xi concatenated with the
  lower bits of xi1.

20
Mersenne Twister How does It Work?

• Step 3
• y A.
• A is chosen carefully so it can be easily
  calculated with right shift and exor.
• Step 4
• Multiply with T (tempering matrix) for better
  equidistribution and good acurracy.
• Steps 5 6
• Increment i by 1 and
  repeat the process

21
Mersenne Twister

• We use an improved version of the Mersenne
  Twister from its original authors (released
  2/2002).
• Some disadvantages of the original version
• Most significant bits only affect most
  significant bits of the array state.
• Non-overlapping subsequences on successive runs
  is not possible.
• Not as fast as the improved version.

22
Sample Output from the Mersenne Twister

• Run 1
• MERSENNE TWISTER INITIALIZED.
• idat 1, xr 0.814723692
• idat 2, xr 0.135477004
• idat 3, xr 0.905791934
• idat 4, xr 0.835008590
• extra xr 0.126986812
• Run 2
• MERSENNE TWISTER CONTINUATION.
• idat 0, xr 0.126986812
• idat 1, xr 0.968867771
• idat 2, xr 0.913375856
• idat 3, xr 0.221034043
• extra xr 0.632359250

23
Mersenne Twister in a302_01

• Output for a chain of 20 000 events
• RANMAR INITIALIZED.
• Mersenne Twister Generation of Gaussian random
  numbers
• Ndat20000, a3.0000
• Acceptance rate found 0.497850
• The acceptance rate is 0.4917000 with the
  Marsaglia generator
• CPU Time for the Metropolis Program with Mersenne
  Twister 0.124 s (Average of 10 runs)
• Marsaglia 0.126 s

24
Mersenne Twister in a302_01

• Comparison between the empirical peaked
  distribution and the exact peaked normal
  distribution

25
Results from Marsaglia and Mersenne Twister
Combined
26
Mersenne Twister and Marsaglia Execution Times

• The CPU times for running the Metropolis program
  with Marsaglia and Mersenne Twister are very
  close.
• Does it mean both generators have roughly the
  same speed?
• Mersenne Twister is faster than
  Marsaglia.
• I wrote a program to only generate 4,294,967,295
  Marsaglia numbers and another one that only
  generates 4,294,967,295 Mersenne Twister numbers.
  
• Mersenne Twister 244.710u
• Marsaglia 312.580u
• MT is approximately 22 faster.

27
Gaussian Difference Test

• Data from the Metropolis program
• Marsaglia
• Mean 0.003543
• Error bar 0.007169
• Mersenne Twister
• Mean -0.008677
• Error bar 0.006998
• Q 0.222553

28
GNUPLOT FOR DISTRIBUTION FUNCTION.

• STMC_DF_GNU(Ndat, data1)
• STMC_DF_GNU(Ndat2, data2)

29
The Two-Sided Kolmogorov Test

• N N1 N2 50000
• KOLM2_AS2(N1, N2, DAT1, DAT2, DEL, Q)
• DAT1 Marsaglia, DAT2 Mersenne Twister
• DEL 0.001172, Q 0.882084
• DAT1 Marsaglia1, DAT2 Marsaglia2
• DEL 0.001148, Q 0.896565
• DAT1 MT1, DAT2 MT2
• DEL 0.002050, Q 0.243914
• DEL is the maximum deviation from the exact
  distribution function
• Q is the approximation to the probability that
  this deviation is due to chance.

30
Conclusions

• The Mersenne Twister generator is more suitable
  for Monte Carlo Simulations than the Marsaglia
  random number generator.
• Mersenne Twister is faster.
• Mersenne Twister more accurate.
• Mersenne Twister has a linear discretization.
• Mersenne Twister is probably the best generator
  at the present time.

31
References

• Berg, Bernd. Markov Chain Monte Carlo
  Simulations and Their Statistical Analysis.
  World Scientific, 2004.
• Karaoglu, Hihat. Mersenne Twister A Study on
  Random Number Generators. Project Prestudy.
  Vrije Universiteit Brussel, Apr 2004.
• Korab, Holly. Best of the Random Number
  Generators. The National Center for
  Supercomputing Applications. Featured story, Feb
  1999.
• LEcuyer Pierre. Uniform Random Number
  Generation. The Annals of Operations Research
  (1994).
• Lee, J.C. Random Number and Statistical Tests.
  Talk at Chuo University, Japan. Jan, 2003.
• Matsumoto, Makoto and Takuji Nishimura. Mersenne
  Twister A 623-Dimensionally Equidistributed
  Uniform Pseudo-Random Number Generator. ACM
  Transactions on Modeling and Computer
  Simulations, Vol. 8, No. 1, Jan 1998, Pages 3-30.

32
References (continued)

• 7. Matsumoto, Makoto and Takuji Nishmura.
  rng_afm.c. An Improved version of the Mersenne
  Twister. http//www.math.keio.ac.jp/matumoto/emt.h
  tml
• 8. Mascagni, Michael, and Srinivasan, Ashok.
  Algorithm 806 SPRNG A Scalable Library for
  Pseudorandom Number Generation. ACM Transactions
  on Mathematical Software, Vol 26, No. 3,
  September 2000, pages 436-461.
• RSA Laboratories Bulletin News and Advice on
  Data Security and Cryptography. Number 8. Sept
  3, 1998.
• Sinclair, Bart. Permutation Test. Class notes
  for Elec 428 Computer Systems Performance.
  http//www.owlnet.rice.edu/elec428/rng/perexp.htm
  l
• The Scalable Parallel Random Number Generator
  Library (SPRNG) for ASCI Monte Carlo
  Computations. http//sprng.cs.fsu.edu/

33
Questions