如何获得scipy.stats.gamma.fit中拟合参数的错误估计值？

Question

我有一些适合伽马分布使用scipy.stats。我能够提取形状，诠释和比例参数，他们看起来合理与我期望的数据范围。

我的问题是：有没有办法让参数中的错误呢？像curve_fit的输出。注：我不直接使用曲线拟合，因为它不能正常工作，大部分时间都无法计算伽玛分布的参数。另一方面，scipy.stats.gamma.fit可以正常工作。

这是我正在做的一个例子（没有我的实际数据）。

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? 

在此先感谢

Answer 1

编辑 警告：下面举例说明以下问题中的例子使用GenericLikelihoodModel。 然而，在伽马分布的情况下，位置参数改变了对最大似然估计的一般假设排除的分布的支持。更标准的用法将修复支持，floc = 0，所以它只是一个双参数分布。在那种情况下，标准的MLE理论适用。

Statsmodels有一个最大似然估计的泛型类，GenericLikelihoodModel。它不是直接为这种情况设计的，但可以使用一些帮助（定义属性并提供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()) 

这将打印基于最大似然估计以下

(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 
============================================================================== 

其他结果然后也可用，例如一个Z测试，所述第一参数为10可以像通过指定执行一个限制矩阵，或通过使用一个字符串表达式与等式：

>>> 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 
==============================================================================