user3266890 - 9 months ago 200

Python Question

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?

Answer

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 `fit`

method.

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