Intervalle de prédiction pour la variable aléatoire binomiale

14

Quelle est la formule (approximative ou exacte) d'un intervalle de prédiction pour une variable aléatoire binomiale?

Supposons que , et nous observons (tiré de ). Le est connu.YBinom(n,p)yYn

Notre objectif est d'obtenir un intervalle de prédiction à 95% pour un nouveau tirage de .Y

L'estimation ponctuelle est , où \ chapeau {p} = \ frac {y} {n} . Un intervalle de confiance pour \ hat {p} est simple, mais je ne peux pas trouver une formule pour un intervalle de prédiction pour Y . Si nous connaissions p (plutôt que \ hat {p} ), alors un intervalle de prédiction à 95% implique simplement de trouver les quantiles d'un binôme. Y a-t-il quelque chose d'évident que j'oublie?p = ynp^p Yp pp^=ynp^Ypp^

Statseeker
la source
1
Voir Quelles méthodes non bayésiennes existe-t-il pour l'inférence prédictive? . Dans ce cas, la méthode utilisant des pivots n'est pas disponible (je ne pense pas) mais vous pouvez utiliser l'une des probabilités prédictives. Ou bien sûr, une approche bayésienne.
Scortchi - Réintégrer Monica
1
Salut les gars, je voudrais prendre un moment pour répondre aux préoccupations qui ont été soulevées. - concernant la confiance pour p: je ne suis pas intéressé pour ça. - en ce qui concerne les prédictions à 95% de la distribution: oui, c'est exactement ce que sont les intervalles de prédiction quel que soit le contexte (dans la régression, vous devez supposer des erreurs normales, alors que les intervalles de confiance dépendent du CLT - oui, l'exemple de prédiction du nombre de têtes dans un tirage au sort est correct. Ce qui rend ce problème difficile, c'est que nous n'avons pas maintenant "p", jut ont une estimation.
Statseeker
3
@Addison Lisez le livre Statistical Intervals de G. Hahn et W. Meeker. Ils expliquent la différence entre les intervalles de confiance, les intervalles de prédiction, les intervalles de tolérance et les intervalles crédibles bayésiens. Un intervalle de prédiction à 95% ne contient pas 95% de la distribution. Il fait ce que font les intervalles les plus fréquentistes. Si vous échantillonnez à plusieurs reprises à partir du B (n, p) et utilisez la même méthode à chaque fois pour produire un intervalle de prédiction de 95% pour p, puis 95% des intervalles de prédiction, vous contiendrez la vraie valeur de p. Si vous voulez couvrir 95% de la distribution, construisez un intervalle de tolérance.
Michael R. Chernick
Les intervalles de tolérance couvrent un pourcentage de la distribution. Pour un intervalle de tolérance de 95% pour 90% de la distribution, vous répétez à nouveau le processus plusieurs fois et utilisez la même méthode pour générer l'intervalle à chaque fois, puis dans environ 95% des cas, au moins 90% de la distribution tombera dans l'intervalle. et 5% du temps, moins de 90% de la distribution sera contenue dans l'intervalle.
Michael R. Chernick
3
Lawless & Fredette (2005), "Frequentist Prediction Intervals and Predictive Distributions", Biometrika , 92 , 3 est une autre bonne référence, en plus de celles du lien que j'ai donné.
Scortchi - Réintégrer Monica

Réponses:

24

Ok, essayons ça. Je donnerai deux réponses - celle bayésienne, qui est à mon avis simple et naturelle, et l'une des possibles fréquentistes.

Solution bayésienne

