如何获得 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 有一个通用的 class 用于最大似然估计 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
==============================================================================

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

How to get error estimates for fit parameters in scipy.stats.gamma.fit?

python

numpy

curve-fitting

scipy

data-fitting