Error in maximum likelihood estimation using R

Question

I am currently finishing my undergraduate thesis about parameter estimation on Exponentiated Modified Weibull Extension (EMWE) distribution introduced by Sarhan and Apaloo (2013) with the following pdf:

f(x,theta)=theta[1]*theta[2]*theta[3]*((x/theta[4])^(theta[2]-1))*(exp(((x/theta[4])^theta[2])+(theta[1]*theta[4]*(1-(exp(x/theta[4])^theta[2])))))*(1-(exp(theta[1]*theta[4]*(1-(exp(x/theta[4])^theta[2])))))^(theta[3]-1)

This distribution has four parameters to be estimated using maximum likelihood estimation. Due to estimate parameters implicitly then I tried to continue with the Newton-Raphson iteration method. For my computation process, I'm using statistical software R Language. Package that I use is "maxLik" with an initial value for Newton-Raphson method is (theta [1] = 0.00007181, theta [2] = 3,148, theta [3] = 0.145, theta [4] = 49.05).

This is the loglikelihood function:

l(theta)=n*(log(theta[1])+log(theta[2])+log(theta[3])+(1-theta[2])*log(theta[4])+theta[1]*theta[4])+(theta[2]-1)*sum(log(xi))+(1/(theta[4]^theta[2]))*sum(xi^theta[2])-(theta[1]*theta[4])*sum(exp((xi/theta[4])^theta[2]))+(theta[3]-1)*sum(1-(exp((theta[1]*theta[4])*(1-(exp((xi/theta[4])^theta[2]))))))

But in this process of parameter estimation with the help of R Language, I impasse because of the results I obtained are not similar with the estimation results in a reference paper I use. This is the following R Language syntax I use:

xi<-c(0.1,0.2,1,1,1,1,1,2,3,6,7,11,12,18,18,18,18,18,21,32,36,40,45,46,47,50,55,60,63,63,67,67,67,67,72,75,79,82,82,83,84,84,84,85,85,85,85,85,86,86);
n <-length (xi);
parameter <-function (theta, xi) {
logL<-(n*(log(theta[1])+log(theta[2])+log(theta[3])+(1-theta[2])*log(theta[4])+theta[1]*theta[4])+(theta[2]-1)*sum(log(xi))+(1/(theta[4]^theta[2]))*sum(xi^theta[2])-(theta[1]*theta[4])*sum(exp((xi/theta[4])^theta[2]))+(theta[3]-1)*sum(1-(exp((theta[1]*theta[4])*(1-(exp((xi/theta[4])^theta[2])))))))
return (-logL)
};
library(maxLik);
output <-maxLik (parameter, start = c (0.00007181,3.148,0.145,49.05), xi = xi);

Based on the syntax, the result of parameter estimation I get is:

theta [1] = 4.785855e-03
theta [2] = 1.759048e-04
theta [3] = 2.983679e + 04
theta [4] = 9.139192e + 02

While on paper belongs Sarhan and Apaloo (2013), the result should be as follows:

theta [1] = 2.506924e-06
theta [2] = 3.148000e + 00
theta [3] = 1.450000e-01
theta [4] = 4.905000e + 01

I am confused where my error in the above program. Previously, I apologize if I bothered you all. I really appreciate your help to me for completing this undergraduate thesis. Soon I will present this thesis and I found a lot of deadlocks. I really expect any help from you, no matter how small the help it would be greatly appreciated. Thank you so much

To be honest, I feel so sorry for my bad english grammar. I don't speak English

Your code doesn't work because it has `Return` and not `return`, but that makes me think you may have made other errors putting the code here. Please try and paste working code (and if you can do it without + signs so we can paste it, that's even better). I've just fixed these things in an edit, but I can't reproduce your results, so something else must be going on... — Spacedman, Jan 22 '16 at 16:30
Is your function supposed to be equation 13 in the cited paper? — Spacedman, Jan 22 '16 at 16:39
@Spacedman Yes in equation 13. I tried 'nlm' package and got the expected results. However, when trying to find the value of the log-likelihood ratio, the Kolmogorov-Smirnov test, AIC test, and the p-value, I got nothing. My lecturer then suggested using another package, but so far I have not succeeded — crhburn, Jan 22 '16 at 23:35

score 1 · Answer 1 · answered Jan 26 '16 at 04:20

Something was wrong with your likelihood function. I was unable to read it but note that maxLik maximizes the objective function, hence you have to return loglik, not -loglik. I re-wrote it in a more readable form (see Sarhan & Apaloo 2013) (sorry, but please name parameters, add space, split long equation into multiple lines...), and I also did not want to use name "parameter" for the log-likelihood function...

loglik <-function(theta, xi) {
   lambda <- theta[1]
   beta <- theta[2]
   gamma <- theta[3]
   alpha <- theta[4]
   xi.a <- xi/alpha
   A <- log(lambda) + log(beta) + log(gamma) + (beta - 1)*log(xi.a)
   LA1 <- lambda*alpha*(1 - exp(xi.a^beta))
   B <- xi.a^beta + LA1
   C <- log(1 - exp(LA1))
   logL <- A + B + (gamma - 1)*C
   return(logL)
}

library(maxLik);
start <- c(2.506924e-6, 3.148,0.145,49.05)
m <- maxLik(loglik, start=start, xi = xi);

It somewhat works. The main issue seems to be numerical instability. Play with different optimization methods, in particular BFGS seems to get you rather close:

summary(maxLik(loglik, start = start, method="bfgs", xi = xi))
--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 337 iterations
Return code 0: successful convergence 
Log-Likelihood: -213.3168 
4  free parameters
Estimates:
      Estimate Std. error t value  Pr(> t)    
[1,] 1.213e-05  5.976e-06   2.030   0.0423 *  
[2,] 3.133e+00  3.709e-02  84.456  < 2e-16 ***
[3,] 1.255e-01  1.897e-02   6.615 3.72e-11 ***
[4,] 4.496e+01         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Also, if you fix the first parameter, then you get the exact value with BHHH:

summary(maxLik(loglik, start = start, method="bhhh", xi = xi, fixed=1))
--------------------------------------------
Maximum Likelihood estimation
BHHH maximisation, 13 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -213.5116 
3  free parameters
Estimates:
      Estimate Std. error t value  Pr(> t)    
[1,] 2.507e-06  0.000e+00      NA       NA    
[2,] 3.091e+00  9.577e-03  322.74  < 2e-16 ***
[3,] 1.153e-01  2.234e-02    5.16 2.46e-07 ***
[4,] 4.189e+01  1.592e+00   26.32  < 2e-16 ***

It hints that the remaining problem is related to numerical instabilities originating from the first component (lambda). I can suggest two remedies:

supply analytic gradient to the maxLik function. I know, it is a hell to do, but maybe even supplying it for just lambda would suffice.
re-parameterize the problem. Even specifying lambda <- theta[1]/1e6, and adjusting the start value accordingly, seems to improve the convergence.

Note also that I removed summing from the likelihood function: now you can also use BHHH method which is often more robust than NR.

Error in maximum likelihood estimation using R

1 Answers1