[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
Ravi Varadhan
rvaradhan at jhmi.edu
Sat Sep 18 22:03:41 CEST 2010
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.
More information about the R-help
mailing list