Maximum Likelihood Method for Gamma Distribution

Question

I used Method of Moments to fit negative binomial, poisson and geometric distributions.

Basic idea: get empirical first, second, etc. moments, then derive distribution parameters from these moments.

You can see the details in this question: Fitting Distributions with Maximum Likelihood Method

And now i want to implement this method for gamma distribution;

For Gamma distribution i applied this;

import pandas as pd
from scipy.stats import gamma

x = pd.Series(x)
mean = x.mean()
var = x.var()
likelihoods = {}

alpha = (mean**2)/var
beta = alpha / mean
likelihoods['gamma'] = x.map(lambda val: gamma.pdf(val, alpha)).prod()

However, the likelihood value is infinite in the results for Gamma Distribution. So, I'm not sure I can apply it correctly this method for Gamma.

Is there anyone who can help for that?

Answer 1

You're using definition of the Gamma distribution with \alphaand \beta, while NumPy and SciPy are using shape and scale parameters, which are k and \theta.

Basically, you have to reciprocate \beta to get scale back.

import numpy as np
import pandas as pd
from scipy.stats import gamma

shape = 2.0
scale = 1.3

x = pd.Series(gamma.rvs(shape, loc = 0.0, scale=scale, size = 100000))
#x = pd.Series(np.random.gamma(shape, scale=scale, size=100000))

mean = x.mean()
var  = x.var()

alpha = (mean**2)/var
beta  = alpha / mean

print((alpha, 1.0/beta)) # NB! scale is 1/beta

likelihoods = {}
likelihoods['gamma'] = x.map(lambda val: gamma.pdf(val, alpha, loc = 0.0, scale = 1.0/beta)).prod()
print(likelihoods['gamma'])