Are Wikipedia and Wolfram Mathworld wrong about the CDF of the beta binomial?

Question

Wikipedia and Mathworld claim that the CDF of the beta-binomial distribution is: $$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)}$$ where $_3F_2(\boldsymbol{a},\boldsymbol{b};k)$ is the generalized hypergeometric function $$_3F_2(1,\! \alpha\! +\! k\!+ \! 1,k\! -\! n\! +\! 1;k\! +\! 2,k\! +\! 2\! -\! \beta\! -\! n;1)\!$$

This is dubious because $_3F_2(\cdot,\cdot,\cdot;\cdot,\cdot;1)$ is a singularity. Using Sympy, I evaluated that generalized hypergeometric function with the last argument changed to different values near $1$, and here's what I got:

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

So it even changes sign. What is the CDF of the beta binomial?

Answer 1

Based on what gammatester said, I visited http://reference.wolfram.com/language/ref/BetaBinomialDistribution.html. There the CDF is defined piecewise with the $k \geq n$ case being defined separately from the $0\leq k<n$ case. In particular, the GHGF should not be used when $k=n$. I edited the Wikipedia article to give the complete formula. That formula is $$\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} $$

If $k$ is not an integer then it's necessary to replace $k$ with $\lfloor k \rfloor$ on the RHS.

---

The following is empirical evidence that the given formula works. It computes the survival function (i.e. $\Pr[K \geq k]$) instead of the CDF (i.e. $\Pr[K\leq k]$).

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

The clever formula using special functions

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

The less clever formula that sums the probability masses

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

The following computes the maximum error in using the clever formula in place of the not-so-clever formula

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