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@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.
