[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
Hey Sky
heyskywalker at yahoo.com
Sat Sep 18 17:36:49 CEST 2010
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: mydata9x.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20100918/5f7187d7/attachment.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: generate 0 hessian matrix.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20100918/5f7187d7/attachment-0001.txt>
More information about the R-help
mailing list