Hey Sky
2010-Sep-18 15:36 UTC
[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
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>
Ravi Varadhan
2010-Sep-18 20:03 UTC
[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
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: rvaradhan at jhmi.edu ----- Original Message ----- From: Hey Sky <heyskywalker at 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 at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.
Hey Sky
2010-Sep-19 00:23 UTC
[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
Dear Ravi thanks for your reply. I think the number of obs may do be the problem for a zero hessian matrix. since it is simulated data, I have increased the sample size to 500 obs and tried for around 20 times, the zero hessian did not appear, compared with the fact that it happens too often with only 50 obs. my real data is a confendential one so I only?be able to try?the simulated?one for a?discussion here. but what?may?lead to the NaN in std.err, or the negative variance when I inverse the hessian matrix? the R code is translated from a Fortran code, which has been used for a long time,?thus the model should be ok. besides a more reasonable initial value, whatelse do you think it might be the reasons and what do you suggest? thanks for your time Nan ? ----- Original Message ---- From: Ravi Varadhan <rvaradhan at jhmi.edu> To: Hey Sky <heyskywalker at yahoo.com> Cc: R <r-help at r-project.org> Sent: Sat, September 18, 2010 4:03:41 PM Subject: Re: [R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data 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: rvaradhan at jhmi.edu ----- Original Message ----- From: Hey Sky <heyskywalker at 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 at 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 at r-project.org mailing list >? >? PLEASE do read the posting guide >? and provide commented, minimal, self-contained, reproducible code.