[R] Hessian Matrix Issue

[Uwe Ligges]

I have not really looked into the details of the lengthy and almost 
unreadable code below. In any case, there are good reasons why numerics 
software typically uses Fisher scoring / IWLS in order to fit GLMs.....

And if your matrix is that "singular", even the common numerical tricks 
may not save the day anymore. 7e-21 is very close to exact singularity!

Uwe Ligges






On 02.09.2011 21:33, tzaihra at alcor.concordia.ca wrote:
> 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.