Beta binomial function in Python
I would like to calculate the probability obtained by the binomial distribution for predefined x (success), n (trial) and p (probability), the later of which is given by the probability mass function Beta (a, b).
I know scipy.stats.binom.pmf(x,n,p)
- but I'm not sure how I can replace p with a probability function. I'm also wondering if I can use an argument loc
scipy.stats.binom.pmf
to emulate this behavior.
source to share
The wiki says the complex distribution function is given
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 binomial. k here x and p disappears because you are integrating over the p values to get this (convolution). That is, you won't find it in scipy, but it is one liner if you have a beta function from scipy .
source to share
If your values n
(# of samples in total) and x
(# of successes) are large, then a more robust way to calculate the beta-binomial probability is to work with logs. Using the gamma extension 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 argument loc
just shifts the distribution left or right, which you don't like here.)
source to share