[R] Can't get a reasonable the maximum likelihood estimator

[Maram Salem]

Sent from my iPhone

Begin forwarded message:

> From: "Maram Salem" <marammagdysalem at gmx.com>
> Date: July 6, 2015 at 2:29:56 AM GMT+2
> To: "r-help at r-project.org" <r-help at r-project.org>
> Subject: [R] NaN produced from log() with positive input
> 
> 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 functions using 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]]
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

	[[alternative HTML version deleted]]