repetir - ¿Libc generador de números aleatorios defectuoso?
srand en c (3)
Considere un algoritmo para probar la probabilidad de que un cierto número se elija de un conjunto de N números únicos después de un número específico de intentos (por ejemplo, con N = 2, ¿cuál es la probabilidad en la ruleta (sin 0) de que se necesita X intentos para Negro para ganar?).
La distribución correcta para esto es pow (1-1 / N, X-1) * (1 / N).
Sin embargo, cuando pruebo esto usando el siguiente código, siempre hay una zanja profunda en X = 31, independientemente de N, e independientemente de la semilla.
¿Se trata de un error intrínseco que no se puede evitar debido a los aspectos específicos de la implementación del PRNG en uso, es un error real o estoy pasando por alto algo obvio?
// C
#include <sys/times.h>
#include <math.h>
#include <stdio.h>
int array[101];
void main(){
int nsamples=10000000;
double breakVal,diffVal;
int i,cnt;
// seed, but doesn''t change anything
struct tms time;
srandom(times(&time));
// sample
for(i=0;i<nsamples;i++){
cnt=1;
do{
if((random()%36)==0) // break if 0 is chosen
break;
cnt++;
}while(cnt<100);
array[cnt]++;
}
// show distribution
for(i=1;i<100;i++){
breakVal=array[i]/(double)nsamples; // normalize
diffVal=breakVal-pow(1-1/36.,i-1)*1/36.; // difference to expected value
printf("%d %.12g %.12g/n",i,breakVal,diffVal);
}
}
Probado en un Xubuntu 12.10 actualizado con el paquete libc6 2.15-0ubuntu20 e Intel Core i5-2500 SandyBridge, pero descubrí esto hace algunos años en una máquina Ubuntu más antigua.
También probé esto en Windows 7 usando Unity3D / Mono (aunque no estoy seguro de qué versión Mono), y aquí la zanja ocurre en X = 55 cuando usa System.Random, mientras que Unity.Random incorporado de Unity no tiene una zanja visible (al menos no para X <100).
La distribución:
Las diferencias:
Como han señalado otros, random() no es lo suficientemente aleatorio.
Usar los bits más altos en lugar de los más bajos no ayuda en este caso. Según el manual ( man 3 rand ), las implementaciones antiguas de rand() tenían un problema en los bits inferiores. Es por eso que se recomienda random() su lugar. Sin embargo, la implementación actual de rand() usa el mismo generador que random() .
Probé el uso correcto recomendado del viejo rand() :
if ((int)(rand()/(RAND_MAX+1.0)*36)==0)
... y consiguió la misma zanja profunda en X = 31
De manera interesante, si mezclo los números de rand() con otra secuencia, me deshago de la zanja:
unsigned x=0;
//...
x = (179*x + 79) % 997;
if(((rand()+x)%36)==0)
Estoy usando un viejo generador lineal congruente . Elegí al azar 79, 179 y 997 de una tabla de números primos. Esto debería generar una secuencia de repetición de longitud 997.
Dicho esto, este truco probablemente introdujo algo de no aleatoriedad, algo de huella ... La secuencia mixta resultante seguramente fallará en otras pruebas estadísticas. x nunca toma el mismo valor en iteraciones consecutivas. De hecho, se necesitan exactamente 997 iteraciones para repetir cada valor.
'''' [..] no se deben generar números aleatorios con un método elegido al azar. Debería usarse alguna teoría ". (DEKnuth," El arte de la programación de computadoras ", vol.2)
Para las simulaciones, si quiere estar seguro, use el Mersenne Twister
Esto se debe a que la función random() glibc no es lo suficientemente aleatoria. Según esta página , para los números aleatorios devueltos por random() , tenemos:
o i = (o i-3 + o i-31 ) % 2^31
o:
o i = (o i-3 + o i-31 + 1) % 2^31 .
Ahora tome x i = o i % 36 , y suponga que la primera ecuación de arriba es la que se usa (esto ocurre con un 50% de probabilidad para cada número). Ahora si x i-31 =0 y x i-3 !=0 , entonces la probabilidad de que x i =0 sea menor que 1/36. Esto se debe a que el 50% del tiempo o i-31 + o i-3 será menor que 2 ^ 31, y cuando eso suceda,
x i = o i % 36 = (o i-3 + o i-31 ) % 36 = o i-3 % 36 = x i-3 ,
que es distinto de cero. Esto hace que la zanja que ve 31 muestras después de una muestra 0.
Lo que se mide en este experimento es el intervalo entre las pruebas exitosas de un experimento de Bernoulli, donde el éxito se define como random() mod k == 0 para algunos k (36 en el OP). Desafortunadamente, se ve empañado por el hecho de que la implementación de random() significa que los ensayos de Bernoulli no son estadísticamente independientes.
Escribiremos rnd i para la salida i th de `random () ''y notamos que:
rnd i = rnd i-31 + rnd i-3 con probabilidad 0.75
rnd i = rnd i-31 + rnd i-3 + 1 con probabilidad 0.25
(Vea a continuación un resumen de la prueba).
Supongamos que rnd i-31 mod k == 0 y actualmente estamos viendo rnd i . Entonces debe ser el caso de que rnd i-3 mod k ≠ 0 , porque de lo contrario habríamos contado el ciclo como una longitud k-3 .
Pero (la mayoría de las veces) (mod k): rnd i = rnd i-31 + rnd i-3 = rnd i-3 ≠ 0 .
Por lo tanto, el ensayo actual no es estadísticamente independiente de los ensayos anteriores, y el 31º ensayo después de un éxito es mucho menos probable que tenga éxito que en una serie imparcial de ensayos de Bernoulli.
El consejo habitual en el uso de generadores lineales congruentes, que en realidad no se aplica al algoritmo random() , es usar los bits de orden superior en lugar de los bits de orden inferior, porque los bits de orden superior son "más aleatorios" ( es decir, menos correlacionado con valores sucesivos). Pero tampoco funcionará en este caso, porque las identidades anteriores son igualmente válidas para la función de high log k bits como para la función mod k == low log k bits .
De hecho, podríamos esperar que un generador lineal congruente funcione mejor, particularmente si usamos los bits de alto orden de la salida, porque aunque el LCG no es particularmente bueno en las simulaciones de Monte Carlo, no sufre la retroalimentación lineal de random() .
Algoritmo random , para el caso por defecto:
Que el state sea un vector de largos sin firmar. Inicialice el state 0 ...state 30 utilizando una semilla, algunos valores fijos y un algoritmo de mezcla. Para simplificar, podemos considerar que el vector de estado es infinito, aunque solo se usan los últimos 31 valores, por lo que en realidad se implementa como un búfer en anillo.
Para generar rnd i : (Note: ⊕ is addition mod 2 32 .)
state i = state i-31 ⊕ state i-3
rnd i = (state i - (state i mod 2)) / 2
Ahora, tenga en cuenta que:
(i + j) mod 2 = i mod 2 + j mod 2 si i mod 2 == 0 o j mod 2 == 0
(i + j) mod 2 = i mod 2 + j mod 2 - 2 si i mod 2 == 1 y j mod 2 == 1
Si i y j se distribuyen uniformemente, el primer caso ocurrirá el 75% del tiempo y el segundo el 25%.
Entonces, por sustitución en la fórmula de generación:
rnd i = (state i-31 ⊕ state i-3 - ((state i-31 + state i-3 ) mod 2)) / 2
= ((state i-31 - (state i-31 mod 2)) ⊕ (state i-3 - (state i-3 mod 2))) / 2 o
= ((state i-31 - (state i-31 mod 2)) ⊕ (state i-3 - (state i-3 mod 2)) + 2) / 2
Los dos casos pueden reducirse aún más a:
rnd i = rnd i-31 ⊕ rnd i-3
rnd i = rnd i-31 rnd i-3 + 1
Como en el caso anterior, el primer caso ocurre el 75% del tiempo, asumiendo que rnd i-31 y rnd i-3 se extraen independientemente de una distribución uniforme (que no lo son, pero es una primera aproximación razonable).