2016-12-12 181 views
2

我有以下Poisson分佈:泊松分佈的MLE引導

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值。

回答

2

的幾點:

(1)對於數值的穩定性,你可能願意考慮負對數似然,而不是消極的可能性。

(2)根據宋哲元的建議下,則需要在1-d參數空間使用優化代替的的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)您可以使用bbmle包中的mle2來計算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