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

277

He visto muchas veces esta pregunta pero nunca he visto una respuesta concreta verdadera. Así que voy a publicar uno aquí que, con suerte, ayudará a las personas a entender por qué exactamente hay un "sesgo de módulo" cuando se usa un generador de números aleatorios, como rand()en C ++.

usuario1413793
fuente

Respuestas:

394

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

¿Qué sucede si quieres generar un número aleatorio entre digamos 0 y 2? En aras de la explicación, digamos que RAND_MAXes 10 y decido generar un número aleatorio entre 0 y 2 llamando rand()%3. ¡Sin embargo, rand()%3no produce los números entre 0 y 2 con la misma 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 sesgar la distribución, sesgando los números más pequeños.

Entonces, ¿cuándo rand()%ndevuelve 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_MAXcon 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 valores bajos de n, ya que solo tiene la n/RAND_MAXposibilidad de obtener un valor en su rango, por lo que deberá realizar RAND_MAX/nllamadas rand()en promedio.

Un enfoque de fórmula más eficiente sería tomar un rango grande con una longitud divisible por n, por ejemplo RAND_MAX - RAND_MAX % n, seguir generando números aleatorios hasta obtener uno que se encuentre en el 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:


usuario1413793
fuente
66
Otra forma de pensar sobre RAND_MAX%n == n - 1_ es (RAND_MAX + 1) % n == 0. Cuando leo el código, tiendo a entenderlo % something == 0como "divisible por igual" más fácilmente que otras formas de calcularlo. Por supuesto, si su C ++ stdlib tiene RAND_MAXel mismo valor que INT_MAX, (RAND_MAX + 1)seguramente no funcionaría; entonces el cálculo de Mark sigue siendo la implementación más segura.
Slipp D. Thompson
muy buena respuesta!
Sayali Sonawane
Es posible que esté haciendo trampas, pero si el objetivo es reducir los bits desperdiciados, podríamos mejorar esto ligeramente para la condición de borde donde RAND_MAX (RM) es solo 1 menos que ser igualmente divisible por N. En este escenario, no es necesario desperdiciar ningún bit por haciendo X> = (RM - RM% N)) que es de poco valor para valores pequeños de N, pero se vuelve de mayor valor para valores grandes de N. Como se menciona por Slipp D. Thompson, hay una solución que funcionará solo cuando INT_MAX (IM)> RAND_MAX pero se rompe cuando son iguales. Sin embargo, hay una solución simple para esto: podemos modificar el cálculo X> = (RM - RM% N) de la siguiente manera:
Ben Personick el
X> = RM - (((RM% N) + 1)% N)
Ben Personick
Publiqué una respuesta adicional explicando el problema en detalle y dando la solución de código de ejemplo.
Ben Personick
36

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

Actualizar

Podríamos hacer que el código sea 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.

Nick Dandoulakis
fuente
2
Yuck :-P convertir a un doble, luego multiplicar por MAX_UPPER_LIMIT / RAND_MAX es mucho más limpio y funciona mejor.
boycy
22
@boycy: te has perdido el punto. Si el número de valores que rand()puede devolver no es un múltiplo de n, entonces, haga lo que haga, inevitablemente obtendrá un "sesgo de módulo", a menos que descarte algunos de esos valores. user1413793 lo explica muy bien (aunque la solución propuesta en esa respuesta es realmente asquerosa).
TonyK
44
@TonyK mis disculpas, no entendí el punto. No pensé lo suficiente y pensé que el sesgo solo se aplicaría con métodos que usan una operación de módulo explícito. Gracias por arreglarme :-)
boycy
La precedencia del operador hace que el RAND_MAX+1 - (RAND_MAX+1) % ntrabajo funcione correctamente, pero sigo pensando que debería escribirse RAND_MAX+1 - ((RAND_MAX+1) % n)en aras de la claridad.
Linus Arver
44
Esto no funcionará si RAND_MAX == INT_MAX (como lo hace en la mayoría de los sistemas) . Vea mi segundo comentario a @ user1413793 arriba.
BlueRaja - Danny Pflughoeft
19

