Simular una variable Bernoulli con probabilidad ${a\over b}$ utilizando una moneda sesgada

Question

¿Puede alguien decirme cómo simular $\mathrm{Bernoulli}\left({a\over b}\right)$ donde $a,b\in \mathbb{N}$ mediante el lanzamiento de una moneda (tantas veces como se desee) con $P(H)=p$ ?

Estaba pensando en utilizar el muestreo de rechazo, pero no pude concretarlo.

Answer 1

Porque hay incontables soluciones, encontremos una eficiente uno.

La idea detrás de ésta parte de una forma estándar de implementar una variable Bernoulli: comparar una variable aleatoria uniforme $U$ al parámetro $a/b$ . Cuando $U \lt a/b$ devolución $1$ ; en caso contrario, devuelve $0$ .

Podemos utilizar el $p$ -coin como generador uniforme de números aleatorios . Para generar un número $U$ uniformemente dentro de cualquier intervalo $[x, y)$ Lanza la moneda. Si sale cara, genera recursivamente un valor uniforme $X$ en la primera $p$ parte del intervalo; cuando es cruz, genera recursivamente $X$ del último $1-p$ parte del intervalo. En algún momento, el intervalo objetivo será tan pequeño que no importará cómo se elija un número de él: así es como se inicia la recursión. Es obvio que este procedimiento genera variantes uniformes (hasta cualquier precisión deseada), como se demuestra fácilmente por inducción.

Esta idea no es eficiente, pero conduce a un método eficiente. Dado que en cada etapa vas a extraer un número de un intervalo dado $[x,y)$ ¿Por qué no comprueba primero si necesita dibujarlo? Si su valor objetivo se encuentra fuera de este intervalo, ya conoce el resultado de la comparación entre el valor aleatorio y el objetivo. Por lo tanto, este algoritmo tiende a terminar rápidamente. (Esto podría interpretarse como la muestreo de rechazo procedimiento solicitado en la pregunta).

Podemos optimizar aún más este algoritmo. En cualquier etapa, en realidad tenemos dos monedas que podemos utilizar: reetiquetando nuestra moneda podemos convertirla en una que salga cara con el azar $1-p$ . Por lo tanto, como precomputación podemos elegir recursivamente cualquier reetiquetado que conduzca al menor número esperado de vueltas necesarias para la terminación. (Este cálculo puede ser un paso costoso).

Por ejemplo, es ineficiente utilizar una moneda con $p=0.9$ para emular un Bernoulli $(0.01)$ directamente: se necesitan casi diez vueltas de media. Pero si utilizamos una $p=1-0.0=0.1$ moneda, entonces en sólo dos tiradas seguro que habremos terminado y el número esperado de tiradas es de sólo $1.2$ .

He aquí los detalles.

Partición de cualquier intervalo semiabierto $I = [x, y)$ en los intervalos

$$[x,y) = [x, x + (y-x)p) \cup [x + (y-x)p, y) = s(I,H) \cup s(I,T).$$

Esto define las dos transformaciones $s(*,H)$ y $s(*,T)$ que funcionan a intervalos semiabiertos.

Como cuestión de terminología, si $I$ es cualquier conjunto de números reales dejemos que la expresión

$$t \lt I$$

significa que $t$ es un límite inferior para $I$ : $t \lt x$ para todos $x \in I$ . Del mismo modo, $t \gt I$ significa $t$ es un límite superior para $I$ .

Escriba a $a/b = t$ . (De hecho, no habrá diferencia si $t$ sea real en lugar de racional; sólo exigimos que $0 \le t \le 1$ .)

Este es el algoritmo para producir una variante $Z$ con el parámetro Bernoulli deseado:

1. Establecer $n=0$ y $I_n = I_0 = [0,1)$ .

2. En $(t\in I_{n})$ {Lanza la moneda para producir $X_{n+1}$ . Establecer $I_{n+1} = S(I_n, X_{n+1}).$ Incremento $n$ .}

3. Si $t \gt I_{n+1}$ a continuación, establezca $Z=1$ . En caso contrario, fije $Z=0$ .

---

Aplicación

Para ilustrarlo, he aquí un R implementación del aloritmo como la función draw . Sus argumentos son el valor objetivo $t$ y el intervalo $[x,y)$ inicialmente $[0,1)$ . Utiliza la función auxiliar s aplicación de $s$ . Aunque no lo necesita, también hace un seguimiento del número de lanzamientos de la moneda. Devuelve la variable aleatoria, el recuento de lanzamientos y el último intervalo que inspeccionó.

s <- function(x, ab, p) {
  d <- diff(ab) * p
  if (x == 1) c(ab[1], ab[1] + d) else c(ab[1] + d, ab[2])
}
draw <- function(target, p) {
  between <- function(z, ab) prod(z - ab) <= 0
  ab <- c(0,1)
  n <- 0
  while(between(target, ab)) {
    n <- n+1; ab <- s(runif(1) < p, ab, p)
  }
  return(c(target > ab[2], n, ab))
}

