# Initial guess for Newton-Raphson iteration method in Maximum Likelihood Estimation

I want to estimate the four parameters of Exponentiated Modified Weibull Extension (EMWE) distribution introduced by Sarhan and Apaloo (2013) with the Maximum Likelihood Estimation (MLE). This distribution is used in reliability and survival analysis to analyze the dataset with bathtub hazard function.

Because the first derivative of the log-likelihood function of the four parameters give implicit solutions then I tried to continue with the Newton-Raphson iteration method. It's just that I am not familiar in determining the good initial guess for this method because my main focus is how to choose the right initial guess to obtain the bathtub hazard function. I'm using RStudio and here is the code:

#Parameter Estimation
xi<-c(3,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,9,9,9,9,9,10,10,11,12,12,14)
n<-length(xi)
log_likelihood<-function(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]))))))
}

result<-nlm(log_likelihood,theta<-c(0.000030763,2.2011,0.11878,6),hessian=TRUE)


I realized that in any optimization problem, there is no systematic choice of an initial guess. In Maximum Likelihood problems, though, a common choice is to choose initial parameters that match the empirical moments (mean, variance, etc...) of the observed data, so I choose:

theta[1]=0.0000022011
theta[2]=1.1878
theta[3]=0.21669
theta[4]=6.4941


where theta[1] and theta[4] are the scale parameter, theta[2] and theta[3] as shape parameter, and here is the result:

result
____________________
$minimum [1] -1167.797$estimate
[1] 1.540561e-06 1.187800e+00 2.166900e-01 6.494100e+00

$gradient [1] 4.251830e+07 7.312058e+01 3.922662e+02 -1.815774e+01$hessian
[,1]          [,2]          [,3]          [,4]
[1,] -2.977361e+10 -7.455449e+02  1.093871e+03  3.396184e+02
[2,] -7.455449e+02 -5.215852e+01  6.316999e-04 -2.725559e+01
[3,]  1.093871e+03  6.316999e-04 -1.808591e+03 -2.906025e-04
[4,]  3.396184e+02 -2.725559e+01 -2.906025e-04  5.665718e+00

$code [1] 2$iterations
[1] 1


The estimated parameters are:

theta[1]=1.540561e-06
theta[2]=1.1878
theta[3]=2.166900e-01
theta[4]=6.4941


To obtain the hazard function plot, I use this code:

#The estimated parameters
p1=theta[4]
p2=theta[2]
p3=theta[3]
p4=theta[1]

#Hazard Function
hazard_EMWE=(p4*p2*p3*((xi/p1)^(p2-1))*(exp(((xi/p1)^(p2))+(p4*p1*(1-(exp((xi/p1)^(p2)))))))*((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3-1)))/(1-(((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3))))

#Cumulative Distribution Function
cum_EMWE=((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3))

#Hazard Function Plot
plot(xi,hazard_EMWE,type="l" ,main="Hazard Function Plot",xlab="x",ylab="Hazard Function",col = "black",lty=1,lwd=3,ylim=c(0.005,0.008))


Here is the plot:

Actually, the plot exhibits the bathtub hazard curve, but I thought that I had chosen the bad initial guess so that the cumulative distribution function is badly different with the empirical distribution function. It means that the data doesn't follow the EMWE distribution whereas my expectation is the data have to follow the EMWE distribution.

So my problem here is related to the initial guess that I use. It could be the poor initial value. Does anyone here has solution in choosing the good initial guess for this dataset?

• Use E-M. Just a few E-M iterations are guaranteed to move the candidate solution towards a (local?) minimum. Then you can pass that suboptimal solution to the N-R to reach the (local?) minimum faster. – usεr11852 Mar 3 '16 at 6:46
• @usεr11852 How to construct EM iteration in R? – crhburn Mar 4 '16 at 12:39
• It will be particular to your problem. See here for an introduction. – usεr11852 Mar 4 '16 at 17:45