[R] Fwd: Warning message with maxLik()
marammagdysalem at gmail.com
marammagdysalem at gmail.com
Wed Jul 22 14:40:23 CEST 2015
Sent from my iPhone
Begin forwarded message:
> From: Maram SAlem <marammagdysalem at gmail.com>
> Date: July 21, 2015 at 11:40:56 PM GMT+2
> To: Arne Henningsen <arne.henningsen at gmail.com>
> Cc: "r-help at r-project.org" <r-help at r-project.org>
> Subject: Re: [R] Warning message with maxLik()
>
> 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]]
More information about the R-help
mailing list