2

I have a model with three parameters, namely: beta > 0, gamma > 0 and phi > 0. Given the data T, I am able to evaluate the MLEs for the 3 parameters.

However, it is very important that I also verify the cases where (i) gamma = 1; (ii) beta = 1 and (iii) gamma = 1, beta = 1.

Therefore, I want to be able to evaluate the MLE under cases (i), (ii) and (iii).

I am certain that there is a way of easily evaluating it, without having to write 4 different classes for the GenericLikelihoodModel in StatsModels, but I just cannot figure it out.

The code I am using for the general model with 3 parameters is as follows, and it is working perfectly, given that I am using the data from this paper to validate the implementation.

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")

Thanks in advance! Cheers, Alexandre

  • you can hardcode the restrictions in the model, or use expand_params method. There is an old example for optionally fixing df in a t-distribution model. https://github.com/statsmodels/statsmodels/blob/master/statsmodels/miscmodels/tmodel.py#L50 Examples for it are in the unit tests – Josef May 12 '20 at 16:48

0 Answers0