Alpha And Beta Estimates For Beta Binomial And Beta Distributions

I am trying to fit my data to a beta-binomial distribution and estimate the alpha and beta shape parameters. For this distribution, the prior is taken from a beta distribution. Pyt

Solution 1:

Your problem is that fitting a beta-binomial model is just not the same as fitting a Beta model with the values equal to the ratios. I'm going to illustrate here with the bbmle package, which will fit similar models to VGAM (but with which I'm more familiar).

Preliminaries:

library("VGAM")  ## for dbetabinom.ab
x <- c(222,909,918,814,970,346,746,419,610,737,
       201,865,573,188,450,229,629,708,250,508)
y <- c(2,18,45,11,41,38,22,7,40,24,34,21,49,35,31,44,20,28,39,17)

library("bbmle")

Fit beta-binomial model:

mle2(y~dbetabinom.ab(size=x+y,shape1,shape2),
     data=data.frame(x,y),
     start=list(shape1=2,shape2=30))
## Coefficients:##    shape1    shape2##  1.736046 26.871526

This agrees more or less perfectly with the VGAM result you quote.

Now use the same framework to fit a Beta model instead:

mle2(y/(x+y) ~ dbeta(shape1,shape2),
     data=data.frame(x,y),
     start=list(shape1=2,shape2=30))
## Coefficients:##    shape1    shape2## 1.582021 24.060570

This fits your Python, beta-fit result. (I'm sure if you used VGAM to fit the Beta you'd get the same answer too.)

Solution 2:

You can use the conjugate_prior package for python

See code for a coin-flip example:

from conjugate_prior import BetaBinomial
heads = 95
tails = 105
prior_model = BetaBinomial() #Uninformative prior
updated_model = prior_model.update(heads, tails)
credible_interval = updated_model.posterior(0.45, 0.55)
print ("There's {p:.2f}% chance that the coin is fair".format(p=credible_interval*100))
predictive = updated_model.predict(50, 50)
print ("The chance of flipping 50 Heads and 50 Tails in 100 trials is {p:.2f}%".format(p=predictive*100))

Code taken from here