[R] NaN produced from log() with positive input

[Maram Salem]

Dear All
I'm trying to find the maximum likelihood estimator  of a certain distribution based on the newton raphson method using maxLik package. I wrote the log-likelihood , gradient, and hessian functionsusing the following code.
 
#Step 1: Creating the theta vector
 theta <-vector(mode = "numeric", length = 3)
# Step 2: Setting the values of r and n
r<- 17
n <-30
 # Step 3: Creating the T vector
T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
# Step 4: Creating the C vector
C<- c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
# The  loglik. func.
loglik <- function(param) {
 theta[1]<- param[1]
 theta[2]<- param[2]
 theta[3]<- param[3]
 l<-(r*log(theta[3]))+(r*log(theta[1]+theta[2]))+(n*theta[3]*log(theta[1]))+(n*theta[3]*log(theta[2]))+ (-1*(theta[3]+1))*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+ (-1*theta[3]*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2]))))
return(l)
 }
# Step 5: Creating the gradient vector and calculating its inputs
U <- vector(mode="numeric",length=3)
gradlik<-function(param = theta,n, T,C)
 {
U <- vector(mode="numeric",length=3)
theta[1] <- param[1]
theta[2] <- param[2]
theta[3] <- param[3]
r<- 17
n <-30
T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
C<- c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
 U[1]<- (r/(theta[1]+theta[2]))+((n*theta[3])/theta[1])+( -1*(theta[3]+1))*sum((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+ (-1*(theta[3]))*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
U[2]<-(r/(theta[1]+theta[2]))+((n*theta[3])/theta[2])+ (-1*(theta[3]+1))*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+ (-1*(theta[3]))*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
U[3]<-(r/theta[3])+(n*log(theta[1]*theta[2]))+ (-1)*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+(-1)*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2])))
return(U)
}
# Step 6: Creating the G (Hessian) matrix and Calculating its inputs
hesslik<-function(param=theta,n,T,C)
{
theta[1] <- param[1]
theta[2] <- param[2]
theta[3] <- param[3]
r<- 17
n <-30
T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
C<- c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
G<- matrix(nrow=3,ncol=3)
G[1,1]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[1])^2)+ (theta[3]+1)*sum(((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+( theta[3])*sum(((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[1,2]<-((-1*r)/((theta[1]+theta[2])^2))+ (theta[3]+1)*sum(((T)/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+ (theta[3])*sum(((C)/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,1]<-G[1,2]
G[1,3]<-(n/theta[1])+(-1)*sum( (T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,1]<-G[1,3] 
G[2,2]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[2])^2)+ (theta[3]+1)*sum(((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+( theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]<-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]<-G[2,3]                 
G[3,3]<-((-1*r)/(theta[3])^2)
return(G)
}
mle<-maxLik(loglik, grad = gradlik, hess = hesslik, start=c(40,50,2))
There were 50 or more warnings (use warnings() to see the first 50)
 
warnings ()
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs produced
 and so on .......
 
Although when I evaluate, for example, log(theta[3])  it gives me a number. and the same applies for the other warnings.
 
Then when I used summary (mle), I got
 
 
Maximum Likelihood estimation
Newton-Raphson maximisation, 7 iterations
Return code 1: gradient close to zero
Log-Likelihood: -55.89012 
3  free parameters
Estimates:
     Estimate Std. error t value Pr(> t)
[1,]   11.132        Inf       0       1
[2,]   47.618        Inf       0       1
[3,]    1.293        Inf       0       1
--------------------------------------------


Where the estimates are far away from the starting values and they have infinite standard errors. I think there is a problem with my gradlik or hesslik functions, but I can't figure it out.
Any help?
Thank you in advance.

Maram 



	[[alternative HTML version deleted]]