[R] Hessian matrix issue

[John C Nash]

Unless you are supplying analytic hessian code, you are almost certainly getting an
approximation. Worse, if you do not provide gradients, these are the result of two levels
of differencing, so you should expect some loss of precision in the approximate Hessian.

Moreover, if your estimate of the optimum is a little bit off, or the optimizer has
terminated (algorithms converge, programs terminate) to a point that is not an optimum,
there is no reason the Hessian should be positive definite.

Package optimx() uses the Jacobian of the gradient if the analytic gradient is available.
This drops the differencing to 1 level. Even better is to code the Hessian, but that is
messy and tedious in most cases.

Best, JN


On 09/03/2011 06:00 AM, r-help-request at r-project.org wrote:
> Message: 59
> Date: Fri, 2 Sep 2011 15:33:13 -0400
> From: tzaihra at alcor.concordia.ca
> To: r-help at r-project.org
> Subject: [R] Hessian Matrix Issue
> Message-ID:
> 	<e6dc43b4487eb4a4055e1ab485f015f0.squirrel at webmail.concordia.ca>
> Content-Type: text/plain;charset=iso-8859-1
> 
> Dear All,
> 
> I am running a simulation to obtain coverage probability of Wald type
> confidence intervals for my parameter d in a function of two parameters
> (mu,d).
> 
> I am optimizing it using "optim" method "L-BFGS-B" to obtain MLE. As, I
> want to invert the Hessian matrix to get Standard errors of the two
> parameter estimates. However, my Hessian matrix at times becomes
> non-invertible that is it is no more positive definite and I get the
> following error msg:
> 
> "Error in solve.default(ac$hessian) : system is computationally singular:
> reciprocal condition number = 6.89585e-21"
> Thank you
> 
> Following is the code I am running I would really appreciate your comments
> and suggestions:
> 
> #Start Code
> #option to trace /recover error
> #options(error = recover)
> 
> #Sample Size
> n<-30
> mu<-5
> size<- 2
> 
> #true values of parameter d
> d.true<-1+mu/size
> d.true
> 
> #true value of  zero inflation index phi= 1+log(d)/(1-d)
> z.true<-1+(log(d.true)/(1-d.true))
> z.true
> 
> # Allocating space for simulation vectors and setting counters for simulation
> counter<-0
> iter<-10000
> lower.d<-numeric(iter)
> upper.d<-numeric(iter)
> 
> #set.seed(987654321)
> 
> #begining of simulation loop########
> 
> for (i in 1:iter){
> r.NB<-rnbinom(n, mu = mu, size = size)
> y<-sort(r.NB)
> iter.num<-i
> print(y)
> print(iter.num)
> #empirical estimates or sample moments
> xbar<-mean(y)
> variance<-(sum((y-xbar)2))/length(y)
> dbar<-variance/xbar
> #sample estimate of proportion of zeros and zero inflation index
> pbar<-length(y[y==0])/length(y)
> 
> ### Simplified function #############################################
> 
> NegBin<-function(th){
> mu<-th[1]
> d<-th[2]
> n<-length(y)
> 
> arg1<-n*mean(y)*ifelse(mu >= 0, log(mu),0)
> #arg1<-n*mean(y)*log(mu)
> 
> #arg2<-n*log(d)*((mean(y))+mu/(d-1))
> arg2<-n*ifelse(d>=0, log(d), 0)*((mean(y))+mu/ifelse((d-1)>= 0, (d-1),
> 0.0000001))
> 
> aa<-numeric(length(max(y)))
> a<-numeric(length(y))
> for (i in 1:n)
> {
> for (j in 1:y[i]){
> aa[j]<-ifelse(((j-1)*(d-1))/mu >0,log(1+((j-1)*(d-1))/mu),0)
> #aa[j]<-log(1+((j-1)*(d-1))/mu)
> #print(aa[j])
> }
> 
> a[i]<-sum(aa)
> #print(a[i])
> }
> a
> arg3<-sum(a)
> llh<-arg1+arg2+arg3
> if(! is.finite(llh))
> llh<-1e+20
> -llh
> }
> ac<-optim(NegBin,par=c(xbar,dbar),method="L-BFGS-B",hessian=TRUE,lower=
> c(0,1) )
> ac
> print(ac$hessian)
> muhat<-ac$par[1]
> dhat<-ac$par[2]
> zhat<- 1+(log(dhat)/(1-dhat))
> infor<-solve(ac$hessian)
> var.dhat<-infor[2,2]
> se.dhat<-sqrt(var.dhat)
> var.muhat<-infor[1,1]
> se.muhat<-sqrt(var.muhat)
> var.func<-dhat*muhat
> var.func
> d.prime<-cbind(dhat,muhat)
> 
> se.var.func<-d.prime%*%infor%*%t(d.prime)
> se.var.func
> lower.d[i]<-dhat-1.96*se.dhat
> upper.d[i]<-dhat+1.96*se.dhat
> 
> if(lower.d[i]  <= d.true & d.true<= upper.d[i])
> counter <-counter+1
> }
> counter
> covg.prob<-counter/iter
> covg.prob
> 
> 
>