Natural example of bad results with a Lehmer Random Number Generator

I tell my students every year that there’s some correlation between successive draws in a Lehmer RNG, so they should use a Mersenne Twister or Marsaglia’s MWC256... but I am unable to provide a natural example where Lehmer would fail. Of course there are specially designed tests that Lehmer generators fail, but can someone provide a natural situation where the autocorrelation results in aberrant results?

• Did you look into chapter 1 of Brian Ripley:"Stochastis Simulation"? Nov 17 '14 at 10:29
• No. Can you provide some link or a quick summary? Nov 17 '14 at 10:34
• amazon.com/Stochastic-Simulation-Brian-D-Ripley/dp/0470009608/… This is THE standard reference for whay you asked about (you can find a copy here: gen.lib.rus.ec) Chapter one has some very interesting plots showing what can happen! Nov 17 '14 at 10:42
• I had a look. This is a nice reference. However, I should have precised that I was aware of the papers of Marsaglia on the topic. I was looking for an example that send a ''real life'' LCG off the road, not a toy example — or an historical example like RANDU. The point is ''exhibit coefficients that leads to this bad behaviour''. I think I partially succeeded with my example on Wichman-Hill below... Nov 17 '14 at 10:47
• I guess you can use that book to make an alghorithm for showing the problems! The point is that those random numbers will always accumulate on some collection of planes. One can make sure there are many of those planes, so the distances between them are short, but one cannot get rid of them. For most applications that do not matter much. It it matters for your application, you better find some better generator. Nov 17 '14 at 10:55

I finally came up with a question that send Wichman-Hill’s generator off the road. It is not as natural as one may wish but I hope it’s spectacular enough.

Here’s the problem: study the distribution of $$X = -10U_1 - 22U_2 + 38U_3 - 3U_4 + U_5 + 4U_6 - 38 U_7$$ with the $U_i$ iid uniform on $(0,1)$.

We will just draw an histogram.

x <- c(-10, -22, 38, -3, 1, 4, -38)
RNGkind(kind="Mersenne")
U <- matrix( runif(7*1e6), nrow= 7 )
hist( colSums(x * U), breaks = seq(-70,50,by=0.25), col="black" )


Now try with Wichman-Hill:

RNGkind(kind="Wichman")
U <- matrix( runif(7*1e6), nrow= 7 )
hist( colSums(x * U), breaks = seq(-70,50,by=0.25), col="black" )


Yep, all generated values are integer:

> head(colSums(x*U), 20)
[1]   3  -9 -52 -21 -23 -14 -18   8 -23  12   5   4 -17 -16 -19 -44   5   4 -23
[20] -15


If some people show interest, I may explain briefly how I constructed this example.

Here is a sketch of the construction of the example. Linear Congruential Generators rely on a sequence in $[1,\dots m-1]$ defined by $x_{n+1} = \alpha x_n \ [m]$ (which means modulo $m$), with $\gcd(\alpha,m)=1$. The pseudo random numbers $u_n = {1\over m} x_n$ behave roughly as uniform random numbers in $(0,1)$.

A known issue of these generators is that you can find lots of coefficients $a_0, a_1, \dots, a_k\in \mathbb{Z}$ such that $$a_0 + a_1 \alpha + \cdots + a_k \alpha^k = 0 \ [m].$$ This results in $a_0 x_n + \cdots + a_k x_{n+k} = 0 \ [m]$ for all $n$, which in turns results in $a_0 u_n + \cdots + a_k u_{n+k} \in \mathbb{Z}$. This means that all $(k+1)$-tuples $(u_n, \dots, u_{n+k})$ fall in planes orthogonal to $(a_0, \dots, a_k)'$.

Of course with $k=1$, $a_0=\alpha$ and $a_1 = -1$, you have such an example. But as $\alpha$ is usually big, this won’t give a nice looking example as above. The point is to find "small" $a_0, \dots, a_k$ values.

Let $$L = \{ (a_0, \dots, a_k) \ :\ a_0 + a_1 \alpha + \cdots + a_k \alpha^k = 0 \ [m] \}.$$ This is a $\mathbb{Z}$-lattice.

Using a little modular algebra, one can check that $f_0 = (m,0,\dots,0)'$, $f_1 = (\alpha,-1,0,\dots,0)'$, $f_2 = (0,\alpha,-1,0,\dots,0)'$, $\dots$, $f_k = (0,\dots,0,\alpha,-1)'$ is a base of this lattice (this is the key result here...).

