Wie kann ich (numerisch) Werte für eine Beta-Verteilung mit großem Alpha & Beta approximieren?

11

Gibt es eine numerisch stabile Methode zur Berechnung der Werte einer Beta-Verteilung für Alpha, Beta (z. B. Alpha, Beta> 1000000)?

Eigentlich brauche ich nur ein 99% -Konfidenzintervall um den Modus, wenn das das Problem irgendwie einfacher macht.

Hinzufügen : Es tut mir leid, meine Frage war nicht so klar gestellt, wie ich dachte. Was ich tun möchte, ist Folgendes: Ich habe eine Maschine, die Produkte auf einem Förderband inspiziert. Ein Teil dieser Produkte wird von der Maschine abgelehnt. Wenn der Maschinenbediener nun eine Inspektionseinstellung ändert, möchte ich ihm die geschätzte Ausschussrate und einen Hinweis darauf geben, wie zuverlässig die aktuelle Schätzung ist.

Also dachte ich, ich behandle die tatsächliche Ablehnungsrate als Zufallsvariable X und berechne die Wahrscheinlichkeitsverteilung für diese Zufallsvariable basierend auf der Anzahl der zurückgewiesenen Objekte N und akzeptierten Objekte M. Wenn ich eine einheitliche vorherige Verteilung für X annehme, ist dies a Beta-Verteilung in Abhängigkeit von N und M. Ich kann diese Verteilung entweder direkt dem Benutzer anzeigen oder ein Intervall [l, r] finden, so dass die tatsächliche Ablehnungsrate in diesem Intervall mit p> = 0,99 (unter Verwendung der Shabbychef-Terminologie) liegt, und dieses anzeigen Intervall. Für kleines M, N (dh unmittelbar nach der Parameteränderung) kann ich die Verteilung direkt berechnen und das Intervall [l, r] approximieren. Für großes M, N führt dieser naive Ansatz jedoch zu Unterlauffehlern, da x ^ N * (1-x) ^ M zu klein ist, um als Float mit doppelter Genauigkeit dargestellt zu werden.

Ich denke, meine beste Wette ist es, meine naive Beta-Verteilung für kleine M, N zu verwenden und zu einer Normalverteilung mit demselben Mittelwert und derselben Varianz zu wechseln, sobald M, N einen bestimmten Schwellenwert überschreitet. Ist das sinnvoll?

Nikie
quelle
1
Möchten Sie die Mathematik oder einfach eine Codelösung in R oder so kennenlernen?
John
Ich muss dies in C # implementieren, damit die Mathematik gut ist. Ein Codebeispiel wäre auch in Ordnung, wenn es nicht auf einer eingebauten R / Matlab / Mathematica-Funktion beruht, die ich nicht in C # übersetzen kann.
Nikie
PDF, CDF oder inverse CDF?
JM ist kein Statistiker
Wenn Sie nicht auf Beta bestehen, können Sie die Kumaraswamy-Distribution verwenden, die sehr ähnlich ist und eine viel einfachere algebraische Form hat: en.wikipedia.org/wiki/Kumaraswamy_distribution
Tim

Antworten:

13

Eine normale Annäherung funktioniert sehr gut, besonders in den Schwänzen. Verwenden Sie einen Mittelwert von und eine Varianz von . Zum Beispiel erreicht der absolute relative Fehler in der Schwanzwahrscheinlichkeit in einer schwierigen Situation (in der die Schiefe von Belang sein könnte) wie Spitzenwert um und ist weniger als wenn Sie sich befinden mehr als 1 SD vom Mittelwert. (Dies liegt nicht daran, dass Beta so groß ist: Mit sind die absoluten relativen Fehler durchα βα/(α+β) α=106,αβ(α+β)2(1+α+β) 0,00026 0,00006 α = β = 10 6 0,0000001α=106,β=1080.000260.00006α=β=1060.0000001.) Somit ist diese Annäherung für im Wesentlichen jeden Zweck, der 99% -Intervalle umfasst, ausgezeichnet.

Beachten Sie angesichts der Änderungen an der Frage, dass man Beta-Integrale nicht durch tatsächliche Integration des Integranden berechnet: Natürlich kommt es zu Unterläufen (obwohl sie nicht wirklich wichtig sind, weil sie nicht wesentlich zum Integral beitragen). . Es gibt viele, viele Möglichkeiten, das Integral zu berechnen oder zu approximieren, wie in Johnson & Kotz (Verteilungen in der Statistik) dokumentiert. Einen Online-Rechner finden Sie unter http://www.danielsoper.com/statcalc/calc37.aspx . Sie brauchen tatsächlich die Umkehrung dieses Integrals. Einige Methoden zur Berechnung der Inversen sind auf der Mathematica-Website unter http://functions.wolfram.com/GammaBetaErf/InverseBetaRegularized/ dokumentiert.. Der Code wird in numerischen Rezepten (www.nr.com) bereitgestellt. Ein wirklich schöner Online-Rechner ist die Wolfram Alpha-Site (www.wolframalpha.com): Geben Sie inverse beta regularized (.005, 1000000, 1000001)für den linken Endpunkt und inverse beta regularized (.995, 1000000, 1000001)für den rechten Endpunkt ein ( , 99% Intervall).α=1000000,β=1000001

whuber
quelle
Perfekt! Ich hatte die ganze Zeit das NR-Buch auf meinem Schreibtisch, hätte aber nie gedacht, dort nachzuschauen. Danke vielmals.
Nikie
3