Nous supposons un Beta prior sur , i, e., , car le modèle Beta-Binomial est conjugué, ce qui signifie que la distribution postérieure est également une distribution Beta avec les paramètres , (j'utilise pour dénoter le nombre de succès dans essais, au lieu de ). Ainsi, l'inférence est grandement simplifiée. Maintenant, si vous avez des connaissances préalables sur les valeurs probables de , vous pouvez l'utiliser pour définir les valeurs de et , c'est-à-dire pour définir votre préalable Beta, sinon vous pourriez supposer un préalable uniforme (non informatif), avecp ~ B e t un ( α , β ) α = α + k , β = β + n - k k n y p α β α = β = 1ppBeta(α,β)α^=α+k,β^=β+nkknypαβα=β=1, ou d'autres priors non informatifs (voir par exemple ici ). En tout cas, votre postérieur est

Pr(p|n,k)=Beta(α+k,β+nk)

Dans l'inférence bayésienne, tout ce qui compte est la probabilité postérieure, ce qui signifie qu'une fois que vous le savez, vous pouvez faire des inférences pour toutes les autres quantités de votre modèle. Vous voulez faire une inférence sur les observables : en particulier, sur un vecteur de nouveaux résultats , où n'est pas nécessairement égal à . Plus précisément, pour chaque , nous voulons calculer la probabilité d'avoir exactement succès dans les prochains essais, étant donné que nous avons obtenu succès dans les essais précédents ; la fonction de masse prédictive postérieure:y = y 1 , , y m m n j = 0 , , m j m k nyy=y1,,ymmnj=0,,mjmkn

Pr(j|m,y)=Pr(j|m,n,k)=01Pr(j,p|m,n,k)dp=01Pr(j|p,m,n,k)Pr(p|n,k)dp

Cependant, notre modèle binomial pour signifie que, conditionnellement à ayant une certaine valeur, la probabilité d'avoir succès dans essais ne dépend pas des résultats passés: c'est simplementp j mYpjm

f(j|m,p)=(jm)pj(1p)j

Ainsi, l'expression devient

Pr(j|m,n,k)=01(jm)pj(1p)jPr(p|n,k)dp=01(jm)pj(1p)jBeta(α+k,β+nk)dp

Le résultat de cette intégrale est une distribution bien connue appelée la distribution bêta-binomiale: en sautant les passages, nous obtenons l'expression horrible

Pr(j|m,n,k)=m!j!(mj)!Γ(α+β+n)Γ(α+k)Γ(β+nk)Γ(α+k+j)Γ(β+n+mkj)Γ(α+β+n+m)

Notre estimation ponctuelle pour , étant donné la perte quadratique, est bien sûr la moyenne de cette distribution, c.-à-d.j

μ=m(α+k)(α+β+n)

Maintenant, recherchons un intervalle de prédiction. Comme il s'agit d'une distribution discrète, nous n'avons pas d'expression de forme fermée pour , telle que . La raison en est que, selon la façon dont vous définissez un quantile, pour une distribution discrète, la fonction quantile n'est pas une fonction ou est une fonction discontinue. Mais ce n'est pas un gros problème: pour les petits , vous pouvez simplement écrire les probabilités et d'ici trouver tel que[j1,j2]Pr(j1jj2)=0.95mmPr(j=0|m,n,k),Pr(j1|m,n,k),,Pr(jm1|m,n,k)j1,j2

Pr(j1jj2)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)0.95

Bien sûr, vous trouveriez plus d'un couple, vous devriez donc idéalement rechercher le plus petit sorte que ce qui précède soit satisfait. Notez que[j1,j2]

Pr(j=0|m,n,k)=p0,Pr(j1|m,n,k)=p1,,Pr(jm1|m,n,k)=pm1

ne sont que les valeurs de la fonction CMF (fonction de masse cumulée) de la distribution bêta-binomiale, et en tant que telle, il existe une expression de forme fermée , mais c'est en termes de fonction hypergéométrique généralisée et donc assez compliquée. Je préfère simplement installer le package R extraDistret appeler pbbinompour calculer le CMF de la distribution bêta-binomiale. Plus précisément, si vous voulez calculer toutes les probabilités en une seule fois, écrivez simplement:p0,,pm1

library(extraDistr)  
jvec <- seq(0, m-1, by = 1) 
probs <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

alphaet betasont les valeurs des paramètres de votre a priori bêta, c'est-à-dire et (donc 1 si vous utilisez un a priori uniforme sur ). Bien sûr, tout serait beaucoup plus simple si R fournissait une fonction quantile pour la distribution bêta-binomiale, mais malheureusement ce n'est pas le cas.αβp

Exemple pratique avec la solution bayésienne

Soit , (ainsi nous avons initialement observé 70 succès dans 100 essais). Nous voulons une estimation ponctuelle et un intervalle de prédiction à 95% pour le nombre de succès dans les prochains essais . alorsn=100k=70jm=20

n <- 100
k <- 70
m <- 20
alpha <- 1
beta  <- 1

où j'ai supposé un a priori uniforme sur : selon les connaissances préalables pour votre application spécifique, cela peut ou non être un bon a priori. Doncp

bayesian_point_estimate <- m * (alpha + k)/(alpha + beta + n) #13.92157

De toute évidence, une estimation non entière pour n'a pas de sens, donc nous pourrions simplement arrondir à l'entier le plus proche (14). Ensuite, pour l'intervalle de prédiction:j

jvec <- seq(0, m-1, by = 1)
library(extraDistr)
probabilities <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

Les probabilités sont

