Avec quelle probabilité une pièce est-elle meilleure que l'autre?

8

Disons que nous avons deux pièces biaisées C1et que les C2deux ont une probabilité différente de tourner la tête.

Nous jetons des C1 n1fois et obtenons des H1têtes, des C2 n2fois et obtenons des H2têtes. Et nous constatons que le rapport des têtes pour une pièce est plus élevé que l'autre.

Quelle est la probabilité avec laquelle nous pouvons dire qu'une pièce est meilleure que l'autre? (mieux signifie ici une probabilité réelle plus élevée de tourner la tête).

Thirupathi Thangavel
la source
2
N'est-ce pas un problème de test d'hypothèse où vous devez tester si les estimations de la probabilité de relever les têtes pour les pièces sont différentes.
naïf
3
Si vous voulez une probabilité qu'une pièce soit meilleure qu'une autre, vous pouvez le faire avec une approche bayésienne; tout ce dont vous avez besoin est une distribution préalable de la probabilité. C'est un calcul assez simple.
Glen_b -Reinstate Monica
2
Il n'y a aucun moyen par lequel ces données vous donneront une probabilité qu'une pièce soit meilleure que l'autre . Vous avez besoin d'informations / d'hypothèses supplémentaires sur une distribution antérieure de la façon dont les pièces sont distribuées. Disons que dans une approche bayésienne, le résultat différera beaucoup en fonction de vos hypothèses antérieures. Par exemple, si ces pièces sont des pièces ordinaires, elles sont (a priori) très susceptibles d'être égales et vous auriez besoin d'un grand nombre d'essais pour trouver une probabilité décente que l'une est `` meilleure '' que l'autre (car la probabilité est élevée que vous a obtenu deux pièces similaires qui testent accidentellement différentes)
Sextus Empiricus
1
@katosh, c'est aussi exact que le précédent. Supposons que vous ayez des pièces entre p = 0,49 et p = 0,51. Si dans ce cas vous utilisez l'hypothèse que les pièces sont uniformes réparties entre 0 et 1, alors les probabilités que vous attribuez sont plus comme de mauvaises hypothèses mises à jour plutôt que des probabilités mises à jour. Vous assignerez trop souvent une forte probabilité d'être «meilleur» à la mauvaise pièce. Il serait faux de considérer le résultat comme «la» probabilité que la pièce soit «meilleure» que l'autre. Données erronées = mauvais résultats. Dans ce problème, nous ne devrions pas travailler à résoudre les mathématiques mais à résoudre les connaissances.
Sextus Empiricus
1
De cette façon, il est indiqué correctement. C'est un «croire» et ce n'est pas si facilement assimilé à une «probabilité».
Sextus Empiricus

Réponses:

7

Il est facile de calculer la probabilité de faire cette observation, étant donné que les deux pièces sont égales. Cela peut être fait par le test exact des pêcheurs . Compte tenu de ces observations

coin 1coin 2headsH1H2tailsn1H1n2H2

la probabilité d'observer ces nombres alors que les pièces sont égales étant donné le nombre d'essais n1, n2 et le nombre total de têtes H1+H2 est

p(H1,H2|n1,n2,H1+H2)=(H1+H2)!(n1+n2H1H2)!n1!n2!H1!H2!(n1H1)!(n2H2)!(n1+n2)!.

Mais ce que vous demandez, c'est la probabilité qu'une pièce soit meilleure. Puisque nous discutons de la croyance sur le biais des pièces, nous devons utiliser une approche bayésienne pour calculer le résultat. Veuillez noter que dans l'inférence bayésienne, le terme croyance est modélisé comme probabilité et que les deux termes sont utilisés de manière interchangeable (s. Probabilité bayésienne ). Nous appelons la probabilité que la piècei lance des têtes pi. La distribution postérieure après observation, pour cettepiest donné par le théorème de Bayes :

f(pi|Hi,ni)=f(Hi|pi,ni)f(pi)f(ni,Hi)
La fonction de densité de probabilité (pdf) f(Hi|pi,ni)est donnée par la probabilité binomiale, car les essais individuels sont des expériences de Bernoulli: je suppose la connaissance antérieure sur est que pourrait se situer n'importe où entre et avec une probabilité égale, d'où . Le nominateur est donc .
f(Hi|pi,ni)=(niHi)piHi(1pi)niHi
f(pi)pi01f(pi)=1f(Hi|pi,ni)f(pi)=f(Hi|pi,ni)

Pour calculer nous utilisons le fait que l'intégrale sur un pdf doit être un . Le dénominateur sera donc un facteur constant pour y parvenir. Il existe un pdf connu qui ne diffère du nominateur que par un facteur constant, qui est la distribution bêta . D'où f(ni,Hi)01f(p|Hi,ni)dp=1

f(pi|Hi,ni)=1B(Hi+1,niHi+1)piHi(1pi)niHi.

Le pdf de la paire de probabilités de pièces indépendantes est

f(p1,p2|H1,n1,H2,n2)=f(p1|H1,n1)f(p2|H2,n2).

