[R] NaN produced from log() with positive input
Arne Henningsen
arne.henningsen at gmail.com
Tue Jul 7 21:55:08 CEST 2015
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