Ein kurzes grafisches Experiment legt nahe, dass die Beta-Verteilung einer Normalverteilung sehr ähnlich sieht, wenn sowohl Alpha als auch Beta sehr groß sind. Durch Googeln von "Beta Distribution Limit Normal" fand ich http://nrich.maths.org/discus/messages/117730/143065.html?1200700623 , was einen handwinkenden "Beweis" liefert.

Die Wikipedia-Seite für die Beta-Verteilung gibt den Mittelwert, den Modus (v nahe am Mittelwert für großes Alpha und Beta) und die Varianz an, sodass Sie eine Normalverteilung mit demselben Mittelwert und derselben Varianz verwenden können, um eine Annäherung zu erhalten. Ob es eine ausreichende Annäherung für Ihre Zwecke ist, hängt von Ihren Zwecken ab.

ein Stop
quelle
Dumme Frage: Wie haben Sie dieses grafische Experiment durchgeführt? Ich habe versucht, die Verteilung für Alpha / Beta um 100 zu zeichnen, konnte aber aufgrund von Unterlauffehlern nichts sehen.
Nikie
Sie möchten den Integranden nicht zeichnen: Sie möchten das Integral zeichnen. Sie können den Integranden jedoch auf viele Arten erhalten. Eine besteht darin, "Diagramm D (Beta (x, 1000000, 2000000), x) / Beta (1, 1000000, 2000000) von 0,3325 bis 0,334" am Wolfram Alpha-Standort einzugeben. Das Integral selbst wird mit "Plot Beta (x, 1000000, 2000000) / Beta (1, 1000000, 2000000) von 0,3325 bis 0,334" gesehen.
whuber
Ich habe den Integranden, dh das PDF der Beta-Distribution, in Stata aufgezeichnet - es hat eine eingebaute Funktion für das PDF. Für großes Alpha und Beta müssen Sie den Bereich des Diagramms einschränken, um zu sehen, dass es nahezu normal ist. Wenn ich es selbst programmieren würde, würde ich seinen Logarithmus berechnen und dann am Ende potenzieren. Das sollte bei den Unterlaufproblemen helfen. Die Beta-Funktion im Nenner wird als Gammafunktion definiert, die Fakultäten für ganzzahliges Alpha und Beta entspricht, und viele Pakete / Bibliotheken enthalten stattdessen lngamma () oder lnfactorial () sowie gamma () und factorial ().
Onestop
2

Ich werde daraus schließen, dass Sie ein Intervall wünschen so dass die Wahrscheinlichkeit, dass eine zufällige Ziehung aus dem Beta-RV erfolgt, im Intervall mit einer Wahrscheinlichkeit von 0,99 liegt, wobei die Bonuspunkte für und im Modus symmetrisch sind. Durch Gaußsche Ungleichung oder die Vysochanskii-Petunin-Ungleichung können Sie Intervalle konstruieren, die das Intervall und ziemlich anständige Näherungen wären. Für hinreichend große , werden Sie numerische Unterlauf Probleme haben in selbst darstellen und als unterschiedliche Zahlen, so dass diese Strecke gut genug sein kann.l r [ l , r ] α , β l r[l,r]lr[l,r]α,β lr

shabbychef
quelle
Wenn Alpha und Beta nicht zu weit voneinander entfernt sind (dh Alpha / Beta sind oben und unten begrenzt), ist die SD von Beta [Alpha, Beta] proportional zu 1 / Sqrt (Alpha). Zum Beispiel liegt die SD für alpha = beta = 10 ^ 6 sehr nahe bei 1 / Sqrt (8) / 1000. Ich denke, es wird kein Problem mit der Darstellung von l und r geben, selbst wenn Sie nur Floats mit einfacher Genauigkeit verwenden .
whuber
das heißt, dass nicht 'ausreichend groß' ist;)106
shabbychef
1
Ja, es ist eine verrückte Zahl für eine Beta-Anwendung. Übrigens erzeugen diese Ungleichungen überhaupt keine guten Intervalle, da sie über alle Verteilungen hinweg extrem sind (bestimmte Einschränkungen erfüllen).
whuber
@whuber: Du hast recht, das sind verrückte Zahlen. Mit meinem naiven Algorithmus waren die "gesunden" Zahlen einfach und funktionierten gut, aber ich konnte mir nicht vorstellen, wie ich sie für "verrückte" Parameter berechnen sollte. Daher die Frage.
Nikie
2
OK, Sie haben Recht: Sobald Alpha + Beta 10 ^ 30 überschreitet, werden Sie Schwierigkeiten mit Doppel haben :-). (Wenn Sie jedoch l und r als Unterschiede zum Mittelwert von Alpha / (Alpha + Beta) darstellen, ist alles in Ordnung, bis Alpha oder Beta etwa 10 ^ 303 überschreiten.)
whuber
1

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

Beispielsweise

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

erzeugt normalerweise eine Ausgabe wie

Zusammenfassung (replizieren (50, f (10000, 100, 1000000))) min. 1. Qu. Median Mean 3rd Qu. Max. 0,01205 0,10870 0,18680 0,24810 0,36170 0,68730

dh typische p-Werte liegen bei 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)))

produziert so etwas wie

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 

mit typischen p-Werten um 0,01

Die R- qqnormFunktion bietet auch eine hilfreiche Visualisierung, die ein sehr geradlinig aussehendes Diagramm für die Log-Odds-Verteilung erzeugt, das die ungefähre Normalität anzeigt. Die Verteilung der Beta-Dsitribute-Variablen erzeugt eine charakteristische Kurve, die die Nicht-Normalität anzeigt

Daher ist es sinnvoll, eine Gaußsche Näherung im Log-Odds-Bereich zu verwenden, selbst für stark verzerrte Werte, solange beide über 100 liegen.α,β

Daniel Mahler
quelle