I have some which I am fitting to the gamma distribution using scipy.stats. I am able to extract the shape, loc, and scale params and they look reasonable with the data ranges I expect.
My question is: is there a way to also get the errors in the parameters? something like the output of curve_fit. NOTE: I don't use curve fit directly because it is not working properly and most of the time is not able to compute the parameters of the gamma distribution. On the other hand, scipy.stats.gamma.fit works ok.
This is an example (no my actual data) of what I am doing.
from scipy.stats import gamma
shape = 12; loc = 0.71; scale = 0.0166
data = gamma.rvs(shape, loc=loc, scale=scale, size=1000)
params = gamma.fit(data) # params close to but not the same as (shape, loc, scale)
# HOW TO ESTIMATE/GET ERRORS FOR EACH PARAM?
Thanks in advance
edit
Warning: The following illustrates the use of GenericLikelihoodModel following the example in the question.
However, in the case of the gamma distribution the location parameter shifts the support of the distribution which is ruled out by the general assumptions for maximum likelihood estimation. The more standard usage would fix the support, floc=0, so it is just a two-parameter distribution. In that case, standard MLE theory applies.
Statsmodels has a generic class for maximum likelihood estimation, GenericLikelihoodModel. It's not directly designed for this case, but can be used with some help (defining attributes and providing start_params).
import numpy as np
from statsmodels.base.model import GenericLikelihoodModel
from scipy.stats import gamma
shape = 12; loc = 0.71; scale = 0.0166
data = gamma.rvs(shape, loc=loc, scale=scale, size=1000)
params = gamma.fit(data) # params close to but not the same as (shape, loc, scale)
# HOW TO ESTIMATE/GET ERRORS FOR EACH PARAM?
print(params)
print('\n')
class Gamma(GenericLikelihoodModel):
nparams = 3
def loglike(self, params):
return gamma.logpdf(self.endog, *params).sum()
res = Gamma(data).fit(start_params=params)
res.df_model = len(params)
res.df_resid = len(data) - len(params)
print(res.summary())
This prints the following
(10.31888758604304, 0.71645502437403186, 0.018447479022445423)
Optimization terminated successfully.
Current function value: -1.439996
Iterations: 69
Function evaluations: 119
Gamma Results
==============================================================================
Dep. Variable: y Log-Likelihood: 1440.0
Model: Gamma AIC: -2872.
Method: Maximum Likelihood BIC: -2852.
Date: Sun, 12 Jul 2015
Time: 04:00:05
No. Observations: 1000
Df Residuals: 997
Df Model: 3
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
par0 10.3187 2.242 4.603 0.000 5.925 14.712
par1 0.7165 0.019 37.957 0.000 0.679 0.753
par2 0.0184 0.002 8.183 0.000 0.014 0.023
==============================================================================
Other results based on the maximum likelihood estimates are then also available, for example a z-test that the first parameter is 10 can be performed like either by specifying a restriction matrix or by using a string expression with the equality:
>>> res.t_test(([1, 0, 0], [10]))
<class 'statsmodels.stats.contrast.ContrastResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
c0 10.3187 2.242 0.142 0.887 5.925 14.712
==============================================================================
>>> res.t_test('par0=10')
<class 'statsmodels.stats.contrast.ContrastResults'>
Test for Constraints
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
c0 10.3187 2.242 0.142 0.887 5.925 14.712
==============================================================================