Ejemplo natural de malos resultados con un generador de números aleatorios de Lehmer

Question

Le digo a mis alumnos todos los años que hay cierta correlación entre empates sucesivos en un RNG de Lehmer, por lo que debe utilizar un Mersenne Twister o MWC256 de Marsaglia... pero soy incapaz de dar un ejemplo natural que no Lehmer. ¿Por supuesto existen pruebas especialmente diseñadas que Lehmer generadores fallan, pero alguien pueden dar una situación natural donde la autocorrelación da lugar a resultados aberrantes?

Gracias por tus pensamientos.

Answer 1

Finalmente me ocurrió una pregunta que enviar Wichman-Hill generador fuera de la carretera. No es tan natural como se puede desear, pero espero que sea suficientemente espectacular.

Aquí está el problema: el estudio de la distribución de $$ X = -10U_1 - 22U_2 + 38U_3 - 3U_4 + U_5 + 4U_6 - 38 U_7$$ con la $U_i$ iid uniforme en $(0,1)$.

Vamos a dibujar un histograma.

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" )

Ahora intenta con 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" )

Sí, todos los valores generados son entero:

> 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

Si algunas personas muestran interés, me puede explicar brevemente cómo he construido este ejemplo.

---

Aquí es un boceto de la construcción del ejemplo. Lineal Congruential Generadores de confiar en una secuencia en la $[1,\dots m-1]$ definido por $x_{n+1} = \alpha x_n \ [m]$ (lo que significa que el modulo $m$ )$\gcd(\alpha,m)=1$. Los números pseudo-aleatorios $u_n = {1\over m} x_n$ se comportan más o menos como los números aleatorios uniformes en $(0,1)$.

Un problema conocido de estos generadores es que usted puede encontrar un montón de coeficientes de $a_0, a_1, \dots, a_k\in \mathbb{Z}$ tal que $$ a_0 + a_1 \alpha + \cdots + a_k \alpha^k = 0 \ [m].$$ Esto se traduce en $a_0 x_n + \cdots + a_k x_{n+k} = 0 \ [m]$ todos los $n$, que se convierte en resultados en $a_0 u_n + \cdots + a_k u_{n+k} \in \mathbb{Z}$. Esto significa que todos los $(k+1)$-tuplas $(u_n, \dots, u_{n+k})$ de caída en los planos ortogonales a $(a_0, \dots, a_k)'$.

Por supuesto, con $k=1$, $a_0=\alpha$ y $a_1 = -1$, usted tiene un ejemplo. Pero como $\alpha$ es generalmente grande, esto no le dará un aspecto agradable en el ejemplo anterior. El punto es encontrar el "pequeño" $a_0, \dots, a_k$ valores.

Vamos $$L = \{ (a_0, \dots, a_k) \ :\ a_0 + a_1 \alpha + \cdots + a_k \alpha^k = 0 \ [m] \}.$$ Esta es una $\mathbb{Z}$-red.

El uso de un poco de álgebra modular, uno puede comprobar que $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)'$ es una base de este entramado (esta es la clave del resultado aquí...).

El problema es, entonces, encontrar una breve vector en $L$. He utilizado LLL algoritmo para este propósito. Algoritmo 2.3 de Brian Ripley Simulación Estocástica en punta (después de mi respuesta) por kjetil b halvorsen podría haber sido utilizado también.

Para Wichman-Hill, el teorema del resto Chino permite comprobar fácilmente que es equivalente a un generador de la anterior clase, con $\alpha = 16555425264690$$m = 30269\times30307\times30323 = 27817185604309$.

Answer 2

La aplicación de la mayoría de que se trate, con pequeñas desviaciones de correlación es la criptografía. La capacidad para predecir un valor pseudo-aleatorio con lo mejor de ordinario la precisión que se puede traducir en las habilidades superiores a romper los esquemas de cifrado.

Esto puede ser ilustrado con un tanto artificial ejemplo, diseñado para ayudar a la intuición. Se muestra cómo un verdadero cambio dramático en la previsibilidad se puede incurrir por arbitrariamente un pequeño correlación serial. Deje $X_1, X_2, \ldots, X_n$ ser iid estándar de variables Normales. Deje $\bar X = (X_1+X_2+\cdots+X_n)/n$ ser su media y designar

$$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}$). Este trade-off entre la hipótesis y el rendimiento es instructivo!