c++ - repetir - srand en c

¿Por qué la gente dice que hay un sesgo de módulo cuando se utiliza un generador de números aleatorios? (9)

He visto que esta pregunta es muy frecuente, pero nunca he visto una verdadera respuesta concreta. Así que voy a publicar uno aquí que espero ayude a la gente a entender por qué hay un "sesgo de módulo" cuando se usa un generador de números aleatorios, como rand() en C ++.

Definición

El sesgo de módulo es el sesgo inherente en el uso de la aritmética de módulo para reducir un conjunto de salida a un subconjunto del conjunto de entrada. En general, existe un sesgo cuando el mapeo entre el conjunto de entrada y salida no se distribuye equitativamente, como en el caso de usar aritmética de módulo cuando el tamaño del conjunto de salida no es un divisor del tamaño del conjunto de entrada.

Este sesgo es particularmente difícil de evitar en la computación, donde los números se representan como cadenas de bits: 0s y 1s. Encontrar fuentes de aleatoriedad verdaderamente aleatorias también es extremadamente difícil, pero está más allá del alcance de esta discusión. Para el resto de esta respuesta, suponga que existe una fuente ilimitada de bits verdaderamente aleatorios.

Ejemplo de problema

Consideremos simular una tirada (0 a 5) utilizando estos bits aleatorios. Hay 6 posibilidades, por lo que necesitamos suficientes bits para representar el número 6, que es de 3 bits. Desafortunadamente, 3 bits aleatorios producen 8 resultados posibles:

000 = 0, 001 = 1, 010 = 2, 011 = 3 100 = 4, 101 = 5, 110 = 6, 111 = 7

Podemos reducir el tamaño del conjunto de resultados a exactamente 6 tomando el valor módulo 6, sin embargo, esto presenta el problema de sesgo de módulo : 110 produce un 0 y 111 produce un 1. Este dado está cargado.

Soluciones potenciales

Enfoque 0:

En lugar de confiar en bits aleatorios, en teoría, se podría contratar un pequeño ejército para lanzar dados todo el día y registrar los resultados en una base de datos, y luego usar cada resultado solo una vez. Esto es tan práctico como parece, y es muy probable que de todos modos no produzca resultados verdaderamente aleatorios (intencionalmente).

Enfoque 1:

En lugar de usar el módulo, una solución ingenua pero matemáticamente correcta es descartar los resultados que arrojan 110 y 111 y simplemente intentar nuevamente con 3 bits nuevos. Desafortunadamente, esto significa que hay un 25% de probabilidad en cada tirada de que se requerirá una nueva tirada, incluyendo cada una de las tiradas . Esto es claramente impráctico para todos, excepto para los usos más triviales.

Enfoque 2:

Use más bits: en lugar de 3 bits, use 4. Esto produce 16 resultados posibles. Por supuesto, si se vuelve a tirar en cualquier momento, el resultado es mayor que 5, lo que empeora las cosas (10/16 = 62.5%), por lo que solo no ayudará.

Tenga en cuenta que 2 * 6 = 12 <16, por lo que podemos tomar de forma segura cualquier resultado inferior a 12 y reducir ese módulo 6 para distribuir uniformemente los resultados. Los otros 4 resultados se deben descartar y luego volver a tirar como en el enfoque anterior.

Suena bien al principio, pero revisemos las matemáticas:

4 discarded results / 16 possibilities = 25%

En este caso, 1 bit extra no ayudó en absoluto!

Ese resultado es desafortunado, pero intentemos nuevamente con 5 bits:

32 % 6 = 2 discarded results; and 2 discarded results / 32 possibilities = 6.25%

Una mejora definitiva, pero no lo suficientemente buena en muchos casos prácticos. La buena noticia es que agregar más bits nunca aumentará las posibilidades de tener que descartar y volver a tirar . Esto es válido no solo para los dados, sino en todos los casos.