@ user1413793 es correcto sobre el problema. No voy a discutir eso más a fondo, excepto para hacer un punto: sí, para valores pequeños de ny 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 calcule un número aleatorio y elija diferentes patrones para diferentes casos. Y si toma la decisión equivocada, los errores que presenta son sutiles y casi imposibles de probar. En comparación 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 cada vez 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 se ha solucionado para abordarlos). Es probable que esto confunda al buscador informal de respuestas, por lo que estoy proporcionando una implementación bien conocida aquí.

Una vez más, la mejor solución es usarlo arc4random_uniformen plataformas que lo proporcionan, o una solución a distancia similar para su plataforma (como Random.nextInten Java). Hará lo correcto sin costo para usted. Esta es casi siempre la llamada correcta para hacer.

Si no lo tiene arc4random_uniform, puede usar el poder de código abierto para ver exactamente cómo se implementa en la parte superior de un RNG de mayor alcance ( ar4randomen este caso, pero un enfoque similar también podría funcionar en la parte superior de 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_boundcomo -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.

Señalado por Jorden Verwer en tech @ ok deraadt; sin objeciones de djm u otto

La implementación de Java también se puede encontrar fácilmente (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;
 }
Rob Napier
fuente
Tenga en cuenta que si arcfour_random() realmente utiliza el algoritmo RC4 real en su implementación, la salida definitivamente tendrá algún sesgo. Esperemos que los autores de su biblioteca hayan cambiado a utilizar un mejor CSPRNG detrás de la misma interfaz. Recuerdo que uno de los BSD ahora usa el algoritmo ChaCha20 para implementar arcfour_random(). Más sobre los sesgos de salida de RC4 que lo hacen inútil para la seguridad u otras aplicaciones críticas como el video póker: blog.cryptographyengineering.com/2013/03/…
rmalayter
2
@rmalayter En iOS y OS X, arc4random lee desde / dev / random, que es la entropía de mayor calidad en el sistema. (El "arc4" en el nombre es histórico y se conserva por compatibilidad.)
Rob Napier
@Rob_Napier es bueno saberlo, pero /dev/randomtambién ha usado RC4 en algunas plataformas en el pasado (Linux usa SHA-1 en modo contador). Desafortunadamente, las páginas de manual que encontré a través de la búsqueda indican que RC4 todavía está en uso en varias plataformas que ofrecen arc4random(aunque el código real puede ser diferente).
rmalayter
1
Estoy confundido. No es -upper_bound % upper_bound == 0??
Jon McClung
1
@JonMcClung -upper_bound % upper_boundserá de hecho 0 si intes más ancho que 32 bits. Debería serlo (u_int32_t)-upper_bound % upper_bound)(suponiendo que u_int32_tes un BSD-ismo para uint32_t).
Ian Abbott
14

Definición

Modulo Bias es el sesgo inherente en el uso del módulo aritmético para reducir un conjunto de salida a un subconjunto del conjunto de entrada. En general, existe un sesgo siempre que el mapeo entre el conjunto de entrada y salida no se distribuye por igual, como en el caso de usar módulo aritmético 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 informática, donde los números se representan como cadenas de bits: 0s y 1s. Encontrar fuentes verdaderamente aleatorias de aleatoriedad 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 de dado (0 a 5) usando 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 exactamente a 6 tomando el valor módulo 6, sin embargo, esto presenta el problema de sesgo de módulo : 110produce un 0 y 111produce un 1. Este dado está cargado.

Soluciones potenciales

Enfoque 0:

En lugar de confiar en bits aleatorios, en teoría uno podría contratar a 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 lo más probable es que no produzca resultados verdaderamente aleatorios de todos modos (juego de palabras).

Enfoque 1:

