¿Por qué rand ()% 6 está sesgado?

Al leer cómo usar std :: rand, encontré este código en cppreference.com

int x = 7;
while(x > 6) 
    x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

¿Qué hay de malo en la expresión de la derecha? Lo probé y funciona perfectamente.

Hay dos problemas con rand() % 6( 1+no afecta a ninguno de los dos).

Primero, como se ha señalado en varias respuestas, si los bits bajos de rand()no son adecuadamente uniformes, el resultado del operador restante tampoco es uniforme.

En segundo lugar, si el número de valores distintos producidos por rand()no es un múltiplo de 6, el resto producirá más valores bajos que valores altos. Eso es cierto incluso si rand()devuelve valores perfectamente distribuidos.

Como ejemplo extremo, imagine que rand()produce valores distribuidos uniformemente en el rango [0..6]. Si observa los restos de esos valores, cuando rand()devuelve un valor en el rango [0..5], el resto produce resultados distribuidos uniformemente en el rango [0..5]. Cuando rand()devuelve 6, rand() % 6devuelve 0, como si rand()hubiera devuelto 0. De modo que obtiene una distribución con el doble de ceros que cualquier otro valor.

El segundo es el verdadero problema con rand() % 6.

La forma de evitar ese problema es descartar los valores que producirían duplicados no uniformes. Calcula el múltiplo más grande de 6 que es menor o igual a RAND_MAX, y siempre que rand()devuelve un valor que es mayor o igual a ese múltiplo, lo rechaza y llama a `rand () nuevamente, tantas veces como sea necesario.

Entonces:

int max = 6 * ((RAND_MAX + 1u) / 6)
int value = rand();
while (value >= max)
    value = rand();

Esa es una implementación diferente del código en cuestión, destinada a mostrar más claramente lo que está sucediendo.

Aquí hay profundidades ocultas:

1. El uso de lo pequeño uen RAND_MAX + 1u. RAND_MAXse define como un inttipo y suele ser el más grande posible int. El comportamiento de RAND_MAX + 1sería indefinido en casos en los que desbordaría un signedtipo. La escritura 1ufuerza la conversión de tipos de RAND_MAXa unsigned, evitando así el desbordamiento.

2. El uso de % 6 can (pero en todas las implementaciones std::randque he visto no lo hace ) introduce un sesgo estadístico adicional más allá de la alternativa presentada. Tales casos en los que % 6es peligroso son casos en los que el generador de números tiene llanuras de correlación en los bits de orden inferior, como una implementación de IBM bastante famosa (en C) de rand, creo, en la década de 1970, que cambió los bits altos y bajos como "un final florecer". Una consideración adicional es que 6 es muy pequeño cf. RAND_MAX, por lo que habrá un efecto mínimo si RAND_MAXno es un múltiplo de 6, que probablemente no lo sea.

En conclusión, en estos días, por su manejabilidad, usaría % 6. No es probable que introduzca anomalías estadísticas más allá de las introducidas por el propio generador. Si aún tiene dudas, pruebe su generador para ver si tiene las propiedades estadísticas adecuadas para su caso de uso.

Este código de ejemplo ilustra que std::randes un caso de legado culto al cargo que debería hacer que sus cejas se eleven cada vez que lo vea.

Hay varios problemas aqui:

El contrato que la gente suele asumir, incluso las pobres almas desventuradas que no saben nada mejor y no pensarán en ello precisamente en estos términos, es que las randmuestras de la distribución uniforme de los números enteros en 0, 1, 2,… RAND_MAX,, y cada llamada produce una muestra independiente .

El primer problema es que el contrato asumido, muestras aleatorias uniformes e independientes en cada llamada, no es realmente lo que dice la documentación y, en la práctica, históricamente las implementaciones no lograron proporcionar ni el más mínimo simulacro de independencia. Por ejemplo, C99 §7.20.2.1 'La randfunción' dice, sin más detalles:

La randfunción calcula una secuencia de enteros pseudoaleatorios en el rango de 0 a RAND_MAX.

Esta es una oración sin sentido, porque la pseudoaleatoriedad es una propiedad de una función (o familia de funciones ), no de un número entero, pero eso no impide que incluso los burócratas de ISO abusen del lenguaje. Después de todo, los únicos lectores a los que les molestaría saber que no deben leer la documentación randpor temor a que sus células cerebrales se descompongan.

Una implementación histórica típica en C funciona así:

static unsigned int seed = 1;

static void
srand(unsigned int s)
{
    seed = s;
}

static unsigned int
rand(void)
{
    seed = (seed*1103515245 + 12345) % ((unsigned long)RAND_MAX + 1);
    return (int)seed;
}

Esto tiene la desafortunada propiedad de que , aunque una sola muestra puede distribuirse uniformemente bajo una semilla aleatoria uniforme (que depende del valor específico de RAND_MAX), alterna entre enteros pares e impares en llamadas consecutivas, después de

int a = rand();
int b = rand();

la expresión (a & 1) ^ (b & 1)da 1 con 100% de probabilidad, lo que no es el caso de muestras aleatorias independientes en cualquier distribución compatible con enteros pares e impares. Por lo tanto, surgió un culto de carga en el que uno debería descartar los bits de bajo orden para perseguir a la escurridiza bestia de la "mejor aleatoriedad". (Alerta de spoiler: este no es un término técnico. Es una señal de que la prosa que estás leyendo no sabe de lo que están hablando o piensa que no tienes ni idea y debes ser condescendiente).

El segundo problema es que incluso si cada llamada muestreó independientemente de una distribución aleatoria uniforme en 0, 1, 2,… RAND_MAX, el resultado derand() % 6 no se distribuiría uniformemente en 0, 1, 2, 3, 4, 5 como un dado. tirar, a menos que RAND_MAXsea ​​congruente con -1 módulo 6. Contraejemplo simple: si RAND_MAX= 6, entonces desde rand(), todos los resultados tienen la misma probabilidad 1/7, pero desde rand() % 6, el resultado 0 tiene probabilidad 2/7 mientras que todos los demás resultados tienen probabilidad 1/7 .

La forma correcta de hacer esto es con muestreo de rechazo: extraer repetidamente una muestra aleatoria uniforme e independiente sde 0, 1, 2,… RAND_MAX, y rechace (por ejemplo) los resultados 0, 1, 2,…, ((RAND_MAX + 1) % 6) - 1—si obtiene uno de esos, empezar de nuevo; de lo contrario, cede s % 6.

unsigned int s;
while ((s = rand()) < ((unsigned long)RAND_MAX + 1) % 6)
    continue;
return s % 6;

De esta manera, el conjunto de resultados rand()que aceptamos es divisible por 6, y cada resultado posible de s % 6se obtiene por el mismo número de resultados aceptados de rand(), por lo que si rand()se distribuye uniformemente, entonces también lo ess . No hay límite en el número de ensayos, pero el número esperado es menor que 2 y la probabilidad de éxito aumenta exponencialmente con el número de ensayos.

La elección de cuál los resultados de la rand()rechaza es irrelevante, siempre y cuando se asigna el mismo número de ellos a cada número entero inferior a 6. El código en cppreference.com hace una diferente opción, debido al primer problema por encima de la que no se garantiza nada acerca de la distribución o independencia de salidas de rand(), y en la práctica, los bits de bajo orden exhibieron patrones que no "parecen lo suficientemente aleatorios" (no importa que la siguiente salida sea una función determinista de la anterior).

Ejercicio para el lector: Demostrar que el código en cppreference.com produce una distribución uniforme sobre rodillos de matriz si rand()los rendimientos de una distribución uniforme en 0, 1, 2, ..., RAND_MAX.

Ejercicio para el lector: ¿Por qué preferiría rechazar uno u otro subconjunto? ¿Qué cálculo se necesita para cada ensayo en los dos casos?

Un tercer problema es que el espacio semilla es tan pequeño que incluso si la semilla se distribuye uniformemente, un adversario armado con el conocimiento de su programa y un resultado, pero no la semilla, puede predecir fácilmente la semilla y los resultados subsiguientes, lo que hace que no parezcan tan al azar después de todo. Así que ni siquiera pienses en usar esto para criptografía.

Puede seguir la elegante ruta de ingeniería excesiva y la std::uniform_int_distributionclase de C ++ 11 con un dispositivo aleatorio apropiado y su motor aleatorio favorito, como el siempre popular tornado Mersenne, std::mt19937para jugar a los dados con su primo de cuatro años, pero incluso eso no va a funcionar. estar en forma para generar material de claves criptográficas, y el tornado de Mersenne también es un terrible acaparador de espacio con un estado de varios kilobytes que causa estragos en la memoria caché de su CPU con un tiempo de configuración obsceno, por lo que es malo incluso para, por ejemplo , simulaciones de Monte Carlo en paralelo con árboles reproducibles de subcomputaciones; su popularidad probablemente se deba principalmente a su pegadizo nombre. ¡Pero puedes usarlo para tirar dados de juguete como este ejemplo!

Otro enfoque es usar un generador de números pseudoaleatorios criptográfico simple con un estado pequeño, como un simple borrado rápido de clave PRNG , o simplemente un cifrado de flujo como AES-CTR o ChaCha20 si está seguro ( por ejemplo , en una simulación de Monte Carlo para investigación en ciencias naturales) que no hay consecuencias adversas para predecir resultados pasados ​​si el estado alguna vez se ve comprometido.

No soy un usuario experimentado de C ++ de ninguna manera, pero estaba interesado en ver si las otras respuestas con respecto a std::rand()/((RAND_MAX + 1u)/6)ser menos sesgadas de lo que 1+std::rand()%6realmente son ciertas. Así que escribí un programa de prueba para tabular los resultados de ambos métodos (no he escrito C ++ en años, verifíquelo). Aquí se encuentra un enlace para ejecutar el código . También se reproduce de la siguiente manera:

// Example program
#include <cstdlib>
#include <iostream>
#include <ctime>
#include <string>

int main()
{
    std::srand(std::time(nullptr)); // use current time as seed for random generator

    // Roll the die 6000000 times using the supposedly unbiased method and keep track of the results

    int results[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

        results[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results[n] << ' ';
    }

    std::cout << "\n";


    // Roll the die 6000000 times using the supposedly biased method and keep track of the results

    int results_bias[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()%6;

        results_bias[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results_bias[n] << ' ';
    }
}

Luego tomé el resultado de esto y usé la chisq.testfunción en R para ejecutar una prueba de Chi-cuadrado para ver si los resultados son significativamente diferentes de lo esperado. Esta pregunta de intercambio de pila explica con más detalle el uso de la prueba de chi-cuadrado para probar la equidad del dado: ¿Cómo puedo probar si un dado es justo? . Estos son los resultados de algunas ejecuciones:

> ?chisq.test
> unbias <- c(100150, 99658, 100319, 99342, 100418, 100113)
> bias <- c(100049, 100040, 100091, 99966, 100188, 99666 )

> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 8.6168, df = 5, p-value = 0.1254

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 1.6034, df = 5, p-value = 0.9008

> unbias <- c(998630, 1001188, 998932, 1001048, 1000968, 999234 )
> bias <- c(1000071, 1000910, 999078, 1000080, 998786, 1001075   )
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.051, df = 5, p-value = 0.2169

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 4.319, df = 5, p-value = 0.5045

> unbias <- c(998630, 999010, 1000736, 999142, 1000631, 1001851)
> bias <- c(999803, 998651, 1000639, 1000735, 1000064,1000108)
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.9592, df = 5, p-value = 0.1585

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 2.8229, df = 5, p-value = 0.7273

En las tres carreras que hice, el valor p para ambos métodos fue siempre mayor que los valores alfa típicos utilizados para probar la significancia (0.05). Esto significa que no consideraríamos a ninguno de ellos como parcial. Curiosamente, el método supuestamente imparcial tiene valores p consistentemente más bajos, lo que indica que en realidad podría estar más sesgado. La advertencia es que solo hice 3 carreras.

ACTUALIZACIÓN: Mientras escribía mi respuesta, Konrad Rudolph publicó una respuesta que adopta el mismo enfoque, pero obtiene un resultado muy diferente. No tengo la reputación de comentar su respuesta, así que lo abordaré aquí. Primero, lo principal es que el código que usa usa la misma semilla para el generador de números aleatorios cada vez que se ejecuta. Si cambia la semilla, en realidad obtiene una variedad de resultados. En segundo lugar, si no cambia la semilla, pero cambia el número de ensayos, también obtendrá una variedad de resultados. Intente aumentar o disminuir en un orden de magnitud para ver a qué me refiero. En tercer lugar, hay un truncamiento o redondeo de enteros en los que los valores esperados no son del todo precisos. Probablemente no sea suficiente para marcar la diferencia, pero está ahí.

Básicamente, en resumen, simplemente obtuvo la semilla correcta y el número de pruebas que podría estar obteniendo un resultado falso.

Se puede pensar en un generador de números aleatorios como si trabajara en una secuencia de dígitos binarios. El generador convierte la transmisión en números dividiéndola en trozos. Si la std:randfunción está trabajando con unRAND_MAX 32767, entonces está usando 15 bits en cada segmento.

Cuando uno toma los módulos de un número entre 0 y 32767 inclusive, se encuentra que 5462 '0's y' 1's pero solo 5461 '2's,' 3's, '4's y' 5's. Por tanto, el resultado está sesgado. Cuanto mayor sea el valor RAND_MAX, menos sesgo habrá, pero es ineludible.

Lo que no está sesgado es un número en el rango [0 .. (2 ^ n) -1]. Puede generar un número (teóricamente) mejor en el rango 0..5 extrayendo 3 bits, convirtiéndolos en un número entero en el rango 0..7 y rechazando 6 y 7.

Se espera que cada bit del flujo de bits tenga la misma probabilidad de ser un '0' o un '1' independientemente de dónde se encuentre en el flujo o de los valores de otros bits. Esto es excepcionalmente difícil en la práctica. Las diferentes implementaciones de software PRNG ofrecen diferentes compromisos entre velocidad y calidad. Un generador congruencial lineal como std::randofrece la velocidad más rápida con la calidad más baja. Un generador criptográfico ofrece la máxima calidad a la menor velocidad.