thr3ads.net - R help - [R] NaN produced from log() with positive input [Jul 2015]

If this information is useful, please help other people find it:
Share via:

Maram Salem

2015-Jul-06 00:29 UTC

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

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]]

Maram Salem

2015-Jul-07 12:46 UTC

head link

[R] Can't get a reasonable the maximum likelihood estimator

Sent from my iPhone

Begin forwarded message:
> From: "Maram Salem" <marammagdysalem at gmx.com>
> Date: July 6, 2015 at 2:29:56 AM GMT+2
> To: "r-help at r-project.org" <r-help at r-project.org>
> Subject: [R] NaN produced from log() with positive input
> 
> 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 functions using 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.
	[[alternative HTML version deleted]]

Arne Henningsen

2015-Jul-07 19:55 UTC

head link

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

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

Maram Salem

2015-Jul-07 20:29 UTC

head link

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

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

Arne Henningsen

2015-Jul-12 06:23 UTC

head link

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

Dear Maram

On 8 July 2015 at 17:52, Maram Salem <marammagdysalem at gmx.com>
wrote:> Dear Arne,
>
> On a second thought, as per your mail "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."
>
> The component of the theta vector are all indeed strictly positive, and the
> initial values I used are c(40,50,2). and this means that " theta[1]
> 0,
> theta[1] + theta[2] > 0, and theta[3] > 0 ".  These initial
values are the
> parameter values that were used for generating the data (the C and T
> vectors). That's why I don't know why the warnings occur in the
first place
> and why the estimates are far away from the initial values.a
> Any suggestions?
- don't send your messages twice;

- do what I suggested in my previous e-mail;

- increase the number of observations;

- check your log-likelihood function;

- use an optimisation algorithm that does not use derivatives;

- use numeric derivatives in the optimisation;

- check your function for returning the derivatives of the
log-likelihood function, e.g. with compareDerivatives();

- check the data generating process

Best regards,
Arne

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

Maram Salem

2015-Jul-12 14:52 UTC

head link

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

Dear Arne,
Will do. Thanks for helping.
Maram

Sent from my iPhone
> On Jul 12, 2015, at 8:23 AM, Arne Henningsen <arne.henningsen at
gmail.com> wrote:
> 
> Dear Maram
> 
>> On 8 July 2015 at 17:52, Maram Salem <marammagdysalem at gmx.com>
wrote:
>> Dear Arne,
>> 
>> On a second thought, as per your mail "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."
>> 
>> The component of the theta vector are all indeed strictly positive, and
the
>> initial values I used are c(40,50,2). and this means that "
theta[1] > 0,
>> theta[1] + theta[2] > 0, and theta[3] > 0 ".  These initial
values are the
>> parameter values that were used for generating the data (the C and T
>> vectors). That's why I don't know why the warnings occur in the
first place
>> and why the estimates are far away from the initial values.a
>> Any suggestions?
> 
> - don't send your messages twice;
> 
> - do what I suggested in my previous e-mail;
> 
> - increase the number of observations;
> 
> - check your log-likelihood function;
> 
> - use an optimisation algorithm that does not use derivatives;
> 
> - use numeric derivatives in the optimisation;
> 
> - check your function for returning the derivatives of the
> log-likelihood function, e.g. with compareDerivatives();
> 
> - check the data generating process
> 
> Best regards,
> Arne
> 
> -- 
> Arne Henningsen
> http://www.arne-henningsen.name

R help - Jul 2015 - NaN produced from log() with positive input

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

[R] Can't get a reasonable the maximum likelihood estimator

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

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

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

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