Como ejemplo de su uso y prueba de su precisión, tomemos el caso $t=1/100$ y $p=0.9$ . Dibujemos $10,000$ utilizando el algoritmo, informar sobre la media (y su error estándar) e indicar el número medio de volteos utilizados.

target <- 0.01
p <- 0.9
set.seed(17)
sim <- replicate(1e4, draw(target, p))

(m <- mean(sim[1, ]))                           # The mean
(m - target) / (sd(sim[1, ]) / sqrt(ncol(sim))) # A Z-score to compare to `target`
mean(sim[2, ])                                  # Average number of flips

En esta simulación $0.0095$ de las tiradas salieron cara. Aunque inferior al objetivo de $0.01$ la puntuación Z de $-0.5154$ no es significativa: esta desviación puede atribuirse al azar. El número medio de lanzamientos fue de $9.886$ --un poco menos de diez. Si hubiéramos utilizado el $1-p$ moneda, la media habría sido $0.0094$ -aún no significativamente diferente del objetivo, pero sólo $1.177$ por término medio.

Answer 2

He aquí una solución (un poco desordenada, pero es mi primer intento). En realidad se puede ignorar el $P(H) = p$ y WLOG asumen $P(H)=1/2$ . ¿Por qué? Existe un algoritmo inteligente para generar una moneda insesgada a partir de dos monedas sesgadas. Así que podemos suponer $P(H)=1/2$ .

Para generar un $\text{Bernoulli}(\frac{a}{b})$ Se me ocurren dos soluciones (la primera no es mía, pero la segunda es una generalización):

Solución 1

Lanza la moneda sin prejuicios $b$ veces. Si $a$ cabezas no están presentes, empezar de nuevo. Si $a$ cabezas son presente, devuelve si la primera moneda sale cara o no (porque $P(\text{first coin is heads | $ a $ heads in $ b $ coins}) = \frac{a}{b}$ )

Solución 2

Esto puede ampliarse a cualquier valor de $\text{Bernoulli}(p)$ . Escriba a $p$ en forma binaria. Por ejemplo, $0.1 = 0.0001100110011001100110011... \text{base 2}$

Crearemos un nuevo número binario lanzando una moneda al aire. Empezaremos con $0.$ y sumar dígitos en función de si aparece cara (1) o cruz (0). En cada tirada, compara el nuevo número binario con la representación binaria de $p$ hasta el mismo dígito . Eventualmente los dos divergirán, y volverán si $bin(p)$ es mayor que su número binario.

En Python:

def simulate(p):
    binary_p = float_to_binary(p)
    binary_string = '0.'
    index = 3
    while True:
        binary_string += '0' if random.random() < 0.5 else '1'
        if binary_string != binary_p[:index]:
            return binary_string < binary_p[:index]
        index += 1

Alguna prueba:

np.mean([simulate(0.4) for i in range(10000)])

es de aproximadamente 0,4 (aunque no es rápido)

Answer 3

Veo una solución sencilla, pero sin duda hay muchas formas de hacerlo, algunas presumiblemente más sencillas que ésta. Este enfoque se puede dividir en dos pasos:

1. Generar a partir de dos sucesos con igual probabilidad dado un procedimiento injusto de lanzamiento de la moneda (la combinación de la moneda concreta y el método por el que se lanza generando una cara con probabilidad $p$ ). Podemos llamar a estos dos sucesos igualmente probables $H^*$ y $T^*$ . [Hay un enfoque simple para esto que requiere tomar pares de lanzamientos $H^*=(H,T)$ y $T^*=(T,H)$ para producir dos resultados igualmente probables, y todos los demás resultados conducen a generar un nuevo par de tiradas para intentarlo de nuevo].

2. Ahora se genera un paseo aleatorio con dos estados absorbentes utilizando la moneda justa simulada. Eligiendo la distancia de los estados absorbentes desde el origen (uno por encima y otro por debajo), puedes establecer la probabilidad de absorción por, digamos, el estado absorbente superior para que sea una proporción deseada de enteros. En concreto, si se coloca la barrera de absorción superior en $a$ y el inferior en $-(b-a)$ (e iniciamos el proceso desde el origen), y recorremos el camino aleatorio hasta la absorción, la probabilidad de absorción en la barrera superior es $\frac{a}{a+(b-a)} = \frac{a}{b}$ .

   (Aquí hay que hacer algunos cálculos para demostrarlo, pero se pueden sacar las probabilidades con bastante facilidad trabajando con relaciones de recurrencia... o se puede hacer sumando series infinitas... o hay otras formas).