7

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.

TheChymera
  • 17,004
  • 14
  • 56
  • 86

3 Answers3

10

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

m00am
  • 5,910
  • 11
  • 53
  • 69
Danny
  • 450
  • 4
  • 10
5

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.

JulienD
  • 7,102
  • 9
  • 50
  • 84
  • is choose supposed to return an integer (I guess so, but [`numpy.choose`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.choose.html) seems to return an array). – TheChymera Nov 15 '14 at 01:18
  • Yes, choose is a combination, look there: http://en.wikipedia.org/wiki/Combination – JulienD Nov 15 '14 at 18:32
  • Ok, so maybe then you mean [comb](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.misc.comb.html), not [choose](http://docs.scipy.org/doc/numpy/reference/generated/numpy.choose.html)? am I right to use [`scipy.misc.comb`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.choose.html) for that? – TheChymera Nov 15 '14 at 19:18
  • Yes `scipy.misc.comb`. But many call it "n choose k". I'll update the answer with `comb`. – JulienD Nov 15 '14 at 20:59
  • Also the beta function B is `scipy.special.beta`. – JulienD Nov 16 '14 at 01:31
  • For future reference, [here is a very nice blog post](http://www.channelgrubb.com/blog/2015/2/27/beta-binomial-in-python) that gives a Python function for the beta-binomial PMF that still works well for large `a` and `b`. – Jake Westfall Oct 07 '17 at 18:47
1

The Beta-Binomial distribution is included in scipy from version 1.4.0 as scipy.stats.betabinom

vahndi
  • 1,055
  • 8
  • 16