c++ - randomize - random() c

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

Question

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 con la expresión de la derecha? Lo probé y funciona perfectamente.

Answer 1

Aquí hay profundidades ocultas:

1. El uso de la pequeña u en RAND_MAX + 1u . RAND_MAX se define como un tipo int y, a menudo, es el int más grande posible. El comportamiento de RAND_MAX + 1 sería indefinido en casos en los que desbordaría 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 implementación de std::rand que he visto no lo hace ) introducir ningún sesgo estadístico adicional más allá de la alternativa presentada. Tales casos en los que % 6 es peligroso son aquellos en los que el generador de números tiene correlaciones simples en los bits de orden bajo, como una implementación de IBM bastante famosa (en C) de rand en, creo, en la década de 1970 que volcó 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, lo que probablemente no sea.

En conclusión, en estos días, debido a su facilidad de uso, 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 apropiadas para su caso de uso.

Answer 2

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

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

Segundo, si el número de valores distintos producidos por rand() no es un múltiplo de 6, entonces 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 para 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() % 6 devuelve 0, como si rand() hubiera devuelto 0. Entonces obtienes una distribución con el doble de 0 que cualquier otro valor.

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

La forma de evitar ese problema es descartar 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 cada vez que rand() devuelve un valor mayor o igual a ese múltiplo, lo rechaza y llama a `rand () nuevamente, tantas veces como sea necesario.

Asi que:

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, con la intención de mostrar más claramente lo que está sucediendo.

Answer 3

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) son menos parciales que 1+std::rand()%6 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 verifíquelo). Un enlace para ejecutar el código se encuentra here . 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 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 intercambio de pila entra en más detalles sobre el uso de la prueba de chi-cuadrado para evaluar la imparcialidad del dado: ¿Cómo puedo probar si un dado es justo? . Aquí están los resultados para 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 ejecuciones que hice, el valor p para ambos métodos siempre fue mayor que los valores alfa típicos utilizados para evaluar la significancia (0.05). Esto significa que no consideraríamos que ninguno de ellos sea 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 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í. 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 obtendrá una variedad de resultados. En segundo lugar, si no cambia la semilla, pero cambia el número de ensayos, también obtiene una variedad de resultados. Intenta aumentar o disminuir en un orden de magnitud para ver a qué me refiero. En tercer lugar, hay algún truncamiento o redondeo de enteros donde los valores esperados no son del todo precisos. Probablemente no sea suficiente para marcar la diferencia, pero está ahí.

Básicamente, en resumen, obtuvo la semilla correcta y la cantidad de pruebas que podría estar obteniendo un resultado falso.

Answer 4

Uno puede pensar en un generador de números aleatorios como trabajando en una secuencia de dígitos binarios. El generador convierte el flujo en números al dividirlo en trozos. Si la función std:rand funciona con un RAND_MAX de 32767, entonces está utilizando 15 bits en cada segmento.

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

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.

Uno espera que cada bit en el flujo de bits tenga la misma posibilidad 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 implementaciones de PRNG de software ofrecen diferentes compromisos entre velocidad y calidad. Un generador congruencial 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 a la velocidad más baja.

Answer 5

Este código de ejemplo ilustra que std::rand es un caso de balderdash de culto de carga heredado que debería hacer que sus cejas se alcen cada vez que lo vea.

Hay varios problemas aqui:

El contrato que la gente suele asumir, incluso las pobres almas desafortunadas que no conocen mejor y no pensarán en eso precisamente en estos términos, es que las muestras rand de la distribución uniforme en los 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 independientes en cada llamada, en realidad no es lo que dice la documentación, y en la práctica, las implementaciones históricamente no lograron proporcionar el más mínimo simulacro de independencia. Por ejemplo, C99 §7.20.2.1 ''La función rand '' dice, sin más detalles:

La función rand 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 ISO abusen del lenguaje. Después de todo, los únicos lectores que estarían molestos por eso 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 , 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

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

la expresión (a & 1) ^ (b & 1) produce 1 con 100% de probabilidad, que no es el caso para muestras aleatorias independientes en cualquier distribución compatible con enteros pares e impares. Por lo tanto, surgió un culto a la carga de que uno debería descartar los bits de bajo orden para perseguir a la bestia esquiva de ''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 qué están hablando o piensa que no tienes ni idea y que debes ser condescendiente).

El segundo problema es que incluso si cada llamada muestreara 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 dado, a menos que RAND_MAX sea ​​congruente con -1 módulo 6. Contraejemplo simple: si RAND_MAX = 6, entonces de rand() , todos los resultados tienen la misma probabilidad 1/7, pero de rand() % 6 , el el resultado 0 tiene una probabilidad de 2/7 mientras que todos los demás resultados tienen una probabilidad de 1/7.

La forma correcta de hacerlo es con el muestreo de rechazo: extraiga repetidamente una muestra aleatoria uniforme e independiente de 0, 1, 2, ..., RAND_MAX y rechace (por ejemplo) los resultados 0, 1, 2, ..., ((RAND_MAX + 1) % 6) - 1 si obtiene uno de esos, comience de nuevo; de lo contrario, rendimiento 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 igualmente divisible por 6, y cada posible resultado de s % 6 se obtiene por el mismo número de resultados aceptados de rand() , por lo que si rand() se distribuye uniformemente entonces también lo es el s . No hay límite en el número de intentos, pero el número esperado es menor que 2, y la probabilidad de éxito crece exponencialmente con el número de intentos.

La elección de los resultados de rand() que rechaza es irrelevante, siempre que asigne un número igual de ellos a cada número entero a continuación 6. El código en cppreference.com toma una decisión diferente , debido al primer problema anterior: que nada es garantizado sobre la distribución o independencia de las salidas de rand() , y en la práctica, los bits de orden inferior exhiben patrones que no ''parecen lo suficientemente aleatorios'' (no importa que la próxima salida sea una función determinista de la anterior).

Ejercicio para el lector: demuestre que el código en cppreference.com produce una distribución uniforme en los dados si rand() produce una distribución uniforme en 0, 1, 2, ..., RAND_MAX .

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

Un tercer problema es que el espacio de semillas 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 parezca que no es así. al azar después de todo. Así que ni siquiera pienses en usar esto para la criptografía.

Puede ir a la ruta elegante de std::uniform_int_distribution clase std::uniform_int_distribution C ++ 11 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 , pero incluso eso no va a ser adecuado para generar material de clave criptográfica, y el tornado Mersenne también es un terrible cerdo espacial con un estado de varios kilobytes que causa estragos en el 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 lanzar dados de juguete como este ejemplo!

Otro enfoque es utilizar un generador de números pseudoaleatorio criptográfico simple con un estado pequeño, como un simple borrado de clave rápida 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 se ve comprometido alguna vez.