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]]
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
Dear Arne, The elements of the theta vector are indeed strictly positive. I've just tried to use instead : lamda = log (theta), which means that theta = exp (lamda), so as to get rid of the log() function that appears in the log-likelihood and is causing the 50 warnings, but still the estimates I got for lamda and then those I got for theta (using theta=exp(lamda)) are irrelvant and their standard errors are infinite, which means that therer is still a problem that I can't yet figure out. Thanks, Maram On 18 July 2015 at 08:01, Arne Henningsen <arne.henningsen at gmail.com> wrote:> 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 >[[alternative HTML version deleted]]
Dear Arne, The elements of the theta vector are indeed strictly positive. I've just tried to use instead : lamda = log (theta), which means that theta = exp (lamda), so as to get rid of the log() function that appears in the log-likelihood and is causing the 50 warnings, but still the estimates I got for lamda and then those I got for theta (using theta=exp(lamda)) are irrelvant and their standard errors are infinite, which means that therer is still a problem that I can't yet figure out. Thanks, Maram> On 18 July 2015 at 08:01, Arne Henningsen <arne.henningsen at gmail.com> wrote: > 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[[alternative HTML version deleted]]