En lugar de utilizar el módulo, una solución ingenua pero matemáticamente correcto es a los resultados de descarte que el rendimiento 110y 111y simplemente tratar de nuevo con 3 nuevos bits. Desafortunadamente, esto significa que hay un 25% de posibilidades en cada lanzamiento de que se requerirá un nuevo lanzamiento, incluyendo cada uno de los mismos. Esto es claramente poco prá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, volver a tirar cada vez que el resultado sea mayor que 5 empeora las cosas (10/16 = 62.5%) para que eso no ayude.

Observe que 2 * 6 = 12 <16, por lo que podemos tomar con seguridad cualquier resultado inferior a 12 y reducir ese módulo 6 para distribuir los resultados de manera uniforme. Los otros 4 resultados deben descartarse y luego volverse 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 extra puede no cambiar nada. De hecho, si aumentamos nuestro rollo a 6 bits, la probabilidad sigue siendo 6.25%.

Esto plantea 2 preguntas adicionales:

  1. Si agregamos suficientes bits, ¿hay una garantía de que disminuirá la probabilidad de un descarte?
  2. ¿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 múltiplos de 2 entre sí, de modo que para un par x> 1,

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

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

Sin ofrecer pruebas adicionales, en general, usar el doble del número de bits requeridos proporcionará una posibilidad menor, generalmente insignificante, de descarte.

Prueba de concepto

Aquí hay un programa de ejemplo que usa libcrypo de OpenSSL para suministrar bytes aleatorios. Al compilar, asegúrese de vincular a la biblioteca con la -lcryptoque casi todos deberían estar 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;
}

Animo a jugar con los valores MODULUSy ROLLSpara ver cuántas repeticiones ocurren en la mayoría de las condiciones. Una persona escéptica también puede desear guardar los valores calculados en el archivo y verificar que la distribución parezca normal.

Jim Wood
fuente
Realmente espero que nadie haya copiado ciegamente su implementación aleatoria uniforme. La randomPool = RAND_bytes(...)línea siempre resultará randomPool == 1debido a la aserción. Esto siempre da como resultado un descarte y una repetición. Creo que querías declarar en una línea separada. En consecuencia, esto hizo que el RNG volviera con 1cada iteración.
Qix - MONICA FUE MAL
Para ser claros, randomPoolsiempre se evaluará de 1acuerdo con la documentación deRAND_bytes() OpenSSL ya que siempre tendrá éxito gracias a la RAND_status()afirmación.
Qix - MONICA FUE MAL
9

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 del módulo generará 0 el doble de veces (cuando el generador genera 0 y 2) como lo hará generar 1 (cuando el generador genera 1). Tenga en cuenta que esto es cierto tan pronto como no elimine los valores, cualquiera que sea la asignación que esté utilizando desde los valores del generador hasta el deseado, uno ocurrirá el doble de veces 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 lo evitará (por ejemplo, la implementación de rand () de muestra en el estándar C usa este tipo de generador, pero descarta 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; 
}

generar un número aleatorio entre 0 yn 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 no sea rand ().

Un programador
fuente
n == RAND_MAX? 1: (RAND_MAX-1) / (n + 1): Entiendo que la idea aquí es dividir primero RAND_MAX en el tamaño de página N igual, luego devolver la desviación dentro de N, pero no puedo asignar el código a esto con precisión.
Zinking
1
La versión ingenua debe ser (RAND_MAX + 1) / (n + 1) ya que hay valores RAND_MAX + 1 para dividir en n + 1 cubos. Si se desea evitar el desbordamiento al calcular RAND_MAX + 1, se puede transformar en 1+ (RAND_MAX-n) / (n + 1). Para evitar el desbordamiento al calcular n + 1, primero se verifica el caso n == RAND_MAX.
Programador
+ más, hacer dividir parece costar más incluso en comparación con los números de regeneración.
Zinking
44
Tomar el módulo y dividir tiene el mismo costo. Algunos ISA incluso proporcionan solo una instrucción que proporciona siempre ambas. El costo de regenerar números dependerá de n y RAND_MAX. Si n es pequeño con respecto a RAND_MAX, puede costar mucho. Y, obviamente, puede decidir que los sesgos no son importantes para su aplicación; Solo doy una manera de evitarlos.
Programador
9

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 16 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 menor que un múltiplo de N(Donde N= el Número de posibles resultados válidos).

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

