Generando números aleatorios distribuidos uniformemente usando una moneda

Tienes una moneda Puedes voltearlo tantas veces como quieras.

Desea generar un número aleatorio tal que donde .$r$$a \leq r < b$$r,a,b\in \mathbb{Z}^+$

La distribución de los números debe ser uniforme.

Es fácil si :$b -a = 2^n$

r = a + binary2dec(flip n times write 0 for heads and 1 for tails) 

¿Qué pasa si ?$b-a \neq 2^n$

Lo que está buscando se basa en el muestreo de rechazo o el método de aceptación-rechazo (tenga en cuenta que la página Wiki es un poco técnica).

Este método es útil en este tipo de situaciones: desea elegir algún objeto aleatorio de un conjunto (un entero aleatorio en el conjunto en su caso), pero no sabe cómo hacerlo, pero puede elegir algún objeto aleatorio de un conjunto más grande que contenga el primer conjunto (en su caso, para algunos modo que ; esto corresponde a lanzamientos de monedas).$[a,b]$$[a, 2^k + a]$$k$$2^k + a \ge b$$k$

En tal escenario, simplemente sigue eligiendo elementos del conjunto más grande hasta que elige aleatoriamente un elemento en el conjunto más pequeño. Si su conjunto más pequeño es lo suficientemente grande en comparación con su conjunto más grande (en su caso, contiene como máximo el doble de enteros que , que es lo suficientemente bueno), esto es eficiente.$[a, 2^k + a]$$[a,b]$

Un ejemplo alternativo: suponga que desea elegir un punto aleatorio dentro de un círculo con radio 1. Ahora, no es realmente fácil encontrar un método directo para esto. Pasamos al método aceptar-rechazar: tomamos muestras de puntos en un cuadrado de 1x1 que abarca el círculo, y probamos si el número que dibujamos se encuentra dentro del círculo.

(tecnicismos: la respuesta se ajusta a la selección del número )$a \le x < b$

Como se le permite lanzar su moneda tantas veces como lo desee, puede obtener su probabilidad tan cerca como desee de uniformarse eligiendo una fracción (usando la raíz binaria: usted lanza el monedas para cada dígito después del punto) y multiplique r por b - a para obtener un número entre 0 y [ba-1] (redondeando a un entero). Añadir a este número para una y ya está.$r\in [0,1]$$r$$b-a$$a$

Ejemplo : diga . 1/3 en binario es 0.0101010101 .... Entonces, si su flip está entre 0 y 0.010101 ... su elección es b . si está entre 0.010101 .. y 0.10101010 ... su selección será un + 1 , y de lo contrario será un + 2 .$b-a=3$$b$$a+1$$a+2$

Si lanza su moneda veces, cada número entre a y b se elegirá con probabilidad $t$$a$$b$.$\frac{1}{b-a}\pm 2^{-(t+1)}$

Elija un número en la próxima potencia más grande del rango 2 y descarte respuestas mayores que .$b-a$

n = b-a;
N = round_to_next_larger_power_of_2(n)
while (1) {
  x = random(0 included to N excluded);
  if (x < n) break;
}
r = a + x;

Nadie mencionó esto, así que vamos a demostrar formalmente que si es un dominio cuyo tamaño no es una potencia de dos, a continuación, un número finito de lanzamientos de moneda razonable no son suficientes para generar un miembro uniformemente aleatoria de D . Supongamos que utiliza k monedas para generar un miembro de D . Para cada d ∈ D , la probabilidad de que su algoritmo generado d es de la forma A / 2 k para algún entero A . El teorema fundamental de la aritmética muestra que A / 2 k ≠ 1 / | D | .$D$$D$$k$$D$$d \in D$$d$$A/2^k$$A$$A/2^k \neq 1/{|D|}$

Si desea generar muestras uniformes independientes de D , entonces el número esperado de lanzamientos de monedas que necesita (bajo el algoritmo óptimo) es aproximadamente n log 2 | D | . En términos más generales, si desea generar a partir de una distribución de entropía H , necesita aproximadamente n H bits aleatorios en expectativa. Un algoritmo que logra esto es la decodificación aritmética, aplicada a una secuencia (aparentemente) infinita de bits aleatorios.$n$$D$$n \log_2 |D|$$H$$nH$

Si ba no es una potencia de 2, entonces puede que tenga que lanzar muchas monedas para obtener un resultado. Incluso puede que nunca obtenga un resultado, pero eso es poco probable en extremo.

Métodos

El método más simple es generar un número en [a, a + 2 ^ n), donde 2 ^ n> = ba, hasta que uno aterrice en [a, b). Desecha mucha entropía con este método.

