¿Están Wikipedia y Wolfram Mathworld equivocados sobre la FDA de la binomial beta?

Question

Wikipedia y Mathworld afirmar que la FDA de la distribución beta-binomial es: $$1- \tfrac{\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)}$$ donde $_3F_2(\boldsymbol{a},\boldsymbol{b};k)$ es la función hipergeométrica generalizada $$_3F_2(1,\! \alpha\! +\! k\!+ \! 1,k\! -\! n\! +\! 1;k\! +\! 2,k\! +\! 2\! -\! \beta\! -\! n;1)\!$$

Esto es dudoso porque $_3F_2(\cdot,\cdot,\cdot;\cdot,\cdot;1)$ es un singularidad . Usando Sympy, evalué esa función hipergeométrica generalizada con el último argumento cambiado a diferentes valores cerca de $1$ y esto es lo que tengo:

In [5]: import mpmath as mp

In [7]: n = 1

In [8]: k = 1

In [9]: alpha = 1

In [10]: beta = 1

In [16]: spec = lambda alpha, beta, n, k: mp.hyp3f2(1, alpha + k + 1, k - n + 1, k + 2, k + 2 - beta 
    ...: - n, '1.0001')

In [17]: spec(alpha, beta, n, k)
Out[17]: mpf('-10000.0000000011')

In [18]: spec = lambda alpha, beta, n, k: mp.hyp3f2(1, alpha + k + 1, k - n + 1, k + 2, k + 2 - beta 
    ...: - n, '0.99999999')

In [19]: spec(alpha, beta, n, k)
Out[19]: mpf('99999999.497524068')

Así que incluso cambia de signo. ¿Cuál es la FDA de la binomial beta?

Answer 1

Basándome en lo que dijo gammatester, visité http://reference.wolfram.com/language/ref/BetaBinomialDistribution.html . Allí la FCD se define a trozos con la $k \geq n$ que se define por separado del $0\leq k<n$ caso. En particular, el GHGF no debe utilizarse cuando $k=n$ . He editado el artículo de Wikipedia para dar la fórmula completa. Esa fórmula es $$\Pr[K\leq k \mid n,\alpha,\beta]=\begin{cases} 0,& k < 0 \\ 1- \tfrac{\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)}, & 0 \leq k < n \\ 1,& k \geq n \end{cases} $$

Si $k$ no es un número entero, entonces es necesario sustituir $k$ con $\lfloor k \rfloor$ en el lado derecho.

---

Lo que sigue son pruebas empíricas de que la fórmula dada funciona. Calcula la función de supervivencia (es decir $\Pr[K \geq k]$ ) en lugar de la FCD (es decir $\Pr[K\leq k]$ ).

python import mpmath as mp from scipy import * from scipy.special import * from sympy.functions import binomial

La fórmula inteligente que utiliza funciones especiales

python def beta_binomial_survival_fn(alpha, beta, n, k): if k <= n: return exp(betaln(beta + n - k, alpha + k) + log(float(mp.hyp3f2(1, alpha + k, k - n, k + 1, k + 1 - beta - n, 1))) - betaln(alpha, beta) - betaln(n - k + 1, k + 1) - log(n+1)) elif k > n: return 0 elif k < 0: return 1 else: raise ValueError("Can't compare k to n")

La fórmula menos inteligente que suma las masas de probabilidad

python def beta_binomial_survival_fn_less_clever(alpha, beta, n, k): return sum(binomial(n, i) * exp(betaln(i+alpha, n-i+beta) - betaln(alpha,beta)) for i in range(k, n+1))

A continuación se calcula el error máximo al utilizar la fórmula inteligente en lugar de la fórmula no tan inteligente

python max(abs(beta_binomial_survival_fn(alpha, beta, n, k) - beta_binomial_survival_fn_less_clever(alpha, beta, n, k)) for alpha in range(1,6) for beta in range(1,6) for n in range(1,6) for k in range(0,6))

8.88178419700125e-16