Maintenant, nous devons l'intégrer dans les cas où afin de savoir comment la pièce est probablement meilleure que la pièce : p1>p212

P(p1>p2)=010p1f(p1,p2|H1,n1,H2,n2)dp2dp1=01B(p1;H2+1,n2H2+1)B(H2+1,n2H2+1)f(p1|H1,n1)dp1

Je ne peux pas résoudre cette dernière intégrale analytiquement mais on peut la résoudre numériquement avec un ordinateur après avoir branché les chiffres. est la fonction bêta et est la fonction bêta incomplète. Notez que car est une variable continue et jamais exactement la même chose que .B(,)B(;,)P(p1=p2)=0p1p2


Concernant l'hypothèse antérieure sur et ses remarques: Une bonne alternative au modèle est considérée par beaucoup comme une distribution bêta . Cela conduirait à une probabilité finale De cette façon, on pourrait modéliser un fort biais vers les pièces de monnaie régulières de grand mais égal , . Cela reviendrait à lancer la pièce des fois supplémentaires et à recevoir des têtes ce qui équivaudrait à avoir simplement plus de données. est la quantité de lancers que nous n'aurions pas à effectuerf(pi)Beta(ai+1,bi+1)

P(p1>p2)=01B(p1;H2+1+a2,n2H2+1+b2)B(H2+1+a2,n2H2+1+b2)f(p1|H1+a1,n1+a1+b1)dp1.
aibiai+biaiai+bi si nous incluons cela avant.

Le PO a déclaré que les deux pièces étaient toutes deux biaisées à un degré inconnu. J'ai donc compris que toutes les connaissances doivent être déduites des observations. C'est pourquoi j'ai opté pour une méthode non informative avant que la dose ne biaise le résultat, par exemple vers des pièces ordinaires.

Toutes les informations peuvent être transmises sous la forme de par pièce. L'absence d'un a priori informatif signifie seulement que plus d'observations sont nécessaires pour décider quelle pièce est meilleure avec une probabilité élevée.(Hi,ni)


Voici le code de R qui fournit une fonction utilisant l'uniforme a priori : P(n1, H1, n2, H2) =P(p1>p2)f(pi)=1

mp <- function(p1, n1, H1, n2, H2) {
    f1 <- pbeta(p1, H2 + 1, n2 - H2 + 1)
    f2 <- dbeta(p1, H1 + 1, n1 - H1 + 1)
    return(f1 * f2)
}

P <- function(n1, H1, n2, H2) {
    return(integrate(mp, 0, 1, n1, H1, n2, H2))
}

Vous pouvez dessiner pour différents résultats expérimentaux et fixer , par exemple avec ce code coupé:P(p1>p2)n1n2n1=n2=4

library(lattice)
n1 <- 4
n2 <- 4
heads <- expand.grid(H1 = 0:n1, H2 = 0:n2)
heads$P <- apply(heads, 1, function(H) P(n1, H[1], n2, H[2])$value)
levelplot(P ~ H1 + H2, heads, main = "P(p1 > p2)")

P (p1> p2) pour n1 = n2 = 4

Vous devrez peut-être d' install.packages("lattice")abord.

On peut voir que même avec un a priori uniforme et un petit échantillon, la probabilité ou croire qu'une pièce est meilleure peut devenir assez solide, lorsque et diffèrent suffisamment. Une différence relative encore plus petite est nécessaire si et sont encore plus grands. Voici un graphique pour et :H1H2n1n2n1=100n2=200

entrez la description de l'image ici


Martijn Weterings a suggéré de calculer la distribution de probabilité postérieure de la différence entre et . Cela peut être fait en intégrant le pdf de la paire sur l'ensemble : p1p2S(d)={(p1,p2)[0,1]2|d=|p1p2|}

f(d|H1,n1,H2,n2)=S(d)f(p1,p2|H1,n1,H2,n2)dγ=01df(p,p+d|H1,n1,H2,n2)dp+d1f(p,pd|H1,n1,H2,n2)dp

Encore une fois, pas une intégrale que je peux résoudre analytiquement mais le code R serait:

d1 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p + d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

d2 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p - d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

fd <- function(d, n1, H1, n2, H2) {
    if (d==1) return(0)
    s1 <- integrate(d1, 0, 1-d, d, n1, H1, n2, H2)
    s2 <- integrate(d2, d, 1, d, n1, H1, n2, H2)
    return(s1$value + s2$value)
}

J'ai tracé pour , , et toutes les valeurs de :f(d|n1,H1,n2,H2)n1=4H1=3n2=4H2

n1 <- 4
n2 <- 4
H1 <- 3
d <- seq(0, 1, length = 500)

get_f <- function(H2) sapply(d, fd, n1, H1, n2, H2)
dat <- sapply(0:n2, get_f)

matplot(d, dat, type = "l", ylab = "Density",
        main = "f(d | 4, 3, 4, H2)")
legend("topright", legend = paste("H2 =", 0:n2),
       col = 1:(n2 + 1), pch = "-")

entrez la description de l'image ici

