Quality of linear congruential generators for random numbers

Question

I'm doing some simulations of the Langevin equation, for various external forces. Being told that C's rand() from stdlib.h can introduce bias in my results, I'm using a Mersenne Twister.

Nevertheless, I would like to know (and see) exactly what kind of errors a linear congruential generator can introduce in my simulation. These are things that I have tried:

• Generating 3D tuples of randoms to try to see hyperplanes. I can't see anything.
• Doing the FFT of a large vector of random numbers. It's almost the same for both the Mersenne Twister and rand().
• Checking the equipartition principle for a particle in Brownian motion. Both integrators agree in the expected value of $\langle \text{KE}\rangle=\frac{1}{2}k_BT$ with the same number of significant digits.
• Seeing how well they bin in a number of bins that is not a power two. Both give the same qualitative results, no one being better.
• Looking at Brownian paths to see clear divergences from $\langle x\rangle = 0$. Again, no luck.
• Distribution of points in a circle. Filled, and only in the perimeter. Between all of them and between nearest neighbours (Shor's answer, below in the comments). Available in this gist, just run it with Julia 0.5.0 after installing the needed libraries (see the gist for instructions).

I would like to emphasize that I am looking for introduced bias in the context of physical simulations. For example, I have seen how rand() fails miserably the dieharder tests while the Mersenne Twister doesn't, but for the moment that doesn't mean too much for me.

Do you have any physical, concrete, examples on how a bad random number generator wrecks a Montecarlo simulation?

Note: I have seen how PRNG's like RANDU can be awful. I'm interested on not obvious examples, of generators that look innocent but ultimately introduce bias.

Answer 1

One interesting reference that describes a failure of a Monte Carlo simulation of a physical system due to inadequate RNG (although they didn't use an LCG) is:

A. Ferrenberg and D. P. Landau. Monte Carlo Simulations: Hidden Errors from "Good" Random Number Generators. Physical Review Letters 63(23):3382-3384, 1992.

The Ising models that Ferrenberg and Landua studied are good tests of RNG's because you can compare with an exact solution (for the 2-D problem) and find errors way out in the digits. These models should show the faults in an old fashioned 32 bit arithmetic PMMLCG without too much difficulty.

Another interesting reference is:

H. Bauke and Stephan Mertens. Pseudo Random Coins Show More Heads Than Tails. arXiv:cond-mat/0307138 [cond-mat.stat-mech]

Bauke and Mertens make a strong case against binary linear feedback shift register style random number generators. Bauke and Mertens have some other papers related to this.

It can be difficult to find the Marsaglia planes in a 3D scatter plot. You can try to rotate the plot to get a better view and sometimes they'll just pop out at you. You can also do 3D tests of statistical uniformity- for the older 32 bit LCG's, these will fail at quite small numbers of bins. e.g. a uniformity test with a 20x20x20 grid of bins in 3 dimensions is sufficient to detect lack of uniformity for the widely used (and previously well regarded) PMMLCG with m=2^31-1, a=7^5.

Answer 2

It is possible to use the TestU01 suite of PRNG tests in order to find out which of those tests rand fails. (See TestU01: A C Library for Empirical Testing of Random Number Generators for an overview of the test suite.) That's easier than coming up with Monte Carlo simulations of one's own. In a way, it's also a question of software composability (and software correctness): given a PRNG that appears to work okay on small, simple tests, how do you know its pathological behaviours won't be triggered by a larger program?

Here is the code:

#include "TestU01.h"

int main() {
  // Same as rand() on my machine
  unif01_Gen* gen = ulcg_CreateLCG(2147483647, 16807, 0, 12345);

  bbattery_SmallCrush(gen);
  bbattery_Crush(gen);

  return 0;
}

For the SmallCrush suite, there are 3 tests failing out of 15 (see guidelongtestu01.pdf in TestU01 for long descriptions and all the references; these are 15 statistics from 10 tests).

• BirthdaySpacings: bin $n$ $t$-dimensional vectors into $d^t$ bins, generating $d^t$ bin counts, $I_1,\ldots$, then the number of collisions among $\{I_{j+1}-I_j\}$ follows approximately a known distribution.

• Collision: generate $n$ $t$-dimensional vectors in $[0,1)^t$, bins them into $d^t$ equal hypercubes, and count the number of collisions.

• MaxOft: generate $n$ groups of $t$ values in $[0,1)$, compute the maximum $X$ in each group, and compare the distribution of the $n$ maxima with the theoretical distribution $\mathbb{P}(X<x) = x^t$; the values are $n=2\times 10^6$, $t=6$. The comparison is done via a chi-square test and an Anderson-Darling test: the $\chi^2$ test fails with p-value $<10^{-300}$, while the AD test passes (this AD test fails in some of the larger similar tests in Crush).

Assuming these are all "typical" Monte Carlo simulations (though they might not be like the problems you had in mind), the conclusion is that rand fails some unknown subset of them. I don't know why it is specifically that subset, though, so it's impossible for me to say whether it will work on your own problem or not.

MaxOft seems particularly suspicious, given how straightforward the description is.

Among the tests in the Crush suite, rand fails 51 out of 140 (140 statistics in 96 tests). Some of the failed tests (like Fourier3) are done on bit strings, so perhaps it's possible they wouldn't be relevant to you. Another curious test that fails is GCD, which tests the distribution of the GCD of two random integers. (Again, I don't know why this particular test fails or whether your simulation will suffer from this.)

P.S.: Yet another thing to note is that rand() is actually slower than some PRNGs that successfully pass all the SmallCrush, Crush, BigCrush tests, such as MRG32k3a (see the L'Ecuyer & Simard paper above).