【问题标题】:MLE bootstrap on poisson distribution泊松分布的 MLE 引导程序
【发布时间】:2017-04-26 21:27:33
【问题描述】:

我有以下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


    【解决方案1】:

    几点:

    (1) 为了数值稳定性,您可能希望考虑负对数似然而不是负似然。

    (2) 根据哲元的建议,一维参数空间的nll最小化需要使用optimize代替optim

    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 
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2018-02-10
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2015-06-13
      • 1970-01-01
      相关资源
      最近更新 更多