Sin embargo, como se demostró , agregar 1 bit adicional puede no cambiar nada. De hecho, si aumentamos nuestra tirada a 6 bits, la probabilidad sigue siendo del 6.25%.

Esto plantea 2 preguntas adicionales:

Si agregamos suficientes bits, ¿existe una garantía de que la probabilidad de un descarte disminuya?
¿Cuántos bits son suficientes en el caso general?

Solución general

Afortunadamente, la respuesta a la primera pregunta es sí. El problema con 6 es que 2 ^ x mod 6 cambia entre 2 y 4, que casualmente son un múltiplo de 2 entre sí, de modo que para un x> 1 par,

[2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1)

Por lo tanto, 6 es una excepción y no la regla. Es posible encontrar módulos más grandes que produzcan potencias consecutivas de 2 de la misma manera, pero eventualmente esto debe envolver, y la probabilidad de un descarte se reducirá.

Sin ofrecer una prueba adicional, en general, utilizar el doble de bits requeridos proporcionará una posibilidad más pequeña, generalmente insignificante, de descartar.

Prueba de concepto

Aquí hay un programa de ejemplo que utiliza el libcrypo de OpenSSL para suministrar bytes aleatorios. Al compilar, asegúrese de enlazar a la biblioteca con -lcrypto que la mayoría de las personas deberían tener disponibles.

#include <iostream> #include <assert.h> #include <limits> #include <openssl/rand.h> volatile uint32_t dummy; uint64_t discardCount; uint32_t uniformRandomUint32(uint32_t upperBound) { assert(RAND_status() == 1); uint64_t discard = (std::numeric_limits<uint64_t>::max() - upperBound) % upperBound; uint64_t randomPool = RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool)); while(randomPool > (std::numeric_limits<uint64_t>::max() - discard)) { RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool)); ++discardCount; } return randomPool % upperBound; } int main() { discardCount = 0; const uint32_t MODULUS = (1ul << 31)-1; const uint32_t ROLLS = 10000000; for(uint32_t i = 0; i < ROLLS; ++i) { dummy = uniformRandomUint32(MODULUS); } std::cout << "Discard count = " << discardCount << std::endl; }

Recomiendo jugar con los valores de MODULUS y ROLLS para ver cuántas repeticiones ocurren en la mayoría de las condiciones. Una persona escéptica también puede desear guardar los valores computados para archivar y verificar que la distribución parece normal.

@ user1413793 es correcto sobre el problema. No lo discutiré más a fondo, excepto para hacer un punto: sí, para valores pequeños de n y valores grandes de RAND_MAX , el sesgo de módulo puede ser muy pequeño. Pero el uso de un patrón inductor de sesgo significa que debe considerar el sesgo cada vez que calcula un número aleatorio y elige patrones diferentes para diferentes casos. Y si toma la decisión equivocada, los errores que introduce son sutiles y casi imposibles de probar por unidad. Comparado con solo usar la herramienta adecuada (como arc4random_uniform ), eso es trabajo extra, no menos trabajo. Hacer más trabajo y obtener una solución peor es una ingeniería terrible, especialmente cuando hacerlo bien es fácil en la mayoría de las plataformas.

Desafortunadamente, las implementaciones de la solución son todas incorrectas o menos eficientes de lo que deberían ser. (Cada solución tiene varios comentarios que explican los problemas, pero ninguna de las soluciones ha sido solucionada para solucionarlos). Es probable que esto confunda al buscador casual de respuestas, por lo que aquí proporciono una buena aplicación.

Una vez más, la mejor solución es usar arc4random_uniform en las plataformas que lo proporcionan, o una solución de rango similar para su plataforma (como Random.nextInt en Java). Hará lo correcto sin costo para usted. Esta es casi siempre la llamada correcta para hacer.

