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