机器算法验证 - 基本自举置信区间的覆盖概率 - 吾爱随笔录

基本自举置信区间的覆盖概率

机器算法验证 r 置信区间自习引导程序蒙特卡洛

2022-03-28 02:40:27

对于我正在学习的课程，我有以下问题：

进行蒙特卡罗研究以估计标准正态自举置信区间和基本自举置信区间的覆盖概率。从正常人群中抽样并检查样本均值的经验覆盖率。

标准正态引导 CI 的覆盖概率很简单：

n = 1000;
alpha = c(0.025, 0.975);
x = rnorm(n, 0, 1);
mu = mean(x);
sqrt.n = sqrt(n);

LNorm = numeric(B);
UNorm = numeric(B);

for(j in 1:B)
{
    smpl = x[sample(1:n, size = n, replace = TRUE)];
    xbar = mean(smpl);
    s = sd(smpl);

    LNorm[j] = xbar + qnorm(alpha[1]) * (s / sqrt.n);
    UNorm[j] = xbar + qnorm(alpha[2]) * (s / sqrt.n);
}

mean(LNorm < 0 & UNorm > 0); # Approximates to 0.95
# NOTE: it is not good enough to look at overall coverage
# Must compute separately for each tail

根据我在本课程中学到的知识，基本的引导置信区间可以这样计算：

# Using x from previous...
R = boot(data = x, R=1000, statistic = function(x, i){ mean(x[i]); });
result = 2 * mu - quantile(R$t, alpha, type=1);

这就说得通了。我不明白的是如何计算基本引导 CI 的覆盖概率。我知道覆盖概率将代表 CI 包含真实值的次数（在这种情况下mu）。我只是boot多次运行该功能吗？

我怎样才能以不同的方式处理这个问题？

1个回答

该术语的使用可能不一致，因此以下仅是我对原始问题的理解。据我了解，您计算的正常 CI 不是所要求的。每组 bootstrap 复制都给你一个置信区间，而不是很多。从一组引导复制的结果中计算不同 CI 类型的方法如下：

B    <- 999                  # number of replicates
muH0 <- 100                  # for generating data: true mean
sdH0 <- 40                   # for generating data: true sd
N    <- 200                  # sample size
DV   <- rnorm(N, muH0, sdH0) # simulated data: original sample

由于我想将计算结果与 package 的结果进行比较boot，我首先定义了一个函数，它将为每个复制调用。它的参数是原始样本和指定单个复制案例的索引向量。它返回的插件估计，以及，平均值方差的插件估计。 -CI需要后者。 $M^{\star}$ $\mu$ $S_{M}^{2\star}$ $\sigma_{M}^{2}$ $t$

> getM <- function(orgDV, idx) {
+     bsM   <- mean(orgDV[idx])                       # M*
+     bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N    # S^2*(M)
+     c(bsM, bsS2M)
+ }

> library(boot)                                       # for boot(), boot.ci()
> bOut <- boot(DV, statistic=getM, R=B)
> boot.ci(bOut, conf=0.95, type=c("basic", "perc", "norm", "stud"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL : 
boot.ci(boot.out = bOut, conf = 0.95, type = c("basic", "perc", "norm", "stud"))

Intervals : 
Level      Normal            Basic         Studentized        Percentile    
95%   ( 95.6, 106.0 )   ( 95.7, 106.2 )  ( 95.4, 106.2 )   ( 95.4, 106.0 )  
Calculations and Intervals on Original Scale

在不使用包的情况boot下，您可以简单地使用replicate()来获取一组引导程序副本。

boots <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))

但是，让我们坚持使用结果boot.ci()以供参考。

boots   <- bOut$t                     # estimates from all replicates
M       <- mean(DV)                   # M from original sample
S2M     <- (((N-1)/N) * var(DV)) / N  # S^2(M) from original sample
Mstar   <- boots[ , 1]                # M* for each replicate
S2Mstar <- boots[ , 2]                # S^2*(M) for each replicate
biasM   <- mean(Mstar) - M            # bias of estimator M

basic、percentile 和 -CI 依赖于 bootstrap 估计的经验分布。为了获得和分位数，我们找到引导估计的排序向量的相应索引（请注意，当索引不是自然数时，将进行更复杂的插值以找到经验分位数） . $t$ $\alpha/2$ $1 - \alpha/2$ boot.ci()

(idx <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2)) # indices for sorted vector of estimates
[1] 25 975

> (ciBasic <- 2*M - sort(Mstar)[idx])          # basic CI
[1] 106.21826  95.65911

> (ciPerc <- sort(Mstar)[idx])                 # percentile CI
[1] 95.42188 105.98103

对于 -CI，我们需要引导估计来计算关键的。对于标准正态 CI，临界值就是标准正态分布的值。 $t$ $t^{\star}$ $t$ $z$

# standard normal CI with bias correction
> zCrit   <- qnorm(c(0.025, 0.975))   # z-quantiles from std-normal distribution
> (ciNorm <- M - biasM + zCrit * sqrt(var(Mstar)))
[1] 95.5566 106.0043

> tStar <- (Mstar-M) / sqrt(S2Mstar)  # t*
> tCrit <- sort(tStar)[idx]           # t-quantiles from empirical t* distribution
> (ciT  <- M - tCrit * sqrt(S2M))     # studentized t-CI
[1] 106.20690  95.44878

为了估计这些 CI 类型的覆盖概率，您必须多次运行此模拟。只需将代码包装到一个函数中，返回一个包含 CI 结果的列表，然后按照这个 gistreplicate()中演示的方式运行它。

其它你可能感兴趣的问题

上一篇如何为留一法交叉验证生成 ROC 曲线？下一篇如何制作一个好的颜色强度刻度？