为什么我的引导间隔覆盖率很差?

机器算法验证 引导程序 诊断
2022-01-21 03:48:20

我想做一个课堂演示,我将 t 间隔与引导间隔进行比较并计算两者的覆盖概率。我希望数据来自偏态分布,因此我选择将数据生成为exp(rnorm(10, 0, 2)) + 1,一个大小为 10 的样本来自移位的对数正态分布。我编写了一个脚本来抽取 1000 个样本,并且对于每个样本,基于 1000 个重复计算 95% 的 t 间隔和 95% 的引导百分位数间隔。

当我运行脚本时,两种方法都给出了非常相似的间隔,并且都具有 50-60% 的覆盖概率。我很惊讶,因为我认为引导间隔会更好。

我的问题是,我有没有

  • 在代码中犯了错误?
  • 在计算间隔时犯了错误?
  • 错误地期望引导间隔具有更好的覆盖属性?

另外,有没有办法在这种情况下构建更可靠的 CI?

 tCI.total <- 0
 bootCI.total <- 0
 m <- 10 # sample size
 true.mean <- exp(2) + 1

for (i in 1:1000){
 samp <- exp(rnorm(m,0,2)) + 1
 tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)

 boot.means <- rep(0,1000)
 for (j in 1:1000) boot.means[j] <- mean(sample(samp,m,replace=T))
 bootCI <- sort(boot.means)[c(0.025*length(boot.means), 0.975*length(boot.means))]

 if (true.mean > min(tCI) & true.mean < max(tCI)) tCI.total <- tCI.total + 1
 if (true.mean > min(bootCI) & true.mean < max(bootCI)) bootCI.total <- bootCI.total + 1 
}
tCI.total/1000     # estimate of t interval coverage probability
bootCI.total/1000  # estimate of bootstrap interval coverage probability
4个回答

Canto, Davison, Hinkley & Ventura (2006) 的Bootstrap 诊断和补救措施似乎是一个合乎逻辑的出发点。他们讨论了引导程序崩溃的多种方式,并且在这里更重要的是提供了诊断和可能的补救措施:

  1. 异常值
  2. 重采样模型不正确
  3. 非关键性
  4. 引导方法不一致

在这种情况下,我认为 1、2 和 4 没有问题。让我们看一下 3。正如@Ben Ogorek 指出的那样(尽管我同意@Glen_b 关于正态性讨论可能是一个红鲱鱼),引导程序的有效性取决于我们感兴趣的统计数据的关键性。

Canty 等人的第 4 节。建议重新采样内重新采样以获得每个引导重采样内参数估计的偏差和方差的度量。这是从 p 复制公式的代码。第15条:

library(boot)
m <- 10 # sample size
n.boot <- 1000
inner.boot <- 1000

set.seed(1)
samp.mean <- bias <- vars <- rep(NA,n.boot)
for ( ii in 1:n.boot ) {
    samp <- exp(rnorm(m,0,2)) + 1
    samp.mean[ii] <- mean(samp)
    foo <- boot(samp,statistic=function(xx,index)mean(xx[index]),R=inner.boot)
    bias[ii] <- mean(foo$t[,1])-foo$t0
    vars[ii] <- var(foo$t[,1])
}

opar <- par(mfrow=c(1,2))
    plot(samp.mean,bias,xlab="Sample means",ylab="Bias",
        main="Bias against sample means",pch=19,log="x")
    abline(h=0)
    plot(samp.mean,vars,xlab="Sample means",ylab="Variance",
        main="Variance against sample means",pch=19,log="xy")
par(opar)

引导诊断

注意对数刻度 - 没有对数,这更加明显。我们很好地看到了 bootstrap 均值估计的方差如何随 bootstrap 样本的均值上升。在我看来,这足以证明非关键性是低置信区间覆盖率的罪魁祸首。

但是,我很高兴地承认,可以通过多种方式跟进。例如,我们可以查看特定引导复制的置信区间是否包括真实均值取决于特定复制的均值。

至于补救措施,Canty 等人。讨论转换,这里会想到对数(例如,引导和建立置信区间不是为了平均值,而是为了记录数据的平均值),但我无法真正让它发挥作用。

坎蒂等人。继续讨论如何通过重要性采样和平滑以及向枢轴图添加置信带来减少内部引导程序的数量和剩余噪声。