Un método más caro le permite mantener toda la entropía, pero se vuelve muy costoso computacionalmente a medida que aumenta el número de lanzamientos de monedas / tiradas de dados. Intuitivamente es como tratar los lanzamientos de monedas como los dígitos de un número binario a la derecha del punto decimal, convirtiendo ese número de la base 2 a la base ab después, y devolviendo los dígitos de ese número a medida que se quedan 'atascados'.

Ejemplo

El siguiente código convierte las tiradas de un dado regular de n lados en tiradas de un dado justo de m lados (en su caso n = 2, m = ab), con un costo marginal creciente a medida que aumenta el número de tiradas. Tenga en cuenta la necesidad de un tipo de número racional con precisión arbitraria. Una buena propiedad es que la conversión de n-cara a m-cara y de vuelta a n-cara devolverá la secuencia original, aunque tal vez se retrase un par de vueltas debido a que los dígitos tienen que atascarse.

public static IEnumerable<BigInteger> DigitConversion(this IEnumerable<BigInteger> inputStream, BigInteger modIn, BigInteger modOut) {
    //note: values are implicitly scaled so the first unfixed digit of the output ranges from 0 to 1
    Rational b = 0; //offset of the chosen range
    Rational d = 1; //size of the chosen range
    foreach (var r in inputStream) {
        //narrow the chosen range towards the real value represented by the input
        d /= modIn;
        b += d * r;
        //check for output digits that have become fixed
        while (true) {
            var i1 = (b * modOut).Floor();
            var i2 = ((b + d) * modOut).Floor(); //note: ideally b+d-epsilon, but another iteration makes that correction unnecessary
            if (i1 != i2) break; //digit became fixed?
            //fix the next output digit (rescale the range to make next digit range from 0 to 1)
            d *= modOut;
            b *= modOut;
            b -= i1;
            yield return i1;
        }
    }
}

Genera un decimal binario. En lugar de almacenarlo explícitamente, solo realice un seguimiento de los valores mínimos y máximos posibles. Una vez que esos valores se encuentran dentro del mismo entero, devuelve ese entero. El bosquejo del código está abajo.

(Editar) Explicación más completa: supongamos que desea generar un número entero aleatorio de 1 a 3 inclusive con 1/3 de probabilidad cada uno. Hacemos esto generando un decimal binario aleatorio real, x en el rango (0, 1). Si x <1/3, devuelve 1, de lo contrario si x <2/3 devuelve 2, de lo contrario devuelve 3. En lugar de generar los dígitos para x explícitamente, solo hacemos un seguimiento de los valores mínimos y máximos posibles de x. Inicialmente, el valor mínimo de x es 0 y el máximo es 1. Si primero gira las cabezas, entonces el primer dígito de x detrás del punto decimal (en binario) es 1. El valor mínimo posible de x (en binario) se convierte en 0.100000 = 1/2 y el máximo es 0.111111111 = 1. Ahora, si su próximo giro es colas, x comienza con 0.10. El valor mínimo posible es 0.1000000 = 1/2 y el máximo es 0.1011111 = 3/4. El valor mínimo posible de x es 1/2 para que sepa que existe ' s no hay posibilidad de devolver 1 ya que eso requiere x <1/3. Aún puede devolver 2 si x termina siendo 1/2 <x <2/3 o 3 si 2/3 <x <3/4. Ahora supongamos que la tercera vuelta es la cola. Entonces x debe comenzar con 0.100. Min = 0.10000000 = 1/2 y max = 0.100111111 = 5/8. Ahora, desde 1/3 <1/2 <5/8 <2/3 sabemos que x debe caer en el intervalo (1/3, 2/3), por lo que podemos dejar de generar dígitos de x y simplemente devolver 2.

El código hace esencialmente esto excepto que en lugar de generar x entre 0 y 1 genera x entre a y b, pero el principio es el mismo.

def gen(a, b):
  min_possible = a
  max_possible = b

  while True:
    floor_min_possible = floor(min_possible)
    floor_max_possible = floor(max_possible)
    if max_possible.is_integer():
      floor_max_possible -= 1

    if floor_max_possible == floor_min_possible:
      return floor_max_possible

    mid = (min_possible + max_possible)/2
    if coin_flip():
      min_possible = mid
    else:
      max_possible = mid

