I have been looking for a way to fit data to a beta binomial distribution and estimate alpha and beta, similar to the way the vglm package in VGAM library does. I have not been able to find how to do this in python. There is a scipy.stats.beta.fit() , but nothing for a beta binomial distribution. Is there a way to do this?
I have not seen estimation for beta-binomial in Python.
If you just want to estimate the parameters, then you can use scipy.optimize to minimize the log-likelihood function which you can write yourself or copy code after a internet search.
You can subclass
rv_discrete in order to use the framework of scipy.stats.distributions, but discrete distributions in scipy do not have a
If you want to use statsmodels, then you could subclass
GenericLikelihoodModel http://statsmodels.sourceforge.net/devel/dev/generated/statsmodels.base.model.GenericLikelihoodModel.html which uses scipy.optimize but defines most of the things we need for Maximum Likelihood estimations.
However, you need to write the code for the log-likelihood function. This would provide the usual maximum likelihood results such as standard errors for the parameters and various tests.
If you need beta-binomial regression, then the mean variance parameterization as used in the R package gamlss would be more common, and can reuse the
link functions to constrain the parameters to be in the valid domain.
As a related example: This is the gist with the GenericLikelihoodModel prototype that lead to a pull request for Beta-Regression for statsmodels: http://gist.github.com/brentp/089c7d6d69d78d26437f