Dear R-helper,
I am trying to do maximum likelihood estimation in R. I use the
"optim" function. Since I have no prior information on the true
values of the parameters, I just randomly select different sets of starting
values to feed into the program. Each time, I get the following error message:
Error in optim(theta0, lf, method = "BFGS", hessian = T, Y = Y, X = X,
: non-finite finite-difference value [6]
In addition: There were 50 or more warnings (use warnings() to see the first 50)
I guess this has something to do with the starting values that I chose. But I
have no idea as to how to pick better starting values in a case like mine or how
to get rid of this error message. Could you give me some suggestion? Thanks a
lot for your kind help.
My R code is as follows. If more information is needed, please let me know.
lf<-function(theta,Y,X,G) {
N<-26*26
l<-matrix(1,26,1)
alpha<-theta[1:4]
betad1<-t(t(theta[5:6]))
betad2<-t(t(theta[7:8]))
betad3<-t(t(theta[9:10]))
betad4<-t(t(theta[11:12]))
betao1<-t(t(theta[13:14]))
betao2<-t(t(theta[15:16]))
betao3<-t(t(theta[17:18]))
betao4<-t(t(theta[19:20]))
gamma<-theta[21:24]
pd<-theta[25]
po<-theta[26]
pw<-theta[27]
E1<- Y-alpha[1]*l%*%t(l)-X%*%betad1%*%t(l)-l%*%t(betao1)%*%t(X)-gamma[1]*G
E2<- WY-alpha[2]*l%*%t(l)-X%*%betad2%*%t(l)-l%*%t(betao2)%*%t(X)-gamma[2]*G
E3<- YW-alpha[3]*l%*%t(l)-X%*%betad3%*%t(l)-l%*%t(betao3)%*%t(X)-gamma[3]*G
E4<-WYW-alpha[4]*l%*%t(l)-X%*%betad4%*%t(l)-l%*%t(betao4)%*%t(X)-gamma[4]*G
Q<-matrix(0,4,4)
Q[1,1]<-sum(E1*t(E1))
Q[1,2]<-sum(E1*t(E2))
Q[1,3]<-sum(E1*t(E3))
Q[1,4]<-sum(E1*t(E4))
Q[2,1]<-sum(E2*t(E1))
Q[2,2]<-sum(E2*t(E2))
Q[2,3]<-sum(E2*t(E3))
Q[2,4]<-sum(E2*t(E4))
Q[3,1]<-sum(E3*t(E1))
Q[3,2]<-sum(E3*t(E2))
Q[3,3]<-sum(E3*t(E3))
Q[3,4]<-sum(E3*t(E4))
Q[4,1]<-sum(E4*t(E1))
Q[4,2]<-sum(E4*t(E2))
Q[4,3]<-sum(E4*t(E3))
Q[4,4]<-sum(E4*t(E4))
T<-matrix(c(1, -pd, -po, -pw),nrow=4)
logl<- log(det(diag(1,N)-pd*wd - po*wo -pw*ww))-N/2*log(t(T)%*%Q%*%T)
return(-logl)
}
theta0<-c(0,0,0,0, 0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0, 0,0,0,0, 0.2, 0.1, 0)
p<-optim(theta0,lf,method="BFGS", hessian=T, Y=Y, X=X, G=G)
OI<-solve(p$hessian)
Hotmail is redefining busy with tools for the New Busy. Get more from your
inbox. See how.
_________________________________________________________________
Hotmail is redefining busy with tools for the New Busy. Get more from your
inbox.
N:WL:en-US:WM_HMP:042010_2
[[alternative HTML version deleted]]
