Beta Binomial Function in Python

If your values of n (total # trials) and x (# successes) are large, then a more stable way to compute the beta-binomial probability is by working with logs. Using the gamma function expansion of the beta-binomial distribution function, the natural log of your desired probability is:

ln(answer) = gammaln(n+1) + gammaln(x+a) + gammaln(n-x+b) + gammaln(a+b) - \
        (gammaln(x+1) + gammaln(n-x+1) + gammaln(a) + gammaln(b) + gammaln(n+a+b))

where gammaln is the natural log of the gamma function, given in scipy.special.

(BTW: The loc argument just shifts the distribution left or right, which is not what you want here.)

Wiki says that the compound distribution function is given by

f(k|n,a,b) = comb(n,k) * B(k+a, n-k+b) / B(a,b)

where B is the beta function, a and b are the original Beta parameters and n is the Binomial one. k here is your x and p disappears because you integrate over the values of p to obtain this (convolution). That is, you won't find it in scipy but it is a one-liner provided you have the beta function from scipy.