Lo que causa esto es que en algún momento Mark pierde de vista la diferencia entre Ny Rand_Max.

Nes un conjunto cuyos miembros válidos están compuestos solo por números enteros positivos, ya que contiene un recuento de respuestas que serían válidas. (por ejemplo: Set N= {1, 2, 3, ... n })

Rand_max Sin embargo, es un conjunto que (como se define para nuestros propósitos) incluye cualquier número de enteros no negativos.

En su forma más genérica, lo que se define aquí Rand Maxes el Conjunto de todos los resultados válidos, que teóricamente podría incluir números negativos o valores no numéricos.

Por Rand_Maxlo tanto, se define mejor como el conjunto de "Posibles respuestas".

Sin embargo, Nopera contra el recuento de los valores dentro del conjunto de respuestas válidas, por lo que incluso según lo definido en nuestro caso específico, Rand_Maxserá un valor uno menos que el número total que contiene.

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 descartaríamos aunque estos cuatro valores comprendan 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 vuelve más pequeño, el porcentaje de repeticiones innecesarias debido a este método aumenta en cada multiplicativo natural. (Cuando RAND_MAX NO es igual a un número primo, esto es una 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 necesario aumenta cuanto más se acerca N 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 del uso del 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 original de Mark:

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, puede establecer N = RAND_MAX +1, a menos que RAND_MAX = INT_MAX.

En cuanto al bucle, podría usar N = 1, y cualquier valor de X será aceptado, sin embargo, y colocará una declaración IF para su multiplicador final. Pero quizás tenga un código que pueda tener una razón válida para devolver un 1 cuando la función se llama con n = 1 ...

Por lo tanto, puede ser mejor usar 0, que normalmente proporcionaría un error Div 0, cuando desee 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 límite 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 mark:

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

 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:

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

int x;

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

    x %= n;
} else {
    x = rand();
}

En algunos idiomas (particularmente los idiomas interpretados) hacer los cálculos de la operación de comparación fuera de la condición while puede conducir a resultados más rápidos, ya que este es un cálculo único, sin importar cuántas repeticiones se requieran. YMMV!

// Assumes:
//  RAND_MAX is a globally defined constant, returned from the environment.
//  int n; // User input, or externally defined, number of valid choices.

int x; // Resulting random number
int y; // One-time calculation of the compare value for x

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

    x %= n;
} else {
    x = rand();
}
Ben Personick
fuente
¿No es seguro decir que el problema con la solución de Mark es que trata a RAND_MAX yn como la misma "unidad de medida" cuando en realidad significan dos cosas diferentes? Mientras que n representa el "número de posibilidades" resultante, RAND_MAX solo representa el valor máximo de la posibilidad original, donde RAND_MAX + 1 sería el número original de posibilidades. Me sorprende que no haya llegado a su conclusión, ya que parecía haber reconocido que ny RAND_MAX no eran lo mismo con la ecuación:RAND_MAX%n = n - 1
Danilo Souza Morães
@ DaniloSouzaMorães Gracias Danilo, has planteado el asunto de manera muy sucinta. Fui a demostrar lo que estaba haciendo junto con el por qué y cómo, pero no creo que alguna vez fui capaz de decir QUÉ estaba haciendo mal elocuentemente, ya que estoy tan absorto en los detalles de la lógica sobre cómo y por qué hay un problema, que no estoy declarando tan claramente lo que está en cuestión. ¿Le importa si modifico mi Respuesta para usar algo de lo que ha escrito aquí como mi propio resumen sobre la cuestión de qué y dónde está haciendo la solución aceptada lo que debe abordarse en la parte superior?
Ben Personick
Que sería increíble. Anímate
Danilo Souza Morães
1

