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.