如何獲得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 
==============================================================================