[R] optim with BFGS--what may lead to this, a strange thing happened

[Hey Sky]

Dear R Users

on a self-written function for calculating maximum likelihood probability (plz 
check function code at the bottom of this message), one value, wden, suddenly 
jump to zero. detail info as following:


w[11]=2.14
lnw            =2.37 2.90 3.76 ...
regw            =1.96 1.77 1.82 ....
wden=0.182 0.178 0.179...

w[11]=2.14
lnw=2.37 2.90 3.76 ...
regw =1.96 1.77 1.82 ....
wden=0.182 0.179 0.180...

w[11]=-0.51
lnw=2.37 2.90 3.76 ...
regw=57.92 57.50 57.60 ...
wden=0 0 0 0 ....1 0 0 0

which I use cat() to get (by the way, how to write from the result from cat() 
directly to a file?). 


the lnw is coming from data. and w[11] is one of the estimated parameters; wden 
and regw is estimated equation. wrk is a 0/1 dummy

wden equation:    wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)

when I use simulated data, it works fine. but when using real data, it does not 
work. when I did not add the following into the function,

func<-ifelse(func<.Machine$double.xmin,0.000001,func)

it will iterate for a while and give a non-finite finite difference in  par[54] 
report. but when I add it, very soon it will give the upper  result. what may 
lead to this and how to fix it?

thanks for any suggestion in advance.

Nan
from Montreal



# the main function
mymatrix<-function(par,data) {

# define the parameter matrix used in following part
vbar2<-matrix(0,n,nt)
vbar3<-matrix(0,n,nt)
v8 <-matrix(0,n,nt)
regw<-matrix(0,n,nt)
wden<-matrix(0,n,nt)
lia<-matrix(0,n,nt)
ccl<-matrix(1,n,ns)
eta<-c(0,0)

# setup the parts for loglikelihood
q1<-exp(par[1])
pr1<-q1/(1+q1)
pr2<-1-pr1

eta[2]<-par[2]

a<-par[3:6]
b<-par[7:11]
w<-par[12:npar]

for(m in 1:ns){

    regw<-exp(w[1]+w[2]*eta[m])*actr+exp(w[3]+ w[4]*eta[m])*acwrk
    
    vbar2=a[1]+     eta[m]+regw*a[2]+acwrk*a[3]+actr*a[4]
     vbar3=b[1]+b[2]*eta[m]+regw*b[3]+acwrk*b[4]+actr*b[5]

    v8=1+exp(vbar2)+exp(vbar3)

     lia<-ifelse(home==1,1/v8,
         ifelse(wrk==1,exp(vbar2)/v8,
            ifelse(tr==1,exp(vbar3)/v8,1)))

    wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)

    ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]*lia[,6]*lia[,7]*lia[,8]*
        wden[,1]*wden[,2]*wden[,3]*wden[,4]*wden[,5]*wden[,6]*wden[,7]*wden[,8]

}

func<-pr1*ccl[,1]+pr2*ccl[,2]
#
func<-ifelse(func<.Machine$double.xmin,0.000001,func)
#
f<-sum(log(func))
return(-f)
}


#*******************************************
mydata<-read.table("C:/Documents and Settings/lciqss/Desktop/R-2.11.1/my 
data_wage_eq.txt", head=F)

nt<<-8    # number of periods
ns<<-2
n<<-50    # number of people
npar<<-16    # number of parameters

tr<-as.matrix(mydata[,1:nt])
wrk<-as.matrix(mydata[,(nt+1):(2*nt)])
home<-as.matrix(mydata[,(2*nt+1):(3*nt)])
actr<-as.matrix(mydata[,(3*nt+1):(4*nt)])
acwrk<-as.matrix(mydata[,(4*nt+1):(5*nt)])
lnw<-as.numeric(mydata[,(5*nt+1)])

guess<-rep(0,times=npar)
guess[npar]<-1.0

system.time(r1<-optim(guess,mymatrix,data=mydata, hessian=F))
guess<-r1$par


system.time(rmat<-optim(guess,mymatrix,data=mydata, method="BFGS",hessian=T, 
        control=list(trace=T, maxit=1000)))

rmat$par
std.err<-sqrt(diag(solve(rmat$hessian)))
result<-cbind(rmat$par,std.err,rmat$par/std.err)
colnames(result)<-c("paramters","std,err","t test")
#print(result)