I would like to calculate the probability given by a binomial distribution for predetermined x(successes), n(trials), and p(probability) - the later of which is given by a probability mass function Beta(a,b).
I am aware of scipy.stats.binom.pmf(x,n,p)
- but I am unsure how I can replace p with a probability function. I am also wondering whether I could use the loc
argument of scipy.stats.binom.pmf
to emulate this behaviour.
If your values of
n
(total # trials) andx
(# 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:where
gammaln
is the natural log of the gamma function, given inscipy.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
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.