I was able to get proper convergence in "BFGS", when I use the starting value from Nelder-Mead with 5000 iterations.
However, the hessian is not positive-definite. This indicates that you have a problem in your model. It seems to me that the model is over-parametrized. You have 20-odd parameters, but only 50 independent data points (I presume these are 50 time-series). In short, there is nothing wrong with optimization algorithms, but there is something not right with your model. Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: Hey Sky <heyskywal...@yahoo.com> Date: Saturday, September 18, 2010 11:38 am Subject: [R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data To: R <r-help@r-project.org> > Hey, R Users > > a few days ago I have met a zero hessian with optim command. I > reproduced it > with simulated data. plz check the code and data at the bottom of the > post. I > also attachment them with this email. hope it can reduce some > workload as > copying and pasting. > > I have simunated data many times and I do get convegence sometime and > hessian > matrix performs good. so, it would not be the problem of code lead to > this (I > may be wrong). > > > the error happens when the optim use a too large step to make some > values in the > optimization way too big and it never come back to normal again. the > values > before and after it happens as following: > > > values in the first part are reasonable. in the second part the W > value jumped > too large and lead to v8=Inf, which has been calculated from vbar2 > and vbar3. > and after that, even W come back to a little reasonable value (due to > simulated > value, I am not picky on it), the v8 is too large and lead to a zero > ccl > value all after that. > > > what may lead to this and any possible way to solve it? any > suggestion are > appreciated. > > **** values that jump ***** > w= 0.3157054 0.3678553 0.7879715 0.2859902 1.290479 > vbar2= -0.04085177 0.1922226 0.1922226 -0.04228498 0.1907894 > -0.0437182 > -0.2782258 -0.2782258 > > vbar3= -0.2226825 -0.2034284 -0.2034284 -0.06159623 -0.04234212 > 0.09949002 > 0.2413222 0.2413222 > > v8= 2.760340 3.027869 3.027869 2.898859 3.168746 3.061831 3.030057 > 3.030057 > lia= 0.289953 0.4002618 0.3302653 0.3243560 0.381919 0.3607669 > 0.4201014 > 0.3300268 > > wden= 1 0.2227371 1 1 0.297258 1 1 1 > lnw= 2.620900 2.1615 3.803036 3.533042 3.519460 2.328614 > regw= 1.260287 1.575992 1.575992 1.943847 2.259553 2.627408 2.995263 > 2.995263 > ccl= 1.546739e-05 1.134482e-05 4.232217e-06 0.0003085958 8.65926e-08 > > 6.858387e-07 1.572476e-05 > > --------------- > w= 71.7346 55.43801 55.13785 9.297906 -14.24756 > vbar2= -13725.15 -14549.52 -14549.52 -15240.92 -16065.29 -16756.70 > -17448.10 > -17448.10 > > vbar3= 15329.20 17870.51 17870.51 19927.61 22468.92 24526.02 26583.12 > 26583.12 > v8= Inf Inf Inf Inf Inf Inf Inf Inf > lia= NaN 0 0 NaN 0 NaN NaN 0 > wden= 1 -6.84084e-33 1 1 -3.222512e-96 1 1 1 > lnw= 2.620900 2.1615 3.803036 3.533042 3.519460 2.328614 > regw= 100.0510 171.7856 171.7856 227.2236 298.9582 354.3962 409.8342 > 409.8342 > ccl= NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN > --------------- > w= 14.59949 11.38189 11.65795 2.088373 -1.816329 > vbar2= -2559.668 -2586.704 -2586.704 -2619.031 -2646.067 -2678.394 > -2710.722 > -2710.722 > > vbar3= 2521.011 2623.554 2623.554 2722.231 2824.774 2923.451 3022.128 > 3022.128 > v8= Inf Inf Inf Inf Inf Inf Inf Inf > lia= NaN 0 0 NaN 0 NaN NaN 0 > wden= 1 -9.797435e-83 1 1 -2.080056e-247 1 1 1 > lnw= 2.620900 2.1615 3.803036 3.533042 3.519460 2.328614 > regw= 23.17728 37.77676 37.77676 49.15865 63.75813 75.14002 86.5219 > 86.5219 > ccl= NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN > --------------- > w= 3.172461 2.570661 2.961967 0.6464669 0.6699175 > vbar2= -503.5519 -503.2665 -503.2665 -505.674 -505.3886 -507.796 > -510.2035 > -510.2035 > > vbar3= 479.2945 483.5891 483.5891 490.9237 495.2183 502.553 509.8877 > 509.8877 > v8= 1.428720e+208 1.047267e+210 1.047267e+210 1.604982e+213 > 1.176469e+215 > 1.802990e+218 > > lia= 1 0 9.548661e-211 1 0 1 1 3.619044e-222 > wden= 1 4.817837e-20 1 1 6.500612e-71 1 1 1 > lnw= 2.620900 2.1615 3.803036 3.533042 3.519460 2.328614 > regw= 5.730039 8.9025 8.9025 11.47316 14.64562 17.21628 19.78695 > 19.78695 > ccl= 0 0 0 0 0 0 0 0 0 0 0 0 > > > > --------------------------------------------- > code and data also attached with this email > > #************** > # the main function > mymat<-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:10] > b<-par[11:19] > w<-par[20:npar] > > for(m in 1:ns){ > regw<-w[1]*acwrk+w[2]*actr+w[3]+w[4]*eta[m] > > vbar2=a[1]+ > > eta[m]+a[2]*acwrk+a[3]*actr+a[4]*edu+a[5]*v_refg+a[6]*v_econ+a[7]*age+regw*a[8] > > > vbar3=b[1]+b[2]*eta[m]+b[3]*acwrk+b[4]*actr+b[5]*edu+b[6]*v_refg+b[7]*v_econ+b[8]*age+regw*b[9] > > > 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] > } > > #**************************** > #cat("w=",w, "\n") > #cat("vbar2=",vbar2[1,], "\n") > #cat("vbar3=",vbar3[1,], "\n") > #cat("v8=",v8[1,], "\n") > #cat("lia=",lia[1,], "\n") > #cat("wden=",wden[1,], "\n") > #cat("lnw=",head(lnw), "\n") > #cat("regw=",regw[1,], "\n") > #cat("ccl=",ccl[1:6,], "\n") > #cat("---------------", "\n") > #**************************** > > func<-pr1*ccl[,1]+pr2*ccl[,2] > func<-ifelse(func<.Machine$double.xmin,0.00001,func) > f<-sum(log(func)) > return(-f) > } > > #********************************* > mydata<-read.table("F:/check the 0 hessian matrix > mistake/mydata9x.txt", head=F) > nt<<-8 # number of periods > ns<<-2 # number of person type > n<<-50 # number of people > npar<<-24 # number of parameters > > id<-as.numeric(mydata[,1]) > tr<-as.matrix(mydata[,2:(nt+1)]) > wrk<-as.matrix(mydata[,(nt+2):(2*nt+1)]) > home<-as.matrix(mydata[,(2*nt+2):(3*nt+1)]) > actr<-as.matrix(mydata[,(3*nt+2):(4*nt+1)]) > acwrk<-as.matrix(mydata[,(4*nt+2):(5*nt+1)]) > lnw<-as.numeric(mydata[,5*nt+2]) > edu<-as.numeric(mydata[,5*nt+3]) > age<-as.numeric(mydata[,5*nt+4]) > v_refg<-as.numeric(mydata[,5*nt+5]) > v_econ<-as.numeric(mydata[,5*nt+6]) > > # the initial guess > guess<-rep(0.5,times=npar) > guess[npar]<-1.0 > > # use "Nelder-Mead" to get the initial value > system.time(r1<-optim(guess,mymat,data=mydata, hessian=F)) > guess<-r1$par > > system.time(r2<-optim(guess,mymat,data=mydata, > method="BFGS",hessian=T, > control=list(trace=T, maxit=1000))) > > std.err<-sqrt(diag(solve(r2$hessian))) > res<-cbind(r2$par,std.err,r2$par/std.err) > colnames(res)<-c("parameter","std.err","t test") > > > > ------------------------------------------------- > the data > "1" 1 0 0 1 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 1 1 1 1 2 2 3 4 4 0 > 1 1 1 2 2 > 2 2 2.62089951476082 16 29 0 0 > "2" 1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 1 0 1 2 2 3 3 3 3 3 0 > 0 0 0 1 1 > 1 2 2.16150014568120 4 19 1 0 > "3" 1 1 1 0 1 0 1 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 2 3 3 4 4 5 6 0 > 0 0 1 1 2 > 2 2 3.80303575377911 16 26 1 0 > "4" 1 0 0 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1 2 3 4 5 6 0 > 1 1 1 1 1 > 1 1 3.53304197313264 16 41 0 1 > "5" 0 0 0 0 1 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 2 3 1 > 2 3 3 3 3 > 3 3 3.51945951068774 3 35 0 0 > "6" 0 0 0 0 1 0 0 0 1 1 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 1 1 1 1 1 > 2 2 3 3 3 > 4 5 2.32861361233518 17 22 0 1 > "7" 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 1 2 3 3 3 4 4 4 0 > 0 0 0 0 0 > 0 1 2.89729301305488 14 26 0 1 > "8" 0 0 1 1 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 2 2 3 3 3 1 > 2 2 2 2 2 > 2 2 2.86090020649135 4 22 0 0 > "9" 0 1 0 1 1 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 1 2 3 4 4 5 0 > 0 1 1 1 1 > 2 2 2.59020843589678 17 23 0 0 > "10" 0 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 2 3 4 4 5 5 > 1 1 1 1 1 1 > 1 1 3.6295328931883 5 22 0 0 > "11" 0 0 0 1 0 0 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1 2 3 > 1 2 3 3 4 4 > 4 4 2.02498448966071 11 26 0 1 > "12" 1 1 1 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1 2 3 3 3 4 4 5 > 0 0 0 1 2 2 > 3 3 3.25450395001099 13 31 0 1 > "13" 1 0 0 1 1 0 0 1 0 0 0 0 0 1 1 0 0 1 1 0 0 0 0 0 1 1 1 2 3 3 3 4 > 0 0 0 0 0 1 > 2 2 2.37046055402607 14 33 0 1 > "14" 0 0 0 0 0 1 1 1 1 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 2 3 > 1 2 2 3 3 3 > 3 3 2.87286716327071 9 23 1 0 > "15" 1 1 0 1 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 1 2 2 3 4 5 5 5 > 0 0 1 1 1 1 > 2 3 2.90179902175441 15 36 0 1 > "16" 1 1 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 0 1 2 2 2 2 2 2 3 > 0 0 1 2 3 4 > 4 4 3.32979543972760 6 25 0 0 > "17" 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1 2 > 1 1 1 2 2 3 > 3 3 2.36153599619865 17 45 0 1 > "18" 1 0 1 0 0 0 1 0 0 1 0 1 0 1 0 1 0 0 0 0 1 0 0 0 1 1 2 2 2 2 3 3 > 0 1 1 2 2 3 > 3 4 3.63236659532413 3 41 1 0 > "19" 1 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 1 2 2 2 3 4 5 5 > 0 0 1 2 2 2 > 2 3 3.69187993369997 2 41 0 0 > "20" 0 1 1 0 0 1 1 0 1 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 1 2 2 2 3 4 4 > 1 1 1 2 3 3 > 3 4 2.01738612353802 8 33 0 0 > "21" 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 3 4 5 6 7 8 > 0 0 0 0 0 0 > 0 0 3.50919563509524 13 22 0 0 > "22" 1 0 0 0 1 1 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 1 1 2 3 4 4 > 0 1 2 2 2 2 > 2 3 3.14363623457029 5 33 0 1 > "23" 1 0 0 0 0 0 1 1 0 1 0 1 1 1 0 0 0 0 1 0 0 0 0 0 1 1 1 1 1 1 2 3 > 0 1 1 2 3 4 > 4 4 2.78580305865034 11 19 1 0 > "24" 1 1 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 2 2 3 4 4 5 6 > 0 0 0 0 0 0 > 0 0 3.91743207862601 9 40 0 1 > "25" 1 1 1 0 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 2 3 3 4 5 5 6 > 0 0 0 1 1 1 > 1 1 3.63302375609055 16 33 0 1 > "26" 1 1 1 0 1 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 2 3 3 4 5 6 7 > 0 0 0 1 1 1 > 1 1 2.28801752673462 3 32 0 1 > "27" 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 0 0 0 1 1 2 2 3 4 > 0 0 0 0 0 1 > 1 1 2.45849566301331 12 18 1 0 > "28" 0 0 0 1 0 0 1 1 1 0 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 1 1 2 3 > 1 1 2 2 2 2 > 2 2 2.74557595746592 8 42 0 0 > "29" 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 1 1 1 2 2 2 2 3 3 > 0 0 0 1 1 2 > 2 2 2.00150080351159 15 32 1 0 > "30" 1 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 1 1 1 1 2 2 2 2 2 2 > 0 1 1 2 2 3 > 3 3 2.72582565387711 14 19 0 1 > "31" 0 0 1 0 1 0 0 1 0 0 0 1 0 1 1 0 1 1 0 0 0 0 0 0 0 0 1 1 2 2 2 3 > 0 0 0 1 1 2 > 3 3 2.88708175066859 10 34 0 0 > "32" 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 2 3 3 3 4 5 6 > 0 0 0 0 0 0 > 0 0 2.24319696752355 6 39 0 0 > "33" 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 1 0 1 1 2 2 2 3 3 3 > 0 0 0 0 1 1 > 1 2 2.6321357563138 16 35 0 1 > "34" 0 0 1 0 1 0 1 1 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 1 2 2 3 4 > 0 1 1 2 2 3 > 3 3 3.26070732064545 17 28 0 1 > "35" 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 1 1 0 1 0 1 0 1 0 0 1 1 1 1 1 1 > 0 0 0 0 1 1 > 2 2 3.4693668698892 7 39 0 0 > "36" 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1 1 2 3 4 5 5 > 0 0 0 0 0 0 > 0 0 2.60646418808028 10 22 0 1 > "37" 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 > 0 0 1 2 2 2 > 2 3 3.45602289121598 15 28 0 1 > "38" 1 1 1 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1 2 3 4 5 5 5 5 > 0 0 0 0 0 1 > 1 2 3.03841971792281 6 41 0 0 > "39" 1 0 1 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 1 1 2 3 3 4 4 5 > 0 0 0 0 1 1 > 1 1 2.90754798706621 4 36 0 0 > "40" 1 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 2 2 2 2 3 4 4 > 0 0 0 1 1 1 > 1 2 3.69572683610022 11 19 0 1 > "41" 1 0 1 0 1 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 2 2 3 4 4 5 > 0 1 1 2 2 2 > 2 2 2.81628963444382 2 36 0 1 > "42" 1 0 0 1 1 0 1 1 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 1 1 2 3 3 4 5 > 0 0 0 0 0 1 > 1 1 2.94380070734769 17 31 0 1 > "43" 0 0 1 0 1 0 1 1 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 1 2 2 3 4 > 0 1 1 1 1 2 > 2 2 2.50514903757721 8 38 0 0 > "44" 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 1 0 0 0 1 1 0 0 0 1 1 1 1 1 2 2 > 0 0 1 2 2 2 > 2 3 3.39924295153469 3 19 0 0 > "45" 1 0 0 0 1 1 0 0 0 1 1 1 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 1 2 3 3 3 > 0 1 2 3 3 3 > 4 5 2.29968624887988 8 32 0 1 > "46" 0 0 0 0 1 1 1 1 0 1 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 2 3 4 > 0 1 2 2 2 2 > 2 2 2.58306567557156 15 27 0 1 > "47" 1 1 1 0 1 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 2 3 3 4 4 5 6 > 0 0 0 0 0 1 > 1 1 3.99967893399298 6 42 0 1 > "48" 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1 1 1 2 3 3 > 0 0 0 0 1 1 > 1 1 3.6599674411118 10 21 0 0 > "49" 0 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 1 2 2 3 3 3 3 > 1 1 1 1 1 1 > 1 1 2.35007652500644 1 30 0 0 > "50" 1 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 1 2 2 3 4 4 5 5 > 0 0 0 0 0 0 > 0 0 2.07408210681751 9 38 1 0 > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list 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.