

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

