Simuler une variable de Bernoulli avec

9

Quelqu'un peut-il me dire comment simuler , où , en utilisant un tirage au sort (autant de fois que vous le souhaitez) avec ?a,bNP(H)=pBernoulli(ab)a,bNP(H)=p

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

Abracadabra
la source
1
Est-ce une question originaire d'un cours ou d'un manuel? Si oui, veuillez ajouter la [self-study]balise et lire son wiki . Notez qu'il n'est pas nécessaire de demander de l'aide à la fin de votre question - nous savons que tous ceux qui publient ici espèrent de l'aide!
Silverfish
1
Il y a un excellent article de @Glen_b ici quelque part (bien que je ne me souvienne pas où) sur pourquoi il n'y a pas de "pièce biaisée avec probabilité ", mais je sais que ce n'est qu'un problème périphérique à votre question! p
Silverfish
2
@dsaxton La question dit "autant que vous le souhaitez"; il sera fini avec la probabilité 1, mais non limité (vous pourriez dépasser un nombre fixe de lancers) mais objecter sur cette base reviendrait à dire "lancer une pièce juste jusqu'à ce que vous obteniez une tête" n'est pas viable comme méthode pour générer des géométries ( nombres aléatoires.12
Glen_b -Reinstate Monica
1
@AbracaDabra Est-ce un exercice pour une classe? Sinon, comment cela se produit-il?
Glen_b -Reinstate Monica
1
@Glen_b: Ce n'est pas un exercice de ma classe. Cela m'est venu à l'esprit dans cette chaîne de pensée ...: Selon la probabilité classique, prenez une pièce équitable, en augmentant le nombre de lancers, rapport de converge vers la moitié. Donc, cela doit être vrai aussi pour ceux biaisés ... Cela signifie que pour faire converger une pièce vers un nombre particulier, vous avez besoin duP(H)pour être ce nombre. Maintenant, je pensais, que se passe-t-il si nous voulons produire un nombre, mais nous avons une pièce avecP(H)un autre nombre (connu ou inconnu)? #Heads#tailsP(H)P(H)
AbracaDabra

Réponses:

8

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 .Ua/bU<a/b10

Nous pouvons utiliser le -coin comme générateur de nombres aléatoires uniformesp . 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 - pU[x,y)XpX1pune 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.)1p

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=10.0=0.11.2

Voici les détails.

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

[x,y)=[x,x+(yx)p)[x+(yx)p,y)=s(I,H)s(I,T).

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'expressionI

t<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 .tIt<xxIt>ItI

É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=tt0t1

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=0In=I0=[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 .}(tIn)Xn+1In+1=S(In,Xn+1).n

  3. Si définissez Z = 1 . Sinon, définissez Z = 0 .t>In+1Z=1Z=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)ss

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/100p=0.910,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.00950.010.51549.8861p0.00941.177

whuber
la source
Je ne peux pas m'empêcher de voir les similitudes entre cette solution et la solution 2 dans ma réponse. Alors que je suppose une pièce non biaisée (PS solution vraiment intéressante au problème des pièces biaisées) et que je fais tous les calculs / comparaisons en base 2, vous faites tous les calculs / comparaisons en base 10. Que pensez-vous?
Cam.Davidson.Pilon
1
@cam Je pense que vous pourriez être trompé par mes exemples: bien qu'ils utilisent de bons nombres dans la base 10, la construction n'a rien à voir avec une base particulière.
whuber
2
(+1) Résolution très nette. L'optimisation se situe dans les limites supérieures et inférieures par des puissances comme p n ( 1 - p ) m et / ou ( n + ma/bpn(1p)m. Il serait intéressant de trouver la dichotomie optimale en termes de nombre de Bernoullis simulés. (n+mm)pn(1p)m
Xi'an
5

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)=pP(H)=1/2P(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):Bernoulli(ab)

Solution 1

baaP(first coin is heads | a heads in b coins)=ab

Solution 2

Bernoulli(p)p0.1=0.0001100110011001100110011...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)

Cam.Davidson.Pilon
la source
Belle réponse, mais pouvez-vous expliquer avec votre méthode 1 comment faire pour p irrationnel?
AbracaDabra
2
p
p(0,1)(1,0)p(1p)1/2
4

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. pHTH=(H,T)T=(T,H)

  2. a(ba)aa+(ba)=ab

    (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.)

Glen_b -Reinstate Monica
la source