对于聪明的学生来说,这可能是一个有趣的论文项目。我会很感激任何指出我出错的地方以及任何其他文献的指针。我会冒昧地将diagnostic标签添加到这个问题中。

虽然我同意 Stephan Kolassa 的分析和结论,

μ^μ
μ^样本均值绝对不是一个近似的支点,让我再补充一句。我调查了使用t-统计
mμ^μσ^
连同自举。结果是覆盖率约为 0.8。不是一个完整的解决方案,而是一个改进。

然后我对整个设置进行了更多思考。只有 10 个观测值和极其偏斜的分布,那么非参数估计均值基本上是不可能的,更不用说构建具有正确覆盖范围的置信区间了?

考虑的对数正态分布具有均值e2+1=8.39. 自从P(X2)=0.84什么时候XN(0,4)平均值是0.84-分布的分位数!这意味着所有 10 个观测值都小于平均值的概率是0.8410=0.178. 因此,在略低于 18% 的情况下,最大观察值小于平均值。为了获得大于 0.82 的覆盖率,我们需要构建超出最大观测值的平均值的置信区间。我很难想象如何在没有事先假设分布极度偏斜的情况下做出这样的构造(和证明是合理的)。但我欢迎任何建议。

计算是正确的,我与著名的包boot进行了交叉检查。此外,我添加了 BCa 间隔(由 Efron 提供),这是百分位引导间隔的偏差校正版本:

for (i in 1:1000) {
  samp <- exp(rnorm(m, 0, 2)) + 1

  boot.out <- boot(samp, function(d, i) sum(d[i]) / m, R=999)
  ci <- boot.ci(boot.out, 0.95, type="all")

  ##tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)
  tCI <- ci$normal[2:3]
      percCI <- ci$perc[4:5]
  bcaCI <- ci$bca[4:5]
      boottCI <- ci$student[4:5]

  if (true.mean > min(tCI) && true.mean < max(tCI)) tCI.total <- tCI.total + 1
  if (true.mean > min(percCI) && true.mean < max(percCI)) percCI.total <- percCI.total + 1 
  if (true.mean > min(bcaCI) && true.mean < max(bcaCI)) bcaCI.total <- bcaCI.total + 1
}

tCI.total/1000     # estimate of t interval coverage probability
0.53
percCI.total/1000  # estimate of percentile interval coverage probability
0.55
bcaCI.total/1000  # estimate of BCa interval coverage probability
0.61

我假设如果原始样本量大于 10,例如 20 或 50,则间隔会好得多。

此外,bootstrap-t方法通常会为偏态统计带来更好的结果。然而,它需要一个嵌套循环,因此需要多 20 倍的计算时间。

对于假设检验,良好的单面覆盖也非常重要。因此,仅查看 2 面覆盖范围通常会产生误导。

我也对此感到困惑,我在 1996 年 DiCiccio 和 Efron 的论文Bootstrap Confidence Intervals 上花了很多时间,但没有多少可以展示。

它实际上使我不再将引导程序视为一种通用方法。我曾经认为它可以在你真正陷入困境时将你从困境中拉出来。但我知道了它肮脏的小秘密:自举置信区间都以某种方式基于正态性。请允许我解释一下。

bootstrap 为您提供了估计器的抽样分布的估计值,这是您所希望的,对吧?但回想一下,抽样分布和置信区间之间的经典联系是基于找到一个关键量对于生锈的人,请考虑以下情况

xN(μ,σ2)
σ是已知的。然后数量
z=xμσN(0,1)
是关键的,即它的分布不依赖于μ. 所以,Pr(1.96xμσ1.96)=0.95剩下的就是历史。

当您考虑什么证明正态分布的百分位数与置信区间相关时,它完全基于这个方便的关键数量。对于任意分布,抽样分布的百分位数和置信区间之间没有理论上的联系,并且采用引导抽样分布的原始比例并不能削减它。

因此,Efron 的 BCa(偏差校正)区间使用转换来接近正态性,而 bootstrap-t 方法依赖于产生的 t 统计量是近似关键的。现在 bootstrap 可以估计出地狱的时刻,并且您始终可以假设正态性并使用标准 +/-2*SE。但是考虑到使用引导程序进入非参数化的所有工作,这似乎不太公平,不是吗?