> probabilities
 [1] 1.335244e-09 3.925617e-08 5.686014e-07 5.398876e-06
 [5] 3.772061e-05 2.063557e-04 9.183707e-04 3.410423e-03
 [9] 1.075618e-02 2.917888e-02 6.872028e-02 1.415124e-01
[13] 2.563000e-01 4.105894e-01 5.857286e-01 7.511380e-01
[17] 8.781487e-01 9.546188e-01 9.886056e-01 9.985556e-01

Pour un intervalle de probabilités égales, nous voulons que le plus petit tel que et le plus grand tel que . De cette façon, nous auronsj2Pr(jj2|m,n,k)0.975j1Pr(j<j1|m,n,k)=Pr(jj11|m,n,k)0.025

Pr(j1jj2|m,n,k)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)0.9750.025=0.95

Ainsi, en regardant les probabilités ci-dessus, nous voyons que et . La probabilité de cet intervalle de prédiction bayésien est de 0,9778494, ce qui est supérieur à 0,95. Nous pourrions trouver des intervalles plus courts tels que , mais dans ce cas, au moins une des deux inégalités pour les probabilités de queue ne serait pas satisfaite.j2=18j1=9Pr(j1jj2|m,n,k)0.95

Solution Frequentist

Je suivrai le traitement de Krishnamoorthy et Peng, 2011 . Soit et distribués indépendamment de façon binominale. Nous voulons un intervalle de prédiction pour , basée sur une observation de . En d'autres termes, nous recherchons tels que:YBinom(m,p)XBinom(n,p)12αYXI=[L(X;n,m,α),U(X;n,m,α)]

PrX,Y(YI)=PrX,Y(L(X;n,m,α)YU(X;n,m,α)]12α

Le " " est dû au fait que nous avons affaire à une variable aléatoire discrète, et donc nous ne pouvons pas nous attendre à obtenir une couverture exacte ... mais nous pouvons chercher un intervalle qui a toujours au moins le couverture nominale, donc un intervalle conservateur. Maintenant, il peut être prouvé que la distribution conditionnelle de étant donné est hypergéométrique avec la taille de l'échantillon , le nombre de succès dans la population et la taille de la population . Ainsi, le pmf conditionnel est12αXX+Y=k+j=ssnn+m

Pr(X=k|X+Y=s,n,n+m)=(nk)(msk)(m+ns)

Le CDF conditionnel de étant donné est doncXX+Y=s

Pr(Xk|s,n,n+m)=H(k;s,n,n+m)=i=0k(ni)(msi)(m+ns)

La première grande chose à propos de ce CDF est qu'il ne dépend pas de , que nous ne connaissons pas. La deuxième grande chose est qu'elle permet de trouver facilement notre PI: en fait, si nous avons observé une valeur de X, alors la limite de prédiction inférieure est le plus petit entier tel quepk1αL

Pr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α

en conséquence, la limite supérieure de prédiction est le plus grand entier tel que1α

Pr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α

Ainsi, est un intervalle de prédiction pour de couverture d'au moins . Notez que lorsque est proche de 0 ou 1, cet intervalle est conservateur même pour les grands , , c'est-à-dire que sa couverture est bien plus grande que .[L,U]Y12αpnm12α

Exemple pratique avec la solution Frequentist

Même réglage que précédemment, mais nous n'avons pas besoin de spécifier et (il n'y a pas de priors dans le framework Frequentist):αβ

n <- 100
k <- 70
m <- 20

L'estimation ponctuelle est maintenant obtenue en utilisant l'estimation MLE pour la probabilité de succès, , ce qui conduit à son tour à l'estimation suivante du nombre de succès dans essais:p^=knm

frequentist_point_estimate <- m * k/n #14

Pour l'intervalle de prédiction, la procédure est un peu différente. Nous recherchons le plus grand tel que , calculons donc l'expression ci-dessus pour tout dans :UPr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>αU[0,m]

jvec <- seq(0, m, by = 1)
probabilities <- phyper(k,n,m,k+jvec)

Nous pouvons voir que le plus grand tel que la probabilité est encore supérieure à 0,025 estU

jvec[which.min(probabilities > 0.025) - 1] # 18

Identique à l'approche bayésienne. La borne de prédiction inférieure est le plus petit entier tel que , DoncLPr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α

probabilities <- 1-phyper(k-1,n,m,k+jvec)
jvec[which.max(probabilities > 0.025) - 1] # 8

Ainsi, notre intervalle de prédiction "exact" fréquentiste est .[L,U]=[8,18]

DeltaIV
la source