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.
