[R] Hessian matrix issue

[Giovanni Petris]

About the numerical calculation of the Hessian matrix, I have found
numDeriv:::hessian to be often more accurate than the Hessian returned
by optim.

Best,
Giovanni Petris

On Sat, 2011-09-03 at 13:39 -0400, John C Nash wrote:
> 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
> > 
> > 
> >
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.