Si no tiene arc4random_uniform , entonces puede usar el poder de opensource para ver exactamente cómo se implementa sobre un RNG de rango más amplio ( ar4random en este caso, pero un enfoque similar también podría funcionar sobre otros RNG) .

Aquí está la implementación de OpenBSD :

/* * Calculate a uniformly distributed random number less than upper_bound * avoiding "modulo bias". * * Uniformity is achieved by generating new random numbers until the one * returned is outside the range [0, 2**32 % upper_bound). This * guarantees the selected random number will be inside * [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound) * after reduction modulo upper_bound. */ u_int32_t arc4random_uniform(u_int32_t upper_bound) { u_int32_t r, min; if (upper_bound < 2) return 0; /* 2**32 % x == (2**32 - x) % x */ min = -upper_bound % upper_bound; /* * This could theoretically loop forever but each retry has * p > 0.5 (worst case, usually far better) of selecting a * number inside the range we need, so it should rarely need * to re-roll. */ for (;;) { r = arc4random(); if (r >= min) break; } return r % upper_bound; }

Vale la pena señalar el último comentario de confirmación sobre este código para aquellos que necesitan implementar cosas similares:

Cambie arc4random_uniform () para calcular 2**32 % upper_bound'''' as -upper_bound% upper_bound ''''. Simplifica el código y lo hace igual en las arquitecturas ILP32 y LP64, y también un poco más rápido en las arquitecturas LP64 mediante el uso de un resto de 32 bits en lugar de un resto de 64 bits.
Destacado por Jorden Verwer en tech @ ok deraadt; no hay objeciones de djm u otto

La implementación de Java también es fácil de encontrar (ver enlace anterior):

