Discussion: View Thread

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 my function 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, I believe 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 as to how to avoid such error msgs.


#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
 

Thank you and enjoy the long weekend.
Best Regards,
-------------------------------------------
[Tasneem] [Zaihra]
[Assistant Professor]
[UNB-SJ]
-------------------------------------------

optim() (in R) can get stuck.  For a concrete example of this (finding an M-estimator) and a detailed analysis of the problem, please read the thread at http://stats.stackexchange.com/questions/7308/can-the-empirical-hessian-of-an-m-estimator-be-indefinite/7629#7629 .

Best,
Bill

-------------------------------------------
William Huber
Quantitative Decisions
Rosemont, PA
-------------------------------------------