Comment puis-je (numériquement) approcher des valeurs pour une distribution bêta avec un grand alpha et bêta

11

Existe-t-il un moyen numériquement stable de calculer les valeurs d'une distribution bêta pour les grands nombres alpha, bêta (par exemple alpha, bêta> 1000000)?

En fait, je n'ai besoin que d'un intervalle de confiance de 99% autour du mode, si cela rend le problème plus facile.

Ajouter : Je suis désolé, ma question n'était pas aussi clairement formulée que je le pensais. Ce que je veux faire, c'est ceci: j'ai une machine qui inspecte les produits sur un tapis roulant. Une partie de ces produits est rejetée par la machine. Maintenant, si l'opérateur de la machine modifie un paramètre d'inspection, je veux lui montrer le taux de rejet estimé et quelques indices sur la fiabilité de l'estimation actuelle.

J'ai donc pensé traiter le taux de rejet réel comme une variable aléatoire X, et calculer la distribution de probabilité pour cette variable aléatoire en fonction du nombre d'objets rejetés N et d'objets acceptés M. Si je suppose une distribution préalable uniforme pour X, ceci est un distribution bêta en fonction de N et M. Je peux soit afficher cette distribution directement à l'utilisateur, soit trouver un intervalle [l, r] afin que le taux de rejet réel soit dans cet intervalle avec p> = 0,99 (en utilisant la terminologie de shabbychef) et afficher ceci intervalle. Pour les petits M, N (c'est-à-dire immédiatement après le changement de paramètre), je peux calculer la distribution directement et approximer l'intervalle [l, r]. Mais pour les grands M, N, cette approche naïve conduit à des erreurs de sous-dépassement, car x ^ N * (1-x) ^ M est trop petit pour être représenté comme un flotteur à double précision.

Je suppose que mon meilleur pari est d'utiliser ma distribution bêta naïve pour les petits M, N et de passer à une distribution normale avec la même moyenne et la même variance dès que M, N dépasse un certain seuil. Cela a-t-il du sens?

nikie
la source
1
Voulez-vous connaître les mathématiques ou simplement une solution de code en R ou quelque chose comme ça?
John
J'ai besoin de l'implémenter en C #, donc les mathématiques seraient bonnes. Un exemple de code serait bien aussi, s'il ne repose pas sur une fonction R / Matlab / Mathematica intégrée que je ne peux pas traduire en C #.
nikie
PDF, CDF ou CDF inverse?
JM n'est pas statisticien
Si vous n'insistez pas sur la bêta, vous pouvez utiliser une distribution Kumaraswamy qui est très similaire et a une forme algébrique beaucoup plus simple: en.wikipedia.org/wiki/Kumaraswamy_distribution
Tim

Réponses:

13

α/(α+β)αβ(α+β)2(1+α+β)α=106,β=1080.000260.00006α=β=1060.0000001.) Ainsi, cette approximation est excellente pour pratiquement tout objectif impliquant des intervalles de 99%.

À la lumière des modifications apportées à la question, notez que l'on ne calcule pas les intégrales bêta en intégrant réellement l'intégrand: bien sûr, vous obtiendrez des sous-flux (bien qu'ils n'aient pas vraiment d'importance, car ils ne contribuent pas de manière appréciable à l'intégrale) . Il existe de très nombreuses façons de calculer l'intégrale ou de l'approcher, comme indiqué dans Johnson & Kotz (Distributions in Statistics). Une calculatrice en ligne se trouve à http://www.danielsoper.com/statcalc/calc37.aspx . Vous avez en fait besoin de l'inverse de cette intégrale. Certaines méthodes de calcul de l'inverse sont documentées sur le site Mathematica à http://functions.wolfram.com/GammaBetaErf/InverseBetaRegularized/inverse beta regularized (.005, 1000000, 1000001)inverse beta regularized (.995, 1000000, 1000001)α=1000000,β=1000001

whuber
la source
Parfait! J'avais le livre NR sur mon bureau tout le temps, mais je n'ai jamais pensé y regarder. Merci beaucoup.
nikie
3

Une expérience graphique rapide suggère que la distribution bêta ressemble beaucoup à une distribution normale lorsque alpha et bêta sont tous deux très importants. En recherchant «limite de distribution bêta normale» sur Google, j'ai trouvé http://nrich.maths.org/discus/messages/117730/143065.html?1200700623 , ce qui donne une «preuve» d'ondulation.

La page wikipedia de la distribution bêta donne sa moyenne, son mode (v proche de la moyenne pour les grands alpha et bêta) et la variance, vous pouvez donc utiliser une distribution normale avec la même moyenne et variance pour obtenir une approximation. Que ce soit une approximation suffisamment bonne pour vos besoins dépend de vos objectifs.

