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

Aquí hay profundidades ocultas: {[18]]}

1. El uso del pequeño u en RAND_MAX + 1u. RAND_MAX se define como un tipo int, y a menudo es el mayor int posible. El comportamiento de RAND_MAX + 1sería undefined en tales casos como usted estaría desbordando un tipo signed. Escribir 1u fuerza la conversión de tipo de RAND_MAX a unsigned, evitando así el desbordamiento.

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

En conclusión, en estos días, debido a su manejabilidad, usaría % 6. No es probable que introduzca ninguna anomalía estadística 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::rand es un caso de balderdash de culto de carga heredada que debería hacer que tus cejas se levanten cada vez que lo veas.

Hay varias cuestiones aquí:

La gente de contrato generalmente asume-incluso las pobres almas desafortunadas que no saben nada mejor y no pensarán en ello precisamente en estos términos-es que rand muestras de la distribución uniforme en los enteros en 0, 1, 2, ..., RAND_MAX, y cada llamada produce un 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, las implementaciones históricamente fallaron en proporcionar incluso el simulacro más simple de independencia. Por ejemplo, C99 §7.20.2.1 'La función rand' dice, sin elaboración:

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

Esta es una oración sin sentido, porque la pseudoaleacion es una propiedad de una función (o familia de funciones), no de un entero, pero eso no impide que incluso los burócratas ISO abusen del lenguaje. Después de todo, los únicos lectores que estarían molestos por ello saben mejor que leer la documentación de rand por 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 a pesar de que 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) produce 1 con 100% de probabilidad, que no es el caso de muestras aleatorias independientes en cualquier distribución soportada en enteros pares e impares. Por lo tanto, un culto de carga surgió que uno debe descartar la bits de orden bajo para perseguir a la bestia escurridiza de 'mejor aleatoriedad'. (Alerta de Spoiler: Este no es un término técnico. Esta es una señal de que cualquiera que sea la prosa que estás leyendo o no sabe de lo que están hablando, o piensa que no tiene ni idea y debe ser condescendiente.)

El segundo problema es que incluso si cada llamada muestra independientemente de una distribución aleatoria uniforme en 0, 1, 2, ..., RAND_MAX, el resultado de rand() % 6 no se distribuiría uniformemente en 0, 1, 2, 3, 4, 5 como un rollo de dados, a menos que RAND_MAX sea congruente con -1 módulo 6. Contraejemplo simple: Si RAND_MAX = 6, entonces desde rand(), todos los resultados tienen igual 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 el muestreo de rechazo: repetidamente dibujar una muestra aleatoria uniforme independiente s de 0, 1, 2, ..., RAND_MAX, y rechazar (por ejemplo) el resultados 0, 1, 2, ..., ((RAND_MAX + 1) % 6) - 1-si obtiene uno de esos, comience de nuevo; de lo contrario, ceda 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 de rand() que aceptamos es uniformemente divisible por 6, y cada posible resultado de s % 6 se obtiene por el mismo número de aceptados resultados de rand(), por lo que si rand() se distribuye uniformemente, entonces también lo es s. No hay ninguna obligado en el número de ensayos, pero el número esperado es menor que 2, y la probabilidad de el éxito crece exponencialmente con el número de ensayos.

La elección de que resultados de rand() rechazas es irrelevante, siempre que mapees un número igual de ellos a cada entero por debajo de 6. El código en cppreference.com hace una elección diferente , debido al primer problema anterior - que nada está garantizado sobre la distribución o independencia de las salidas de rand(), y en la práctica los bits de orden bajo 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 en rollos de troquel si rand() produce una distribución uniforme en 0, 1, 2, ..., RAND_MAX.

Ejercicio para el lector: ¿Por qué preferirías rechazar uno u otro subconjunto? ¿Qué cómputo se necesita para cada juicio en los dos casos?

Un tercer problema es que el espacio de la 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 posteriores, lo que hace que no parezcan tan aleatorios después de todo. Así que ni siquiera pienses en usar esto para criptografía.