The problem is then to find a short vector in $L$. I used LLL algorithm for this purpose. Algorithm 2.3 from Brian Ripley Stochastic Simulation pointed (after my answer) by kjetil b halvorsen could have been used as well.

For Wichman-Hill, the Chinese remainder theorem allows to check easily that it is equivalent to a generator of the above kind, with $\alpha = 16555425264690$ and $m = 30269\times30307\times30323 = 27817185604309$.

• Wow! Yes, please explain. Nov 17 '14 at 10:03

The application most concerned with tiny deviations from correlation is cryptography. The ability to predict a pseudo-random value with better than ordinary accuracy can translate into superior abilities to break encryption schemes.

This can be illustrated with a somewhat artificial example designed to help the intuition. It shows how a truly dramatic change in predictability can be incurred by an arbitrarily tiny serial correlation. Let $X_1, X_2, \ldots, X_n$ be iid standard Normal variables. Let $\bar X = (X_1+X_2+\cdots+X_n)/n$ be their mean and designate

$$Y_i = X_i - \bar X$$

as their residuals. It is elementary to establish that (1) the $Y_i$ have identical Normal distributions but (2) the correlation among $Y_i$ and $Y_j$ (for $i\ne j$) is $-1/(n-1)$. For large $n$ this is negligible, right?

Consider the task of predicting $Y_n$ from $Y_1, Y_2, \ldots, Y_{n-1}$. Under the assumption of independence, the best possible predictor is their mean. In fact, though, the construction of the $Y_i$ guarantees that their sum is zero, whence

$$\hat{Y_n} = -\sum_{i=1}^{n-1}Y_i$$

is a perfect (error-free) predictor of $Y_n$. This shows how even a tiny bit of correlation can enormously improve predictability. In cryptanalysis the details will differ--one studies streams of bits, not Normal variates, and the serial correlations can be tiny indeed--but the potential for equally dramatic results nevertheless exists.

For those who like to play with things, the following simulation in R replicates this hypothetical situation and summarizes the mean and standard deviations of the prediction errors. It also summarizes the first-order serial correlation coefficients of the simulated $Y_i$. For reference, it does the same for the $X_i$.

predict.mean <- function(x) mean(x[-length(x)])
predict.correl <- function(x) -sum(x[-length(x)])
n <- 1e5
set.seed(17)
sim <- replicate(1e3, {
x <-rnorm(n)
y <- x - mean(x)
c(predict.mean(y) - y[n], predict.correl(y) - y[n],
predict.mean(x) - x[n], predict.correl(x) - x[n],
cor(y[-1], y[-n]))
})
rownames(sim) <- c("Mean method", "Exploit",
"Mean (reference)", "Exploit (reference)",
"Autocorrelation")
zapsmall(apply(sim, 1, mean))
zapsmall(apply(sim, 1, sd))


As written, this will take a good fraction of a minute to run: it performs $1000$ simulations of the case $n=10^5$. The timing is proportional to the product of these numbers, so reducing that by an order of magnitude will give satisfactory turnaround for interactive experiments. In this case, because such a large number of simulations were performed, the results will be pretty accurate, and they are

Mean method             Exploit    Mean (reference) Exploit (reference)     Autocorrelation
-0.064697            0.000000           -0.064697            5.502087           -0.000046
Mean method             Exploit    Mean (reference) Exploit (reference)     Autocorrelation
1.0235              0.0000              1.0235            321.8314              0.0031


On average, estimating $Y_n$ with the mean of the preceding values made an error of $-0.06$, but the standard deviation of those errors was close to $1$. The exploit offered by recognition of the correlation was absolutely perfect: it always got the prediction right. However, when applied to the truly independent values $(X_i)$, the exploit performed terribly, with a standard deviation of $321.8$ (essentially equal to $\sqrt{n}$). This trade-off between assumptions and performance is instructive!

• Well... Even Mersenne twister is not suited for cryptography... :/ Oct 9 '14 at 20:38
• "It is elementary to establish that (1) the $Y_i$ are iid Normal" Perhaps I'm reading this wrong, but should that say "identically distributed" rather than "iid"? Jan 1 '15 at 23:43
• @Silver Good catch! You are absolutely right; these variables are manifestly not independent. (I have been trying to avoid jargon like "iid," but somehow it slipped in.) I'll make a suitable change.
– whuber
Jan 2 '15 at 0:10