Wie implementiere ich eine verallgemeinerte hypergeometrische Funktion zur Verwendung in Beta-Binomial-PDF, SF, PFF?

Ich schreibe eine Unterklasse scipy.stats._distn_infrastructure.rv_discretefür die Beta-Binomialverteilung, deren PMF ist

$$P(X=k \mid N, \alpha, \beta){N \choose k} \frac{\mathrm{B}(k+\alpha,N-k+\beta)}{\mathrm{B}(\alpha,\beta)},$$

Dabei ist $\mathrm{B}$ die Beta-Funktion. Meine derzeitige Implementierung von CDF und SF (Überlebensfunktion, entspricht 1 - CDF) ist ungenau; Die Strategie, die ich angewendet habe, berechnet den erwarteten Wert des Binomial-PDF in Bezug auf die Beta-Komponente:

$$P_{BB}(X \le k \mid N, \alpha, \beta) = E_p\left[P_{Binom}(X \le k \mid N, p)\right],$$ wobei . Ich erreiche dies mit der Methode, die nicht von Natur aus vektorisiert ist (sie stürzt auf etwas anderem als einem Float- oder 0d-Array ab).$p \sim \mathrm{Beta}(\alpha, \beta)$scipy.stats.beta.expect

Die PPF ist noch schlimmer - es ist eine Brute-Force-Schleife über die ganzen Zahlen so dass$k=0, \ldots, N$

$$P(X\le k \mid N, \alpha, \beta) \le q.$$

Laut Wikipedia ist die Überlebensfunktion für die Beta-Binomialverteilung

$$P(X > k \mid N, \alpha, \beta) = \frac{\mathrm{B}(\beta+n-k-1,\alpha+k+1)_3F_2(\boldsymbol{a},\boldsymbol{b};k)} {\mathrm{B}(\alpha,\beta)\mathrm{B}(n-k,k+2) (n+1)},$$

Dabei ist die verallgemeinerte hypergeometrische Funktion. Gibt es eine effiziente Möglichkeit, dies in Python zu berechnen, damit ich den Verweis auf entfernen kann ? Wie würde ich diese Funktion auch invertieren, um nach zu lösen, wenn ?${}_3F_2$beta.expect$k$$q=P(X \le k\mid N, \alpha, \beta)$

Dies beantwortet Ihre Frage nicht direkt, aber wenn Sie daran denken, die kumulative Verteilungsfunktion von Beta-Binomial effizienter abzuschätzen, können Sie einen rekursiven Algorithmus verwenden, der etwas effizienter ist als die naive Implementierung.

Beachten Sie die Wahrscheinlichkeitsmassenfunktion der Beta-Binomialverteilung

$$f(x) = {n \choose x} \frac{\mathrm{B}(x+\alpha, n-x+\beta)}{\mathrm{B}(\alpha, \beta)}$$

kann neu geschrieben werden, wenn Sie sich daran erinnern, dass undund dass , so dass es wird$\mathrm{B}(x,y)=\tfrac{\Gamma(x)\,\Gamma(y)}{\Gamma(x+y)}$$\Gamma(x) = (x-1)!$${n \choose k} = \prod_{i=1}^k \tfrac{n+1-i}{i}$

$$f(x) = \left( \prod_{i=1}^x \frac{n+1-i}{i} \right) \frac{\frac{(\alpha+x-1)!\,(\beta+n-x-1)!}{(\alpha+\beta+n-1)!}}{\mathrm{B}(\alpha,\beta)}$$

Dies erleichtert das Aktualisieren von auf$x$$x+1$

$$f(x\color{red}{+1}) = \left( \prod_{i=1}^x \frac{n+1-i}{i} \right) \color{red}{\frac{n+1-x+1}{x+1}} \frac{\frac{(\alpha+x-1)! \,\color{red}{(\alpha+x)}\,(\beta+n-x-1)! \, \color{red}{(\beta+n-x)^{-1}}}{(\alpha+\beta+n-1)!\,\color{red}{(\alpha+\beta+n)}}}{\mathrm{B}(\alpha,\beta)}$$

und damit können Sie die kumulative Verteilungsfunktion als berechnen

$$F(x) = \sum_{k=0}^x f(k)$$

Verwenden Sie nur einfache arithmetische Operationen, anstatt computerintensivere Funktionen zu berechnen.

Nebenbemerkung: Wenn Sie mit großen Zahlen arbeiten, treten Probleme mit der numerischen Genauigkeit auf, sodass robusterer Code mit Logarithmen arbeiten muss, obwohl Sie eine Verbesserung der Effizienz erwarten können (bis zu zwei- bis dreimal schnellerer Code, wenn ich nur wenige Benchmarks ausgeführt habe) C ++ - Code, der es implementiert, im Vergleich zur naiven Implementierung).