0
$\begingroup$

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
        
    return total_nll

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?

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ – stats_model Jan 14 at 17:45
  • $\begingroup$ 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 $\endgroup$ – Jonathan F Jan 15 at 9:12
0
$\begingroup$

The problem is that I'm generating data using a beta distribution, if I use a beta binomial the fit works fine:

from scipy.stats import beta
data = [[round(r), 100] for r in betabinom.rvs(100, 70, 250, loc=0, size=1000, random_state=None)]
fmin(beta_bin_nll, [1, 1], args=(data,))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.