[R] one step just of cliff-- zero hessian matrix in optim, with reproducable code and data
Hey Sky
heyskywalker at yahoo.com
Sun Sep 19 02:23:16 CEST 2010
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.
More information about the R-help
mailing list