# How to estimate beta distribution parameters using a beta binomial with empirical bayes

I would like to estimate parameters for a beta distribution using a maximum likelihood approach in python (as mentioned here). I can do this for a beta:

from scipy.stats import beta

beta.fit(beta.rvs(a=70, b=250, loc=0, size=100), floc=0, fscale=1)


(74.75869456937754, 263.8103868963194, 0, 1)

If I then try to optimise parameters with respect to beta binomial likelihood I get the parameters tending towards infinity (although with the correct mean).

from scipy.stats import betabinom
from scipy.optimize import fmin

def beta_bin_nll(param, *args):
a, b = param

data = args[0]
total_nll = 0
for datum in data:
logpmf = betabinom.logpmf(k=datum[0],n=datum[1],a=a,b=b)
total_nll -= logpmf

data = [[round(r), 100] for r in beta.rvs(a=70, b=250, loc=0, scale=100, size=100)]
fmin(beta_bin_nll, [1, 1], args=(data,))


Optimization terminated successfully.
Current function value: 245.060220
Iterations: 164
Function evaluations: 360
array([4.72558414e+07, 1.72507488e+08])

How can I get this to correctly predict my beta parameters?

• Two points about your code are confusing to me. First, why are you rounding the beta draws to be integers? Second, if you're drawing from a beta distribution, why do you end up fitting an MLE for a beta-binomial distribution? If there is some more context around this, it might be useful to share. – stats_model Jan 14 at 17:45
• I've tidied up the code so that it's more understandable. I'm trying to follow this example, but I can't understand how to use the beta-binomial he mentions in the apendix – Jonathan F Jan 15 at 9:12

from scipy.stats import beta