时间戳

首页 > 解决方案 > 如何在python中获得2参数weibull分布的渐近协方差矩阵

python - 如何在python中获得2参数weibull分布的渐近协方差矩阵

问题描述

我正在尝试用 python 实现一个 matlab 代码。

在 Matlab 中，weibull 分布的置信区间得到如下：

CI = 0.95 % 
p = [scale_hat shape_hat];
[nlogl,avar]=wbllike(p,y_ordered);
[qmid,qlo,qup]=wblinv(1-1/nCount,p(1),p(2),avar,1-CI);

在 scipy 中，我找不到 wbllike 函数的等价物。好像没有实现这个功能。相反，我尝试自己实现 wblike 功能。到目前为止，负对数似然可以很容易地获得如下：

def weibull_neg_loglike(x, shape, scale):
    return -stats.weibull_min.logpdf(x, shape_param, scale=scale_param).sum()

我发现python中有statsmodels包，它提供了GeneralLikelyhoodModel。但不知道如何让事情正常进行。我研究了一点网站，http://rlhick.people.wm.edu/posts/estimating-custom-mle.html#statsmodels但仍然是徒劳的。

任何建议都会有所帮助。

标签： pythonscipyweibull

解决方案

以下是我不久前写的一个笔记本的一部分。它的目标是审查 weibull 回归，因此它需要一个明确指定的 exog，例如一组。

它还包括一个在我的旧 scipy 版本中不起作用的 logsf 解决方法。

    def weibull_min_logsf(x, c, scale=1):
        x = x / scale
        return -np.power(x, c)

    from statsmodels.base.model import GenericLikelihoodModel

    class WeibullModel(GenericLikelihoodModel):

        def __init__(self, endog, exog, censored=None, **kwds):
            super(WeibullModel, self).__init__(endog, exog, **kwds)
            if censored is not None:
                self.censored = censored.astype(int)
            else:
                self.censored = np.zeros(len(self.endog), np.int)

            self.k_params = self.exog.shape[1] + 1


        def loglike(self, params):
            params_ex = params[:-1]
            params_shape = params[-1]
            m = self.exog.dot(params_ex)
            llf = (1 - self.censored) * stats.weibull_min.logpdf(self.endog, params_shape, scale=m)
            #stats.weibull_min.logsf overflows in my older version of scipy, changed in newer versions
            #llf2 = self.censored * stats.weibull_min.logsf(self.endog, params_shape, scale=m)
            llf += self.censored * weibull_min_logsf(self.endog, params_shape, scale=m)
            return llf.sum()

        def _get_distribution(self, params, exog):
            """similar to a predict method
            """
            params_ex = params[:-1]
            params_shape = params[-1]
            m = exog.dot(params_ex)
            return stats.weibull_min(params_shape, scale=m)

    mod = WeibullModel(data['duration'], np.ones(len(data)))
    res = mod.fit(start_params=np.ones(2)) 
    print(res.summary())

推荐阅读