MLE bootstrap on poisson distribution

Question

I have the following Poisson distribution:

Data
3 5 3 1 2 1 2 1 0 2 4 3 1 4 1 2 2 0 4 2 2 4 0 2 1 0 5 2 0 1 
2 1 3 0 2 1 1 2 2 0 3 2 1 1 2 2 5 0 4 3 1 2 3 0 0 0 2 1 2 2 
3 2 4 4 2 1 4 3 2 0 3 1 2 1 3 2 6 0 3 5 1 3 0 1 2 0 1 0 0 1 
1 0 3 1 2 3 3 3 2 1 1 2 3 0 0 1 5 1 1 3 1 2 2 1 0 3 1 0 1 1

I used the following code to find the MLE Θ̂

lik<-function(lam) prod(dpois(data,lambda=lam)) #likelihood function
nlik<- function(lam) -lik(lam) #negative-likelihood function
optim(par=1, nlik) 

What I want to do is create a bootstrap confidence interval to test the null hypothesis that Θ = 1 at level 0.05 and find the p-value using the numerical optimization that I used above. I thought it would be something along the lines of this

n<-length(data)
nboot<-1000
boot.xbar <- rep(NA, nboot)
for (i in 1:nboot) {
data.star <- data[sample(1:n,replace=TRUE)]
boot.xbar[i]<-mean(data.star)
}
quantile(boot.xbar,c(0.025,0.975))

But I don't think this utilizes the optimization, and I am not sure how to get the p-value.

Answer 1

A few points:

(1) For numerical stability you may like to consider negative log likelihood instead of negative likelihood.

(2) As per Zheyuan's suggestion, you need to use optimize in stead of optim for nll minimization in the 1-D parameter space.

lik<-function(lam) sum(log(dpois(data,lambda=lam))) #log likelihood function
nlik<- function(lam) -lik(lam) #negative-log-likelihood function
optimize(nlik, c(0.1, 2), tol = 0.0001)

# $minimum
# [1] 1.816661    
# $objective
# [1] 201.1172

n<-length(data)
nboot<-1000
boot.xbar <- rep(NA, nboot)
for (i in 1:nboot) {
  data.star <- data[sample(1:n,replace=TRUE)]
  boot.xbar[i]<-mean(data.star)
}
quantile(boot.xbar,c(0.025,0.975))
#  2.5%    97.5% 
# 1.575000 2.066667 

(3) You may use mle2 from bbmle package to compute the MLE and construct the confidence intervals at once.

library(bbmle)
res <- mle2(minuslogl = nlik, start = list(lam = 0.1))
res   
# Coefficients:
#     lam 
# 1.816708     
# Log-likelihood: -201.12 

confint(profile(res)) # confint w.r.t. the likelihood profile
# 2.5 %   97.5 % 
# 1.586083 2.068626 

confint(res, method="uniroot") # based on root-finding to find the exact point where the profile crosses the critical level     
#    2.5 %   97.5 % 
# 1.586062 2.068606