Simuler une variable de Bernoulli avec

Quelqu'un peut-il me dire comment simuler , où , en utilisant un tirage au sort (autant de fois que vous le souhaitez) avec ?a,b∈NP(H)=$\mathrm{Bernoulli}\left({a\over b}\right)$$a,b\in \mathbb{N}$$P(H)=p$

Je pensais utiliser un échantillonnage de rejet, mais je ne pouvais pas le déterminer.

Parce qu'il existe d'innombrables solutions, trouvons-en une efficace .

L'idée derrière celle-ci commence par une manière standard d'implémenter une variable de Bernoulli: comparer une variable aléatoire uniforme au paramètre a / b . Lorsque U < a / b , retourne 1 ; sinon, retournez 0 .$U$$a/b$$U \lt a/b$$1$$0$

Nous pouvons utiliser le -coin comme générateur de nombres aléatoires uniformes$p$ . Pour générer un nombre uniformément dans n'importe quel intervalle [ x , y ) , lancez la pièce. Quand c'est tête, générez récursivement une valeur uniforme X dans la première partie p de l'intervalle; quand c'est pile, générer récursivement X à partir du dernier 1 - $U$$[x, y)$$X$$p$$X$$1-p$une partie de l'intervalle. À un moment donné, l'intervalle cible deviendra si petit qu'il importe peu de savoir comment vous en choisissez un nombre: c'est ainsi que la récursion commence. Il est évident que cette procédure génère des variations uniformes (jusqu'à n'importe quelle précision souhaitée), comme cela est facilement prouvé par induction.

Cette idée n'est pas efficace, mais elle conduit à une méthode efficace. Étant donné qu'à chaque étape, vous allez dessiner un nombre à partir d'un intervalle donné , pourquoi ne pas d'abord vérifier si vous devez le dessiner? Si votre valeur cible se situe en dehors de cet intervalle, vous connaissez déjà le résultat de la comparaison entre la valeur aléatoire et la cible. Ainsi, cet algorithme a tendance à se terminer rapidement. (Cela pourrait être interprété comme la procédure d' échantillonnage de rejet demandée dans la question.)$[x,y)$

Nous pouvons optimiser davantage cet algorithme. À tout moment, nous avons en fait deux pièces que nous pouvons utiliser: en réétiquetant notre pièce, nous pouvons en faire une pièce avec des chances de chance . Par conséquent, en tant que précalcul, nous pouvons choisir récursivement le réétiquetage qui conduit au plus petit nombre attendu de retournements nécessaires pour la terminaison. (Ce calcul peut être une étape coûteuse.)$1-p$

Par exemple, il est inefficace d'utiliser une pièce avec pour émuler directement une variable de Bernoulli ( 0,01 ) : cela prend près de dix flips en moyenne. Mais si nous utilisons une pièce p = 1 - 0,0 = 0,1 , alors en seulement deux flips, nous serons sûrs d'avoir terminé et le nombre attendu de flips n'est que de 1,2 .$p=0.9$$(0.01)$$p=1-0.0=0.1$$1.2$

Voici les détails.

Partitionner tout intervalle semi-ouvert donné dans les intervalles$I = [x, y)$

Ceci définit les deux transformations et s ( ∗ , T ) qui opèrent à intervalles semi-ouverts.$s(*,H)$$s(*,T)$

En termes de terminologie, si est un ensemble de nombres réels, laissez l'expression$I$

dire que est une limite inférieure pour I : t < x pour tout x ∈ I . De même, t > I signifie t est une limite supérieure pour I .$t$$I$$t \lt x$$x \in I$$t \gt I$$t$$I$

Écrivez . (En fait, cela ne fera aucune différence si t est réel au lieu de rationnel; nous exigeons seulement que 0 ≤ t ≤ 1. )$a/b = t$$t$$0 \le t \le 1$

Voici l'algorithme pour produire une variable avec le paramètre de Bernoulli souhaité:$Z$

1. Définissez et I n = I 0 = [ 0 , 1 ) .$n=0$$I_n = I_0 = [0,1)$

2. Tandis que {Lancez la pièce pour produire X n + 1 . Réglez I n + 1 = S ( I n , X n + 1 ) . Incrémenter n .}$(t\in I_{n})$$X_{n+1}$$I_{n+1} = S(I_n, X_{n+1}).$$n$

3. Si définissez Z = 1 . Sinon, définissez Z = 0 .$t \gt I_{n+1}$$Z=1$$Z=0$

la mise en oeuvre

Pour illustrer, voici une Rimplémentation de l'alorithme comme fonction draw. Ses arguments sont la valeur cible et l'intervalle [ x , y ) , initialement [ 0 , 1 ) . Il utilise la fonction auxiliaire implémentant s . Bien que cela ne soit pas nécessaire, il suit également le nombre de lancers de pièces. Il renvoie la variable aléatoire, le nombre de lancers et le dernier intervalle qu'il a inspecté.$t$$[x,y)$$[0,1)$s$s$

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))
}

A titre d'exemple de l'utilisation et de test de sa précision, prenons le cas et p = 0,9 . Tirons 10 , 000 valeurs en utilisant l'algorithme, rapport sur la moyenne (et son écart - type), et indiquer le nombre moyen de flips utilisé.$t=1/100$$p=0.9$$10,000$

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

Dans cette simulation, des flips étaient des têtes. Bien que inférieur à l'objectif de 0,01 , le score Z de - 0,5154 n'est pas significatif: cet écart peut être attribué au hasard. Le nombre moyen de flips était de 9,886 - un peu moins de dix. Si nous avions utilisé le 1 - p pièce, la moyenne aurait été 0,0094 --still pas significativement différent de la cible, mais seulement 1,177 flips aurait été nécessaire en moyenne.$0.0095$$0.01$$-0.5154$$9.886$$1-p$$0.0094$$1.177$

Voici une solution (un peu désordonnée, mais c'est mon premier coup de couteau). Vous pouvez réellement ignorer le et WLOG supposer P ( H ) = une / deux . Pourquoi? Il existe un algorithme intelligent pour générer un lancer de pièce non biaisé à partir de deux tours de pièce biaisés. Donc , on peut supposer P ( H ) = 1 / 2 .$P(H) = p$$P(H)=1/2$$P(H)=1/2$

Pour générer un , je peux penser à deux solutions (la première n'est pas la mienne, mais la seconde est une généralisation):$\text{Bernoulli}(\frac{a}{b})$

Solution 1

$b$$a$$a$$P(\text{first coin is heads | $a$ heads in $b$ coins}) = \frac{a}{b}$

Solution 2

$\text{Bernoulli}(p)$$p$$0.1 = 0.0001100110011001100110011... \text{base 2}$

$0.$$p$ $bin(p)$

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

Une preuve:

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

est d'environ 0,4 (pas rapide cependant)

Je vois une solution simple, mais il y a sans aucun doute de nombreuses façons de le faire, certaines plus simples sans doute. Cette approche peut être décomposée en deux étapes:

1. $p$$H^*$$T^*$$H^*=(H,T)$$T^*=(T,H)$

2. $a$$-(b-a)$$\frac{a}{a+(b-a)} = \frac{a}{b}$

   (Il y a quelques calculs à faire ici pour le montrer, mais vous pouvez obtenir les probabilités assez facilement en travaillant avec des relations de récurrence ... ou vous pouvez le faire en sommant des séries infinies ... ou il existe d'autres moyens.)