public int nextInt(int n) { if (n <= 0) throw new IllegalArgumentException("n must be positive"); if ((n & -n) == n) // i.e., n is a power of 2 return (int)((n * (long)next(31)) >> 31); int bits, val; do { bits = next(31); val = bits % n; } while (bits - val + (n-1) < 0); return val; }

Acabo de escribir un código para el Método de cambio de moneda imparcial de Von Neumann, que en teoría debería eliminar cualquier sesgo en el proceso de generación de números aleatorios. Se puede encontrar más información en ( http://en.wikipedia.org/wiki/Fair_coin )

int unbiased_random_bit() { int x1, x2, prev; prev = 2; x1 = rand() % 2; x2 = rand() % 2; for (;; x1 = rand() % 2, x2 = rand() % 2) { if (x1 ^ x2) // 01 -> 1, or 10 -> 0. { return x2; } else if (x1 & x2) { if (!prev) // 0011 return 1; else prev = 1; // 1111 -> continue, bias unresolved } else { if (prev == 1)// 1100 return 0; else // 0000 -> continue, bias unresolved prev = 0; } } }

Como indica la respuesta aceptada , el "sesgo de módulo" tiene sus raíces en el bajo valor de RAND_MAX . Utiliza un valor extremadamente pequeño de RAND_MAX (10) para mostrar que si RAND_MAX fuera 10, entonces intentó generar un número entre 0 y 2 usando%, los siguientes resultados resultarían:

rand() % 3 // if RAND_MAX were only 10, gives output of rand() | rand()%3 0 | 0 1 | 1 2 | 2 3 | 0 4 | 1 5 | 2 6 | 0 7 | 1 8 | 2 9 | 0

Así que hay 4 salidas de 0 (4/10 de probabilidad) y solo 3 salidas de 1 y 2 (3/10 de probabilidad cada una).

Así que es parcial. Los números más bajos tienen una mejor probabilidad de salir.

Pero eso solo aparece tan obviamente cuando RAND_MAX es pequeño . O más específicamente, cuando el número por el que está modificando es grande en comparación con RAND_MAX .

Una solución mucho mejor que el bucle (que es increíblemente ineficiente y que ni siquiera debería sugerirse) es usar un PRNG con un rango de salida mucho mayor. El algoritmo Mersenne Twister tiene una salida máxima de 4,294,967,295. Como tal, MersenneTwister::genrand_int32() % 10 para todos los propósitos y propósitos, se distribuirá por igual y el efecto de sesgo de módulo desaparecerá.

Con un valor RAND_MAX de 3 (en realidad debería ser mucho más alto que eso, pero el sesgo seguiría existiendo), a partir de estos cálculos, tiene sentido que exista un sesgo:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

En este caso, el % 2 es lo que no debe hacer cuando desea un número aleatorio entre 0 y 1 . Sin embargo, podría obtener un número aleatorio entre 0 y 2 haciendo % 3 , porque en este caso: RAND_MAX es un múltiplo de 3 .

Otro método

Es mucho más simple, pero para agregar a otras respuestas, aquí está mi solución para obtener un número aleatorio entre 0 y n - 1 , así que n diferentes posibilidades, sin sesgo.

el número de bits (no bytes) necesarios para codificar el número de posibilidades es el número de bits de datos aleatorios que necesitará
codificar el número de bits aleatorios
Si este número es >= n , reinicie (sin módulo).

Los datos realmente aleatorios no son fáciles de obtener, entonces, ¿por qué usar más bits de los necesarios?

A continuación se muestra un ejemplo en Smalltalk, que utiliza un caché de bits de un generador de números pseudoaleatorios. No soy un experto en seguridad, así que lo uso bajo su propio riesgo.

next: n | bitSize r from to | n < 0 ifTrue: [^0 - (self next: 0 - n)]. n = 0 ifTrue: [^nil]. n = 1 ifTrue: [^0]. cache isNil ifTrue: [cache := OrderedCollection new]. cache size < (self randmax highBit) ifTrue: [ Security.DSSRandom default next asByteArray do: [ :byte | (1 to: 8) do: [ :i | cache add: (byte bitAt: i)] ] ]. r := 0. bitSize := n highBit. to := cache size. from := to - bitSize + 1. (from to: to) do: [ :i | r := r bitAt: i - from + 1 put: (cache at: i) ]. cache removeFrom: from to: to. r >= n ifTrue: [^self next: n]. ^r

Entonces rand() es un generador de números pseudoaleatorios que elige un número natural entre 0 y RAND_MAX , que es una constante definida en cstdlib (consulte este article para obtener una descripción general sobre rand() ).

Ahora, ¿qué pasa si quieres generar un número aleatorio entre 0 y 2? En aras de la explicación, digamos que RAND_MAX es 10 y decido generar un número aleatorio entre 0 y 2 llamando a rand()%3 . Sin embargo, ¡ rand()%3 no produce los números entre 0 y 2 con igual probabilidad!

Cuando rand() devuelve 0, 3, 6 o 9, rand()%3 == 0 . Por lo tanto, P (0) = 4/11

Cuando rand() devuelve 1, 4, 7 o 10, rand()%3 == 1 . Por lo tanto, P (1) = 4/11

Cuando rand() devuelve 2, 5 u 8, rand()%3 == 2 . Por lo tanto, P (2) = 3/11

Esto no genera los números entre 0 y 2 con igual probabilidad. Por supuesto, para rangos pequeños, este podría no ser el mayor problema, pero para un rango más grande podría distorsionar la distribución, sesgando los números más pequeños.

Entonces, ¿cuándo rand()%n devuelve un rango de números de 0 a n-1 con igual probabilidad? Cuando RAND_MAX%n == n - 1 . En este caso, junto con nuestro supuesto anterior rand() devuelve un número entre 0 y RAND_MAX con igual probabilidad, las clases de módulo de n también se distribuirían por igual.

Entonces, ¿cómo resolvemos este problema? Una forma cruda es seguir generando números aleatorios hasta que obtenga un número en el rango deseado:

int x; do { x = rand(); } while (x >= n);

pero eso es ineficiente para los valores bajos de n , ya que solo tiene una probabilidad n/RAND_MAX de obtener un valor en su rango, por lo que deberá realizar llamadas RAND_MAX/n a rand() en promedio.

Un enfoque de fórmula más eficiente sería tomar un rango grande con una longitud divisible por n , como RAND_MAX - RAND_MAX % n , seguir generando números aleatorios hasta que obtenga uno que se encuentre dentro del rango y luego tomar el módulo:

int x; do { x = rand(); } while (x >= (RAND_MAX - RAND_MAX % n)); x %= n;

Para valores pequeños de n , esto rara vez requerirá más de una llamada a rand() .

Obras citadas y lecturas adicionales:

Hay dos quejas habituales con el uso del módulo.

Uno es válido para todos los generadores. Es más fácil de ver en un caso límite. Si su generador tiene un RAND_MAX que es 2 (que no cumple con el estándar C) y desea solo 0 o 1 como valor, el uso de módulo generará 0 dos veces más a menudo (cuando el generador genere 0 y 2) como lo hará generar 1 (cuando el generador genera 1). Tenga en cuenta que esto es así tan pronto como no elimine los valores, independientemente de la asignación que esté utilizando desde los valores del generador hasta el deseado, uno ocurrirá con el doble de frecuencia que el otro.
algún tipo de generador tiene sus bits menos significativos menos aleatorios que el otro, al menos para algunos de sus parámetros, pero lamentablemente esos parámetros tienen otras características interesantes (como tener RAND_MAX uno menos que una potencia de 2). El problema es bien conocido y durante mucho tiempo la implementación de la biblioteca probablemente evite el problema (por ejemplo, la implementación rand () de ejemplo en el estándar C utiliza este tipo de generador, pero elimina los 16 bits menos significativos), pero a algunos les gusta quejarse eso y puede que tengas mala suerte

Usando algo como

int alea(int n){ assert (0 < n && n <= RAND_MAX); int partSize = n == RAND_MAX ? 1 : 1 + (RAND_MAX-n)/(n+1); int maxUsefull = partSize * n + (partSize-1); int draw; do { draw = rand(); } while (draw > maxUsefull); return draw/partSize; }

para generar un número aleatorio entre 0 y n evitará ambos problemas (y evita el desbordamiento con RAND_MAX == INT_MAX)

Por cierto, C ++ 11 introdujo formas estándar para la reducción y otro generador que rand ().

La solución de Mark (la solución aceptada) es casi perfecta.

int x; do { x = rand(); } while (x >= (RAND_MAX - RAND_MAX % n)); x %= n;
Editado el 25 de marzo de 2016 a las 23:16.
Mark Amery 39k21170211

Sin embargo, tiene una advertencia que descarta 1 conjunto válido de resultados en cualquier escenario donde RAND_MAX (RM) es 1 menos que un múltiplo de N (donde N = el número de posibles resultados válidos).

es decir, cuando el ''recuento de valores descartados'' (D) es igual a N, en realidad son un conjunto válido (V), no un conjunto no válido (I).

Usando la solución de Mark, los valores se descartan cuando: X => RM - RM% N

EG: Ran Max Value (RM) = 255 Valid Outcome (N) = 4 When X => 252, Discarded values for X are: 252, 253, 254, 255 So, if Random Value Selected (X) = {252, 253, 254, 255} Number of discarded Values (I) = RM % N + 1 == N IE: I = RM % N + 1 I = 255 % 4 + 1 I = 3 + 1 I = 4 X => ( RM - RM % N ) 255 => (255 - 255 % 4) 255 => (255 - 3) 255 => (252) Discard Returns $True

Como puede ver en el ejemplo anterior, cuando el valor de X (el número aleatorio que obtenemos de la función inicial) es 252, 253, 254 o 255, lo descartaremos incluso si estos cuatro valores comprenden un conjunto válido de valores devueltos .

IE: Cuando el recuento de los valores descartados (I) = N (el número de resultados válidos), la función original descartará un conjunto válido de valores de retorno.

Si describimos la diferencia entre los valores N y RM como D, es decir:

D = (RM - N)

Luego, a medida que el valor de D se hace más pequeño, el Porcentaje de repeticiones innecesarias debidas a este método aumenta con cada multiplicativo natural. (Cuando RAND_MAX NO es igual a un Número Prime, esto es de preocupación válida)

P.EJ:

RM=255 , N=2 Then: D = 253, Lost percentage = 0.78125% RM=255 , N=4 Then: D = 251, Lost percentage = 1.5625% RM=255 , N=8 Then: D = 247, Lost percentage = 3.125% RM=255 , N=16 Then: D = 239, Lost percentage = 6.25% RM=255 , N=32 Then: D = 223, Lost percentage = 12.5% RM=255 , N=64 Then: D = 191, Lost percentage = 25% RM=255 , N= 128 Then D = 127, Lost percentage = 50%

Dado que el porcentaje de Rerolls necesarios aumenta a medida que N se acerca más a RM, esto puede ser una preocupación válida en muchos valores diferentes dependiendo de las restricciones del sistema que ejecuta el código y los valores que se buscan.

Para negar esto podemos hacer una enmienda simple como se muestra aquí:

int x; do { x = rand(); } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ); x %= n;

Esto proporciona una versión más general de la fórmula que explica las peculiaridades adicionales de usar el módulo para definir sus valores máximos.

Ejemplos de uso de un valor pequeño para RAND_MAX que es un multiplicativo de N.

Versión Mark''original:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3. When X >= (RAND_MAX - ( RAND_MAX % n ) ) When X >= 2 the value will be discarded, even though the set is valid.

Versión generalizada 1:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3. When X > (RAND_MAX - ( ( RAND_MAX % n ) + 1 ) % n ) When X > 3 the value would be discarded, but this is not a vlue in the set RAND_MAX so there will be no discard.

Además, en el caso donde N debería ser el número de valores en RAND_MAX; en este caso, podría establecer N = RAND_MAX +1, a menos que RAND_MAX = INT_MAX.

En lo que respecta a los bucles, puede usar N = 1, y se aceptará cualquier valor de X, sin embargo, y poner una declaración IF en su multiplicador final. Pero quizás tenga un código que puede tener una razón válida para devolver un 1 cuando se llama a la función con n = 1 ...

Por lo tanto, puede ser mejor usar 0, que normalmente proporcionaría un error Div 0, cuando desea tener n = RAND_MAX + 1

Versión generalizada 2:

int x; if n != 0 { do { x = rand(); } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ); x %= n; } else { x = rand(); }

Ambas soluciones resuelven el problema con resultados válidos descartados innecesariamente que ocurrirán cuando RM + 1 sea un producto de n.

La segunda versión también cubre el escenario de caso extremo cuando necesita n para igualar el conjunto total posible de valores contenidos en RAND_MAX.

El enfoque modificado en ambos es el mismo y permite una solución más general a la necesidad de proporcionar números aleatorios válidos y minimizar los valores descartados.

Reiterar:

La Solución General Básica que amplía el ejemplo de la marca:

int x; do { x = rand(); } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ); x %= n;

La solución general extendida que permite un escenario adicional de RAND_MAX + 1 = n:

int x; if n != 0 { do { x = rand(); } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) ); x %= n; } else { x = rand(); }

Seguir seleccionando un azar es una buena manera de eliminar el sesgo.

Actualizar

Podríamos hacer el código rápido si buscamos una x en el rango divisible por n .

// Assumptions // rand() in [0, RAND_MAX] // n in (0, RAND_MAX] int x; // Keep searching for an x in a range divisible by n do { x = rand(); } while (x >= RAND_MAX - (RAND_MAX % n)) x %= n;

El bucle anterior debe ser muy rápido, digamos 1 iteración en promedio.