MLE Bootstrap auf Poisson-Verteilung

Question

Ich habe folgende Poisson Verteilung:

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 

habe ich den folgenden Code der MLE Θ & # x0302 zu finden;

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

Was ich tun möchte, ist ein Bootstrap-Konfidenzintervall schafft die Null-Hypothese zu testen, dass Θ = 1 auf Stufe 0,05 und finden Sie den p-Wert die numerische Optimierung mit, die ich oben verwendet. Ich dachte, es etwas entlang der Linien dieser

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)) 

sein würde, aber ich glaube nicht, dass dies die Optimierung nutzt, und ich bin nicht sicher, wie die p-Wert zu erhalten.

Answer 1

Einige Punkte:

(1) Für numerische Stabilität Sie können negative logarithmische Wahrscheinlichkeit anstelle von negativer Wahrscheinlichkeit zu berücksichtigen.

(2) Gemäß Zheyuan Vorschlag, müssen Sie optimize anstelle von optim für NLL Minimierung in dem 1-D Parameterraum verwenden.

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) Sie können mle2 von bbmle Paket verwenden, um die MLE und konstruieren, um die Vertrauensintervalle auf einmal zu berechnen.

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