un arrêt
la source
Question stupide: Comment avez-vous fait cette expérience graphique? J'ai essayé de tracer la distribution de l'alpha / bêta autour de 100, mais je n'ai rien vu en raison d'erreurs de sous-dépassement.
nikie
Vous ne voulez pas tracer l'intégrale: vous voulez tracer l'intégrale. Cependant, vous pouvez obtenir l'intégrale de plusieurs manières. L'une consiste à saisir «tracé D (bêta (x, 1000000, 2000000), x) / bêta (1, 1000000, 2000000) de 0,3325 à 0,334» sur le site Wolfram Alpha. L'intégrale elle-même est vue avec "Plot beta (x, 1000000, 2000000) / beta (1, 1000000, 2000000) de 0,3325 à 0,334".
whuber
J'ai tracé l'intégrande, c'est-à-dire le pdf de la distribution bêta, dans Stata - il a une fonction intégrée pour le pdf. Pour les grands alpha et bêta, vous devez restreindre la plage de l'intrigue pour voir qu'elle est proche de la normale. Si je le programmais moi-même, je calculerais son logarithme, puis j'exposerais à la fin. Cela devrait aider à résoudre les problèmes de sous-dépassement. La fonction bêta du dénominateur est définie en termes de fonctions gamma, équivalentes aux factorielles pour les nombres alpha et bêta entiers, et de nombreux packages / bibliothèques incluent à la place lngamma () ou lnfactorial () / ainsi que les fonctions gamma () et factorielle ().
onestop
2

[l,r]lr[l,r]α,β lr en tant que nombres distincts, donc cette route peut être assez bonne.

shabbychef
la source
Lorsque l'alpha et la bêta ne sont pas trop éloignés l'un de l'autre (c.-à-d. Que l'alpha / bêta sont délimités au-dessus et en dessous), l'écart-type de la bêta [alpha, bêta] est proportionnel à 1 / sqq (alpha). Par exemple, pour alpha = beta = 10 ^ 6, la SD est très proche de 1 / Sqrt (8) / 1000. Je pense qu'il n'y aura pas de problème avec la représentation de l et r même si vous n'utilisez que des flotteurs simple précision .
whuber
106
1
Oui, c'est un chiffre fou pour une application bêta. BTW, ces inégalités ne produiront pas du tout de bons intervalles, car elles sont extrêmes sur toutes les distributions (satisfaisant certaines contraintes).
whuber
@whuber: Vous avez raison, ce sont des chiffres fous. Avec mon algorithme naïf, les nombres "sains" étaient faciles et fonctionnaient bien, mais je ne pouvais pas imaginer comment les calculer pour des paramètres "fous". D'où la question.
nikie
2
OK, vous avez raison: une fois que l'alpha + bêta dépasse 10 ^ 30 environ, vous aurez des difficultés avec les doubles :-). (Mais si vous représentez l et r comme des différences par rapport à la moyenne de alpha / (alpha + beta), tout ira bien jusqu'à ce que alpha ou beta dépasse environ 10 ^
303.
1

pplog(p/(1p))min(α,β)>100

Par exemple

f <- function(n, a, b) {
    p <- rbeta(n, a, b)
    lor <- log(p/(1-p))
    ks.test(lor, 'pnorm', mean(lor), sd(lor))$p.value
}
summary(replicate(50, f(10000, 100, 1000000)))

produit généralement une sortie comme

résumé (réplique (50, f (10000, 100, 1000000))) Min. 1er Qu. Médiane Moyenne 3e Qu. Max. 0,01205 0,10870 0,18680 0,24810 0,36170 0,68730

c'est-à-dire que les valeurs de p typiques sont d'environ 0,2.

α=100,β=100000

p

f2 <- function(n, a, b) {
    p <- rbeta(n, a, b)
    ks.test(p, 'pnorm', mean(p), sd(p))$p.value
}
summary(replicate(50, f2(10000, 100, 1000000)))

produit quelque chose comme

summary(replicate(50, f2(10000, 100, 1000000)))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
2.462e-05 3.156e-03 7.614e-03 1.780e-02 1.699e-02 2.280e-01 

avec des valeurs de p typiques autour de 0,01

La qqnormfonction R donne également une visualisation utile, produisant un tracé très direct pour la distribution log-odds indiquant la normalité approximative la distribution de la variable beta dsitribute produit une courbe distinctive indiquant la non normalité

α,β

Daniel Mahler
la source