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
"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
#*************** #the 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("j:/new folder/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]) # initial value guess<-rep(0.1,times=npar) guess[npar]<-1.0 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") #res
______________________________________________ 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.