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

Maram Salem marammagdysalem at gmx.com
Tue Jul 7 22:29:56 CEST 2015


Dear Arne,

Sorry for posting my mail twice and thanks a lot for your help.

Best regards,
Maram 

Sent from my iPhone

> On Jul 7, 2015, at 9:55 PM, Arne Henningsen <arne.henningsen at gmail.com> wrote:
> 
> Dear Maram
> 
> Please do NOT post your message twice!
> 
> The warning messages occur each time, when maxLik() tries to calculate
> the logLik value for theta[1] <= 0, theta[1] + theta[2] <= 0, theta[3]
> <= 0 or something similar. According to the log-likelihood function,
> it seems that the parameters theta[1], theta[2], and theta[3] must be
> strictly positive. I suggest to re-parameterise your model so that the
> estimated parameters can take any values between minus infinity and
> infinity, e.g. by theta[1] <- exp( param[1] ); theta[2] <- exp(
> param[2] ); theta[3] <- exp( param[3] ) so that your estimated
> parameter vector 'param' consists of log( theta[1] ), log( theta[2] ),
> and log( theta[3] ). After the estimation, you can obtain the
> estimated values of the thetas by exp( param[1] ), exp( param[2] ),
> and exp( param[3] ) .
> 
> Best regards,
> Arne
> 
> 
> 
> 2015-07-06 2:29 GMT+02:00 Maram Salem <marammagdysalem at gmx.com>:
>> 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]]
>> 
>> ______________________________________________
>> 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



More information about the R-help mailing list