如何在scipy.stats.gamma.fit中获取拟合参数的误差估计？

Question

我有一些我适合使用scipy.stats的gamma分布.我能够提取形状,位置和比例参数,它们看起来与我期望的数据范围合理.

我的问题是:有没有办法在参数中得到错误？类似于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
==============================================================================

然后,基于最大似然估计的其他结果也是可用的,例如,可以通过指定限制矩阵或通过使用具有相等性的字符串表达式来执行第一参数为10的z测试:

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