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