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-requ...@r-project.org wrote: > > Message: 59 > > Date: Fri, 2 Sep 2011 15:33:13 -0400 > > From: tzai...@alcor.concordia.ca > > To: r-help@r-project.org > > Subject: [R] Hessian Matrix Issue > > Message-ID: > > <e6dc43b4487eb4a4055e1ab485f015f0.squir...@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@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. ______________________________________________ R-help@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.