Qfyilyi

I am trying to compute observed 95% Wald Confidence intervals centred at the maximum likelihood estimate for realisations of i.i.d Zero Inflated Poisson(p, lambda) random variables (for both p and lambda). I also am trying to use simulations to compute the coverage of these intervals. At the moment I am generating my simulations from (alpha, n, p, lambda)=(0.05, 100, 0.2, 1.5). I am creating multiple functions in r but am wondering if there is an easier way?
This is the code I have so far for the Wald Confidence Intervals

compute.ml.estimates <- function(xs) {





 z <- ifelse(xs>0, 0, 1)

  n0 <- sum(z)

  n <- length(xs)

  g <- function(ab) {

    p <- ab[1]

    lambda <- ab[2]

    n0*log(p + (1-p)*exp(-lambda)) + (n-n0)*(log(1-p)-lambda) + log(lambda)*sum(xs)

}

  optim.out <- optim(c(0.1, 0.5), fn=g, method="L-BFGS-B", lower=c(0,0), upper=c(Inf, Inf), control=list(fnscale=-1))

  p.hat <- optim.out$par[1]

  lambda.hat <- optim.out$par[2]

  return(c(p.hat, lambda.hat))

}

compute.ml.estimates(data)

theta <- compute.ml.estimates(data)



Fisher.info <- function(theta) {



  p <- theta[1]

  lambda <- theta[2]

  phi <- p + (1-p)*exp(-lambda)



Fisher.info.matrix <- matrix(c(((1-exp(-lambda))/((1-p)*phi)), 

                              (-(exp(-lambda)/phi)), 

                              (-(exp(-lambda)/phi)), 

                              ((1-p)*(phi-p*lambda*exp(-lambda))/lambda*phi)),

                              nrow=2, byrow=TRUE)



return(Fisher.info.matrix)

}



phi



Fisher.info(data)



solve(Fisher.info(data))

but these are the answers

[1] 0.04480545 2.42678858

[1] 0.1291682

     [,1] [,2]

[1,]    0   -1

[2,]   -1  Inf

     [,1] [,2]

[1,]  NaN   -1

[2,]  NaN    0

which is obviously not ideal. I can't tell why it gives NaN