Generating normal random numbers, difference between Mersenne Twister and other RNGs

I have a question about the difference in generating normal (pseudo) random numbers with the mt19937ar algorithm and the other RNGs, e.g. the mrg32k3a algorithm.

For some reason the mt19937ar algorithm uses 1.0218 uniform (pseudo) random numbers on average while the mrg32k3a algorithm uses 2.0282.

Both algorithms use the Ziggurat algorithm to tranform the uniform random numbers into normal random numbers.

Can anyone explain the difference?

Test code:

stream0 = RandStream('mt19937ar','Seed',0); RandStream.setDefaultStream(stream0); a=rand(3e6,1);

stream0 = RandStream('mt19937ar','Seed',0); RandStream.setDefaultStream(stream0); b=randn(1e6,1); c=rand; find(a==c)

ans =

1021808

stream1 = RandStream('mrg32k3a','Seed',0); RandStream.setDefaultStream(stream1); a=rand(3e6,1);

stream1 = RandStream('mrg32k3a','Seed',0); RandStream.setDefaultStream(stream1); b=randn(1e6,1); c=rand; find(a==c)

ans =

2028211

[Daniel Shub]

I have no idea what is happening, but that is cool.

Answer 1

[Peter Perkins]

Jan, this requires some details that almost noone needs to think about, but you asked.

The Ziggurat algorithm, most of the time, requires a random integer to determine the level, and a random uniform to determine the position within that level. The original Ziggurat algorithm as published by Marsaglia used and reused 32bits of randomness for both purposes, and that’s what SHR3CONG still does (largely to be backwards compatible). Very fast, but if you look very closely, you can see a “griddiness” that’s a consequence of reusing bits. Almost certainly not a practical problem, though. The MT algorithm as implemented in MATLAB’s RandStream uses 64bits of randomness for those two values, and doesn’t reuse any bits. But the original author of MRG32K3A didn’t do that, L’Ecuyer recommended using two full U(0,1)’s, so that’s what MATLAB does for MRG32K3A, and it uses up 128 bits of randomness.

Of course, something like 1.5% of the time, the RandStream Ziggurat algorithm requires more randomness, so that’s why you don’t see exactly 1 (MT) or 2(MRG32K3A) d.p. uniform per normal.

If you switch to ‘Inversion’, you’ll find that exactly one uniform gets used up for each normal in all cases, but it’s slower to do that computation.