Bootstrap 用于置信区间

Question

我的问题如下：首先，我必须创建 1000 bootstrap 个大小为 100 的 "theta hat" 样本。我有一个随机变量 X，它遵循缩放 t_5-distribution。以下代码创建 1000 bootstrap 个 theta 帽子样本：

library("metRology", lib.loc="~/R/win-library/3.4")
 # Draw some data
data <- rt.scaled(100, df=5, mean=0, sd=2)

thetahatsq <- function(x){(3/500)*sum(x^2)}
sqrt(thetahatsq(data))

n <- 100

thetahat <- function(x){sqrt(thetahatsq(x))}
thetahat(data)

# Draw 1000 samples of size 100 from the fitted distribution, and compute the thetahat
tstar<-replicate(1000,thetahat(rt.scaled(n, df=5, mean=0, sd=thetahat(data)))) 
mean(tstar)

hist(tstar, breaks=20, col="lightgreen")

现在我想比较覆盖概率的准确性和使用百分位数方法构建的 95% bootstrap 置信区间的宽度。我想重复上面的代码1000次，每次都检查参数的真实值是否属于对应的bootstrap置信区间，并计算每个区间的长度。然后平均结果值。

Answer 1

也许 bootstrap 的最佳方法是使用基础包 boot。函数 boot 和 boot.ci 是您想要的，函数 boot.ci 为您提供有关要计算的置信区间类型的选项，包括 type = "perc".

看看以下是否回答了您的问题。

set.seed(402)    # make the results reproducible
data <- rt.scaled(100, df=5, mean=0, sd=2)

stat <- function(data, index) thetahat(data[index])

hans <- function(data, statistic, R){
    b <- boot::boot(data, statistic, R = R)
    ci <- boot::boot.ci(b, type = "perc")
    lower <- ci$percent[4]
    upper <- ci$percent[5]
    belongs <- lower <= true_val && true_val <= upper
    data.frame(lower, upper, belongs)
}

true_val <- sqrt(thetahatsq(data))

df <- do.call(rbind, lapply(seq_len(1000), function(i) hans(data, statistic = stat, R = n)))
head(df)
#     lower    upper belongs
#1 1.614047 2.257732    TRUE
#2 1.592893 2.144660    TRUE
#3 1.669754 2.187214    TRUE
#4 1.625061 2.210883    TRUE
#5 1.628343 2.220374    TRUE
#6 1.633949 2.341693    TRUE

colMeans(df)
#   lower    upper  belongs 
#1.615311 2.227224 1.000000

解释：

函数 stat 是您感兴趣的统计数据的包装器，供 boot 使用。
函数 hans 自动调用 boot::boot 和 boot::boot.ci。
对 hans 的调用是由 lapply 伪装的循环进行的。
返回的结果是data.frames的列表，所以我们需要调用do.call才能rbind将它们df.
其余为标准R代码。

Bootstrap 用于置信区间

Bootstrap for Confidence Intervals

statistics

replication

r

confidence-interval

statistics-bootstrap