Usted puede ir a la ruta de fantasía overengineered y clase de C++11 std::uniform_int_distribution con un dispositivo aleatorio apropiado y su motor aleatorio favorito como el siempre popular Mersenne twister std::mt19937 para jugar a los dados con su primo de cuatro años de edad, pero incluso eso no va a ser apto para generar material de clave criptográfica-y el Mersenne twister es un terrible cerdo espacial también con un estado de varios kilobytes causando estragos en la caché de su CPU con un tiempo de configuración obsceno, por lo que es malo incluso para, por ejemplo, , simulaciones paralelas de Monte Carlo con árboles reproducibles de subcomputaciones; su popularidad probablemente surge principalmente de su nombre pegadizo. Pero puedes usarlo para dados de juguete rodando así ejemplo!

Otro enfoque es utilizar un simple generador de números pseudoaleatorios criptográficos 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 la investigación en ciencias naturales) de 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 sesgado que 1+std::rand()%6 en realidad es cierto. Así que escribí un programa de prueba para tabular los resultados de ambos métodos (no he escrito C++ en años, por favor compruébelo). Un enlace para ejecutar el código se encuentra aquí. También se reproduce como sigue:

// 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é la salida de esto y usé la función chisq.test en R para ejecutar una prueba de Chi-cuadrado para ver si los resultados son significativamente diferentes de lo esperado. Esta pregunta de stackexchange entra en más detalle del uso de la prueba de chi-cuadrado para probar la equidad de la matriz: ¿Cómo puedo probar si una matriz es justa?. Aquí están los resultados de algunas carreras:

> ?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 corridas que hice, el valor de p para ambos métodos siempre fue mayor que los valores alfa típicos utilizados para probar la significación (0.05). Esto significa que no consideraríamos que ninguno de ellos esté sesgado. Curiosamente, el método supuestamente imparcial tiene valores de p consistentemente más bajos, lo que indica que en realidad podría ser más sesgado. La advertencia es que solo hice 3 carreras.

ACTUALIZACIÓN: Mientras escribía mi respuesta, Konrad Rudolph publicó una respuesta que toma el mismo enfoque, pero obtiene un resultado muy diferente. No tengo la reputación de comentar su respuesta, así que voy a abordarla aquí. En primer lugar, lo principal es que el código que utiliza utiliza la misma semilla para el generador de números aleatorios cada vez que se ejecuta. Si cambia el semilla, en realidad se obtiene una variedad de resultados. Segundo, si no cambias la semilla, sino el número de ensayos, también obtienes una variedad de resultados. Trate de aumentar o disminuir en un orden de magnitud para ver lo que quiero decir. En tercer lugar, hay un cierto truncamiento entero o redondeo pasando donde los valores esperados no son muy precisos. Probablemente no sea suficiente para marcar la diferencia, pero está ahí.

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

Uno puede pensar en un generador de números aleatorios como trabajando en un flujo de dígitos binarios. El generador convierte la corriente en números cortándola en trozos. Si la función std:rand está trabajando con un RAND_MAX de 32767, entonces está usando 15 bits en cada segmento.

Cuando uno toma los módulos de un número entre 0 y 32767 inclusive uno encuentra que 5462 '0 'y '1 'pero solo 5461 '2', '3', ' 4 ' y '5'. Cuanto mayor sea el valor RAND_MAX, menor será el sesgo lo será, 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 entero en el rango 0..7 y rechazando 6 y 7.

Uno espera que cada bit en el flujo de bits tenga la misma probabilidad de ser un '0' o un '1' independientemente de dónde esté en el flujo o los valores de otros bits. Esto es excepcionalmente difícil en la práctica. Las diferentes las implementaciones de software PRNGs ofrecen diferentes compromisos entre velocidad y calidad. Un generador congruente lineal como std::rand ofrece la velocidad más rápida para la calidad más baja. Un generador criptográfico ofrece la más alta calidad para la velocidad más baja.