Observación: probé este código con el método de aceptar / rechazar y ambos producen distribuciones uniformes. Este código requiere menos lanzamientos de monedas que aceptar rechazar, excepto cuando b - a está cerca de la siguiente potencia de 2. Ej. Si desea generar a = 0, b = 62, entonces aceptar / rechazar funciona mejor. Pude demostrar que este código puede usar en promedio no más de 2 lanzamientos de monedas que aceptar / rechazar. Según mi lectura, parece que Knuth y Yao (1976) dieron un método para resolver este problema y demostraron que su método es óptimo en el número esperado de lanzamientos de monedas. Además demostraron que el número esperado de cambios debe ser mayor que la entropía de Shannon de la distribución. Sin embargo, no pude encontrar una copia del texto del documento y sería curioso ver cuál es su método. (Actualización: acabo de encontrar una exposición de Knuth Yao 1976 aquí:http://www.nrbook.com/devroye/Devroye_files/chapter_fifteen_1.pdf pero aún no lo he leído). Alguien también mencionó a Han Hoshi en este hilo que parece ser más general y lo resuelve usando una moneda sesgada. Ver también http://paper.ijcsns.org/07_book/200909/20090930.pdf por Pae (2009) para una buena discusión de la literatura.

La respuesta simple?

$b-a$$r$

$r$$b - a = 2^n + 1$$n$

Esta es una solución propuesta para el caso cuando b - a no es igual a 2 ^ k. Se supone que funciona en un número fijo de pasos (no es necesario descartar candidatos que están fuera de su rango esperado).

Sin embargo, no estoy seguro de que esto sea correcto. Critique y ayude a describir la no uniformidad exacta en este generador de números aleatorios (si existe) y cómo medirlo / cuantificarlo.

Primero convierta al problema equivalente de generar números aleatorios distribuidos uniformemente en el rango [0, z-1] donde z = b - a.

Además, sea m = 2 ^ k la potencia más pequeña de 2> = z.

Según la solución anterior, ya tenemos un generador de números aleatorios R (m) distribuido uniformemente en el rango [0, m-1] (se puede hacer arrojando k monedas, una por cada bit).

    Keep a random seed s and initialize with s = R(m).   

    function random [0, z-1] :
        x = R(m) + s 
        while x >= z:
            x -= z
        s = x
        return x

El ciclo while se ejecuta como máximo 3 veces, dando el siguiente número aleatorio en un número fijo de pasos (mejor caso = peor caso).

Vea un programa de prueba para los números [0,2] aquí: http://pastebin.com/zuDD2V6H

Algoritmo teóricamente óptimo

Aquí hay una mejora de la otra respuesta que publiqué. La otra respuesta tiene la ventaja de que es más fácil extenderse al caso más general de generar una distribución discreta a partir de otra. De hecho, la otra respuesta es un caso especial del algoritmo debido a Han y Hoshi.

El algoritmo que describiré aquí se basa en Knuth y Yao (1976). En su artículo, también demostraron que este algoritmo logra el mínimo número esperado posible de lanzamientos de monedas.

Para ilustrarlo, considere el método de muestreo de rechazo descrito por otras respuestas. Como ejemplo, suponga que desea generar uno de 5 números de manera uniforme [0, 4]. La siguiente potencia de 2 es 8, por lo que lanza la moneda 3 veces y genera un número aleatorio hasta 8. Si el número es 0 a 4, entonces lo devuelve. De lo contrario, lo tira y genera otro número hasta el 8 e intenta nuevamente hasta que tenga éxito. Pero cuando arrojas el número, desperdicias algo de entropía. En su lugar, puede condicionar el número que arrojó para reducir el número de futuros lanzamientos de monedas que necesitará con expectativa. Concretamente, una vez que genera el número [0, 7], si es [0, 4], regrese. De lo contrario, son 5, 6 o 7, y haces algo diferente en cada caso. Si es 5, voltea la moneda nuevamente y devuelve 0 o 1 según el lanzamiento. Si son las 6, voltee la moneda y devuelva 2 o 3. Si es 7, voltee la moneda; si es cara, devuelve 4, si es cruz comienza de nuevo.

La entropía sobrante de nuestro intento fallido inicial nos dio 3 casos (5, 6 o 7). Si solo tiramos esto, tiramos los lanzamientos de monedas log2 (3). En cambio, lo guardamos y lo combinamos con el resultado de otro cambio para generar 6 casos posibles (5H, 5T, 6H, 6T, 7H, 7T), lo que nos permite volver a intentarlo de inmediato para generar una respuesta final con probabilidad de éxito 5/6 .

Aquí está el código:

# returns an int from [0, b)
def __gen(b):
  rand_num = 0
  num_choices = 1

  while True:
    num_choices *= 2
    rand_num *= 2
    if coin.flip():
      rand_num += 1

    if num_choices >= b:
      if rand_num < b:
        return rand_num
      num_choices -= b
      rand_num -= b

# returns an int from [a, b)
def gen(a, b):
  return a + __gen(b - a)