Con un RAND_MAXvalor de 3(en realidad debería ser mucho más alto que eso, pero el sesgo aún existiría) tiene sentido a partir de estos cálculos que existe un sesgo:

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

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

Otro método

Hay mucho más simple pero para agregar a otras respuestas, aquí está mi solución para obtener un número aleatorio entre 0y n - 1, por lo tanto, ndiferentes posibilidades, sin sesgos.

  • la cantidad de bits (no bytes) necesarios para codificar la cantidad de posibilidades es la cantidad 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 pseudoaleatorio. No soy un experto en seguridad, así que úselo 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
Rivenfall
fuente
-1

Como indica la respuesta aceptada , "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 trataste de generar un número entre 0 y 2 usando%, se obtendrían los siguientes resultados:

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

Por lo tanto, hay 4 salidas de 0 (probabilidad de 4/10) y solo 3 salidas de 1 y 2 (3/10 posibilidades cada una).

Entonces es parcial. Los números más bajos tienen una mejor oportunidad de salir.

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

Una solución mucho mejor que el bucle (que es increíblemente ineficiente y 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() % 10para todos los efectos, se distribuirá por igual y el efecto de sesgo de módulo casi desaparecerá.

bobobobo
fuente
3
El tuyo es más eficiente y probablemente sea cierto que si RAND_MAX es significativamente mayor que el número por el que estás modificando, sin embargo, el tuyo seguirá siendo parcial. De todos modos, todos estos son generadores de números pseudoaleatorios y eso en sí mismo es un tema diferente, pero si asume un generador de números completamente aleatorio, su forma aún sesga los valores más bajos.
user1413793
Como el valor más alto es impar, MT::genrand_int32()%2elige el 0 (50 + 2.3e-8)% del tiempo y el 1 (50 - 2.3e-8)% del tiempo. A menos que esté construyendo el RGN de ​​un casino (para el que probablemente usaría un RGN de ​​rango mucho mayor), cualquier usuario no notará un 2.3e-8% adicional del tiempo. Estás hablando de números demasiado pequeños para importar aquí.
bobobobo
77
Looping es la mejor solución. No es "increíblemente ineficiente"; requiere menos del doble de las iteraciones en el peor de los casos promedio. El uso de un RAND_MAXvalor alto disminuirá el sesgo del módulo, pero no lo eliminará. Looping will.
Jared Nielsen
55
Si RAND_MAXes lo suficientemente grande como el número por el que está modificando, el número de veces que necesita regenerar el número aleatorio es muy pequeño y no afectará la eficiencia. Le digo que siga el ciclo, siempre que esté probando contra el múltiplo más grande en nlugar de solo ncomo lo propone la respuesta aceptada.
Mark Ransom
-3

Acabo de escribir un código para el método imparcial de volteo de monedas de Von Neumann, que teóricamente 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;
        }
    }
}
Yavuz Koroglu
fuente
Esto no aborda el sesgo de módulo. Este proceso podría usarse para eliminar el sesgo en un flujo de bits. Sin embargo, para pasar de un flujo de bits a una distribución uniforme de 0 a n, donde n es uno menos que una potencia de dos, se requiere un sesgo de módulo de direccionamiento. Por lo tanto, esta solución no puede eliminar ningún sesgo en el proceso de generación de números aleatorios.
Rick
2
@Rick hmm. La extensión lógica del método de Von Neumann para eliminar el sesgo de módulo cuando se genera un número aleatorio entre, digamos, 1 y 100, sería: A) llamar rand() % 100100 veces. B) si todos los resultados son diferentes, tome el primero. C) de lo contrario, GOTO A. Esto funcionará, pero con un número esperado de iteraciones de aproximadamente 10 ^ 42, tendrá que ser bastante paciente. E inmortal.
Mark Amery
@ MarkAmery De hecho, eso debería funcionar. Revisando este algoritmo aunque no está implementado correctamente. El primero debería ser:else if(prev==2) prev= x1; else { if(prev!=x1) return prev; prev=2;}
Rick