Начальная загрузка MLE на распределении Пуассона

У меня есть следующий 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 ̂

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, и найти значение p, используя численную оптимизацию, которую я использовал выше. Я думал, что это будет что-то в этом роде

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

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


r machine-learning statistics poisson statistics-bootstrap
person Jamie Leigh    schedule 12.12.2016    source источник

Ответы (1)


arrow_upward
3
arrow_downward

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

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

(2) Согласно предложению Чжэюаня, вам нужно использовать optimize вместо optim для минимизации nll в одномерном пространстве параметров.

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 для вычисления MLE и одновременного построения доверительных интервалов.

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
person Sandipan Dey    schedule 12.12.2016