MLE самозагрузки по распределению Пуассона

Question

У меня есть следующий Poisson распределение:

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 

Я использовал следующий код, чтобы найти MLE Θ & # x0302;

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

То, что я хочу сделать, это создать загрузочный доверительный интервал, чтобы проверить нулевую гипотезу, что Θ = 1 на уровне 0,05 и найти р-значение с помощью численной оптимизации, которые я использовал выше. я думал, что это будет что-то вдоль линий этого

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

Но я не думаю, что это использует оптимизацию, и я не уверен, как получить значение р.

Answer 1

Несколько моментов:

(1) Для численной устойчивости вы можете, как рассмотреть отрицательный логарифм вероятности вместо отрицательного вероятности.

(2) По предложению Zheyuan, вы должны использовать Оптимизировать вместо из Optim для СДЛ минимизации в пространстве параметров 1-D.

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) Вы можете использовать mle2 из bbmle пакета для вычисления ОМП и построить доверительные интервалы сразу.

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