From: dmurdoch@mast.queensu.ca (Duncan Murdoch)
Subject: Re: Nonrandom Random Numbers?
Organization: Queen's University
Date: Tue, 12 Jan 1993 21:13:55 GMT

In article <1993Jan12.140114.620@hnrc.tufts.edu> jerry@hnrc.tufts.edu (Jerry Dallal) writes:
>For those who do not have the Times, the article says, " . . . the scientists
>showed that five of the most popular computer programs for generating streams
>of random numbers produced errors when they were used in a simple mathematical
>model of the behavior of atoms in a magnetic crystal."  

>Unfortunately (and understandably) the article does not identify the 5 
>algorithms.  However, the article cites the original report from Ferrenberg et
>al. in the Dec. 7 issue of Physical Review Letters.  I would be grateful if
>someone with access to PRL would post the five generators.

The paper is Ferrenberg, A. M., Landau, D.P., and Wong, Y.J. (1992), Monte 
Carlo simulations:  hidden errors from "good" random number generators,
Physical Review Letters 69 no. 23, 3382-3384.

The generators tested were:

CONG:  A linear congruential generator
         X_n = (16807 X_{n-1}) mod (2^{31} - 1)

R250:  A shift register algorithm:
         X_n = X_{n-103} XOR X_{n-250}

R1279: Another shift register algorithm:
         X_n = X_{n-1063} XOR X_{n-1279}

In R250 and R1279, "XOR" is a bitwise exclusive OR; it doesn't say how 
many bits were used, but I'd assume 32 from the context.

SWC:   A subtract with carry generator (Marsaglia et al, Comput Phys Commun
60, p 345, 1990):

         X_n = X_{n-22} - X_{n-43} - C
         if X_n >= 0, C=0,
         if X_n < 0, X_n = X_n + (2^32-5), C=1

SWCW:   A combined subtract with carry Weyl generator (ibid):

         Z_n = X_{n-22} - Z_{n-43} - C
	 if Z_n >= 0, C=0
         if Z_n < 0, Z_n = Z_n+(2^32-5), C=1
	 Y_n = (Y_{n-1} - 362436069) mod 2^32
	 X_n = (Z_n - Y_n) mod 2^32

The simulations they did were of an Ising model on a 16 x 16 lattice with
periodic boundary conditions, where exact results are known.  They estimated 
the internal energy and specific heat using the Wolff cluster-flipping 
algorithm (Phys Rev Lett 62, p 361, 1989) the Swendsen-Wang algorithm 
(Phys Rev Lett 58, p 86, 1987) and a single-spin-flop Metropolis algorithm 
with sequential updating (see e.g. D.P. Landau in Applications of the Monte 
Carlo Method, K. Binder ed., Springer-Verlag 1987).

They used CONG to initialize the other generators.

Between 5 and 10 runs of 10^7 updates were performed.

They found the Wolff algorithm produced estimates which were 
systematically biased for internal energy and/or specific heat for all 
random number generators except CONG; the biases were all small (around 1 
part in 10^-4), but large in terms of standard errors (from 3 to 100 
standard errors from the true value for at least one of those measures).  
The CONG algorithm was within 1 standard error on both measures.

CONG did badly (more than 13 s.e.) when used in the Metropolis algorithm 
but okay in Swendsen-Wang; R250 and SWC were fine in Metropolis but bad in 
S-W  (6 and 4 s.e. respectively).  However, overall the problems weren't as 
severe as with the Wolff algorithm.

The authors conclude (with their italics marked _so_):

"In conclusion, extensive Monte Carlo simulations on an Ising model for 
which the exact answers are known have shown that ostensibly high quality 
random number generators may lead to subtle, but _dramatic_, systematic 
errors for some algorithms, but not others.  Since there is no reason to 
believe that the model which we have investigated has any special 
idiosyncracies, these results offer another stern warning about the need to 
very carefully test the implementation of new algorithms.  In particular, 
this means that a specific algorithm must be tested together with the random 
number generator being used _regardless_ of the tests which the generator 
has passed."

Duncan Murdoch
dmurdoch@mast.queensu.ca
Newsgroups: sci.math.stat