Vous pouvez calculer la probabilité deêtre supérieur à une valeur par . Gardez à l'esprit que la double application de l'intégrale numérique s'accompagne d'une erreur numérique. Par exemple, doit toujours être égal à car prend toujours une valeur comprise entre et . Mais le résultat s'écarte souvent légèrement.|p1p2|dintegrate(fd, d, 1, n1, H1, n2, H2)integrate(fd, 0, 1, n1, H1, n2, H2)1d01

Katosh
la source
Je ne connais pas la valeur réelle de p1
Thirupathi Thangavel
1
Je suis désolé pour ma mauvaise notation 😅 mais vous n'avez pas besoin de brancher une valeur fixe pour . Le (maintenant modifié) dans l'intégrale est la variable sur laquelle vous . Tout comme vous pouvez intégrer sans avoir une valeur fixe pour . p1p101x2dxx
katosh
1
Le test exact de Fisher concerne plus spécifiquement l'hypothèse selon laquelle les pièces ont la même probabilité et les totaux marginaux sont fixes . Ce n'est pas le cas dans ce problème de pièces. Si vous refaites le test, vous pourrez observer une autre quantité totale de têtes.
Sextus Empiricus
@MartijnWeterings dans mon cas, la probabilité de tourner la tête pour une pièce est toujours fixe. N'est-ce pas suffisant?
Thirupathi Thangavel,
1
@ThirupathiThangavel le problème avec le test de Fisher concerne les totaux marginaux non fixes. Le modèle de test exact suppose que la probabilité de têtes sont les mêmes et fixes, mais aussi marginaux sont fixés avant l'expérience. Cette deuxième partie n'est pas le cas. Cela donne une probabilité conditionnelle différente pour les valeurs extrêmes. Dans l'ensemble, le test de Fisher sera conservateur. La probabilité du résultat compte tenu de l'hypothèse VRAIE (c.-à-d. Probabilité fixe et similaire pour les têtes, mais pas nécessairement des totaux marginaux fixes) est inférieure à celle calculée (vous obtenez des valeurs p plus élevées).
Sextus Empiricus
2

J'ai fait une simulation numérique avec R, probablement vous cherchez une réponse analytique, mais j'ai pensé que cela pourrait être intéressant à partager.

set.seed(123)
# coin 1
N1 = 20
theta1 = 0.7

toss1 <- rbinom(n = N1, size = 1, prob = theta1)

# coin 2
N2 = 25
theta2 = 0.5

toss2 <- rbinom(n = N2, size = 1, prob = theta2)

# frequency
sum(toss1)/N1 # [1] 0.65
sum(toss2)/N2 # [1] 0.52

Dans ce premier code, je simule simplement deux lancers de pièces. Ici, vous pouvez voir bien sûr que cela theta1 > theta2, alors bien sûr la fréquence de H1sera supérieure à H2. Notez les différentes N1, N2tailles.

Voyons ce que nous pouvons faire avec différents thetas. Notez que le code n'est pas optimal. Du tout.

simulation <- function(N1, N2, theta1, theta2, nsim = 100) {
  count1 <- count2 <- 0

  for (i in 1:nsim) {
    toss1 <- rbinom(n = N1, size = 1, prob = theta1)
    toss2 <- rbinom(n = N2, size = 1, prob = theta2)

    if (sum(toss1)/N1 > sum(toss2)/N2) {count1 = count1 + 1} 
    #if (sum(toss1)/N1 < sum(toss2)/N2) {count2 = count2 + 1} 
  }

  count1/nsim

}
set.seed(123)
simulation(20, 25, 0.7, 0.5, 100)
#[1] 0.93

Donc 0,93 est la fréquence (sur 100) où la première pièce avait plus de têtes. Cela semble correct, vu theta1et theta2utilisé.

Voyons voir avec deux vecteurs de thetas.

theta1_v <- seq(from = 0.1, to = 0.9, by = 0.1)
theta2_v <- seq(from = 0.9, to = 0.1, by = -0.1)

res_v <- c()
for (i in 1:length(theta1_v)) {

  res <- simulation(1000, 1500, theta1_v[i], theta2_v[i], 100)
  res_v[i] <- res

}

plot(theta1_v, res_v, type = "l")

N'oubliez pas que ce res_vsont les fréquences où H1 > H2, sur 100 simulations.

entrez la description de l'image ici

Donc, à theta1mesure que les augmentations, la probabilité d' H1être plus élevé augmente, bien sûr.

Je l' ai fait d'autres simulations et il semble que les tailles N1, N2sont moins importants.

Si vous connaissez, Rvous pouvez utiliser ce code pour faire la lumière sur le problème. Je suis conscient que ce n'est pas une analyse complète et qu'elle peut être améliorée.

RLave
la source
2
Intéressant de voir comment res_vchange continuellement lorsque les thétas se rencontrent. J'ai compris la question car elle posait des questions sur le biais intrinsèque des pièces après avoir fait une seule observation. Vous semblez répondre aux observations que l'on ferait après avoir connu le biais.
katosh