2

我有一个具有三个参数的模型，即：beta > 0、gamma > 0 和 phi > 0。给定数据 T，我能够评估 3 个参数的 MLE。

但是，我还要验证 (i) gamma = 1; 的情况，这一点非常重要。(ii) beta = 1 和 (iii) gamma = 1，beta = 1。

因此，我希望能够在情况 (i)、(ii) 和 (iii) 下评估 MLE。

我确信有一种方法可以轻松评估它，而无需在 StatsModels 中为 GenericLikelihoodModel 编写 4 个不同的类，但我就是想不通。

我用于具有 3 个参数的通用模型的代码如下，它运行良好，因为我正在使用本文中的数据来验证实现。

import math
import numpy
from statsmodels.base.model import GenericLikelihoodModel

def log_likelihood_WPLP(T, beta, gamma, phi):
    n = len(T) - 1
    ret = n * (numpy.log(phi) + numpy.log(beta) + numpy.log(gamma)) + \
          (beta - 1) * numpy.log(T[1]) + (gamma - 1) * numpy.log( T[1] ** beta) - phi * T[1] ** (beta * gamma) + \
          sum([(beta - 1) * numpy.log(T[i]) + (gamma - 1) * numpy.log(T[i] ** beta - T[i-1] ** beta) - \
          phi * (T[i] ** beta - T[i-1] ** beta) ** gamma for i in range(2, n+1)])
    return ret

class WeibullPowerLawTRP(GenericLikelihoodModel):
    def __init__(self, endog, exog=None, **kwds):
        if exog is None:
            exog = numpy.zeros_like(endog)

        super(WeibullPowerLawTRP, self).__init__(endog, exog, **kwds)

    def nloglikeobs(self, params):
        beta, gamma, phi = numpy.exp(params[0]), numpy.exp(params[1]), numpy.exp(params[2])
        return -log_likelihood_WPLP(self.endog, beta=beta, gamma=gamma, phi=phi)

    def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
        if start_params is None:
            start_params = numpy.array([-.5, -.5, -.5])

        return super(WeibullPowerLawTRP, self).fit(start_params=start_params,
                                                    maxiter=maxiter, maxfun=maxfun, **kwds)

T = [   0,    1,    4,  305,  330,  651,  856,  996, 1016, 1155, 1520, 1597, 
     1729, 1758, 1852, 2070, 2073, 2093, 2213, 3197, 3555, 3558, 3724, 3768, 
     4103, 4124, 4170, 4270, 4336, 4416, 4492, 4534, 4578, 4762, 5474, 5573, 
     5577, 5715, 6424, 6692, 6830, 6999]

model = WeibullPowerLawTRP(T, model='WPLP')
result = model.fit()
print(result.summary())

beta, gamma, phi = numpy.exp(result.params)
alpha = phi ** (1/gamma) / (math.gamma(1+1/gamma))
print(f"α = {alpha:.4f}\n"
      f"β = {beta:.4f}\n"
      f"γ = {gamma:.4f}\n")

提前致谢！干杯，亚历山大

4

0 回答 0