[R] Hessian Matrix Issue

[tzaihra at alcor.concordia.ca]

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