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

Question

我有一些数据，使用scipy.stats拟合成伽马分布。我能够提取形状、位置和缩放参数，并且它们看起来符合我的数据范围。我的问题是：是否有办法获取参数的误差？类似于curve_fit的输出。注意：我不直接使用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?

Thanks in advance

- iluvatar

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

- Josef

---

非常感谢，我一直在使用 http://scipy-central.org/item/36/2/error-estimates-for-fit-parameters-resulting-from-least-squares-fits-using-bootstrap-resampling 中描述的bootstrap方法，但是您提供的答案似乎更完整，可以让我应用更多的测试。不过，我想知道为什么使用这种方法，在运行相同脚本多次后，会得到（有时非常）不同的参数结果。这正常吗？例如，c0的值在8到20之间变化（当然，标准误差也会改变）。再次感谢。 - iluvatar

注意，我使用了scipy拟合结果作为GenericLikelihoodModel的起始参数。一般来说，这可能没有很好的默认起始值，在某些情况下，目标函数的曲率也不是“好的”。我不知道三参数Gamma分布的行为如何。也许“basinhopping”可以为优化提供成功的全局最小值。 - Josef

另一个一般性的评论：我猜通常我们会将“loc”设置为固定值零来模拟正值数据。最大似然经常或一般情况下在估计分布的支持时存在问题。 - Josef

感谢您的评论。在我的特定情况下，我有一个与零不同的loc（实际上，在我的特定理论框架中，这是一个非常重要的参数），所以我不能将其设置为零。尽管如此，我认为您的答案是合适的，并且我会选择它作为最佳答案。 - iluvatar

你好，能帮我解答一下我的问题吗？谢谢。http://stackoverflow.com/questions/44022955/why-are-the-parameters-of-weibull-not-unique-for-a-given-data - yanachen