Fit beta binomial

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 1

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

Answer 2

This python module provides that https://github.com/lfiaschi/fastbetabino