[R] Warning message with maxLik()

[Arne Henningsen]

Dear Maram

- Please do not start a new thread for the same issue but reply to
previous messages in this thread [1].

- Please read my previous responses [1] more carefully, e.g. to use
"theta <- exp( param )" which guarantees that all elements of "theta"
are always positive.

[1] http://r.789695.n4.nabble.com/NaN-produced-from-log-with-positive-input-td4709463.html

Best regards,
Arne



2015-07-18 2:46 GMT+02:00 Maram SAlem <marammagdysalem at gmail.com>:
> Dear All,
> I'm trying to get the MLe for a certain distribution using maxLik ()
> function. I wrote the log-likelihood function as follows:
> theta <-vector(mode = "numeric", length = 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)
> # 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)
>  }
>
> then, I evaluated it at theta<- c(40,50,2)
>
> v<-loglik(param=theta)
> v
> [1] -56.66653
>
> I used this same log-likelihood function, once with analytic gradient and
> another time with numerical one, with the maxLik function, and in both
> cases I got the same 50 warning messages and an MLE which is completely
> unrealistic as per my applied example.
>
> a <- maxLik(loglik, gradlik, hesslik, start=c(40,50,2))
>
> where gradlik and hesslik are the analytic gradient and Hessian matrix,
> respectively, given by:
>
> 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)
> }
> 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)
> }
>
> and using numeric gradient and hessian matrix:
>
> a <- maxLik(loglik, start=c(40,50,2))
> 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
> 5: In log((C * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
> produced
> 6: In log(theta[3]) : NaNs produced
> 7: In log(theta[1] + theta[2]) : NaNs produced
> and so on…..
>
> I don't know why I get these 50 warnings although:
> 1- The inputs of the log() function are strictly positive.
> 2- When I evaluated the log-likelihood fuction at the very begining it gave
> me a number(which is -56.66) and not (NAN).
>
> I've also tried to:
> 1- Reparamtrize my model using lamda(i)= log(theta(i)), for i=1,2,3, so
> that it may solve the problem, but it didn't.
> 2- I've used the comparederivitive() function, and the analytic and numeric
> gradients were so close.
>
> Any help please?
> Maram Salem
>
>         [[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.



-- 
Arne Henningsen
http://www.arne-henningsen.name