R 中的 Bootstrap 重抽样
使用带放回的 sample() 实现 Bootstrap 置信区间
R 中的 Bootstrap 重抽样 是 CoddyKit 上的免费 R Academy 课时。 这是第 4 节课，共 4 节。 你可以在下方免费阅读本课时的完整内容 — 然后在浏览器中使用内置代码编辑器和全天候 AI 导师进行实践。 这是 R Academy 学习路径的一部分，你的进度在网页和 CoddyKit 应用中同步。 R Academy 课程共包含 4 节课。
什么是自助法？
自助法是一种重抽样方法，通过从观测数据中进行有放回抽样，估计统计量的抽样分布。它不需要分布假设——数据本身就是分布。
set.seed(42)
# Original sample of 20 observations
original_data <- c(14, 18, 11, 13, 6, 8, 2, 17,
20, 10, 15, 12, 9, 16, 7, 3,
19, 5, 4, 1)
n <- length(original_data)
# One bootstrap resample: sample WITH replacement
bootstrap_sample1 <- sample(original_data, n, replace = TRUE)
print(bootstrap_sample1)
# Contains repeats! Some original values are missing.
mean(original_data) # original mean
mean(bootstrap_sample1) # bootstrap mean (different)均值的自助法分布
生成大量自助法重抽样样本，分别计算每个样本的统计量，并收集结果。这个自助法分布可以近似统计量的抽样分布。
set.seed(42)
x <- c(14, 18, 11, 13, 6, 8, 2, 17, 20, 10,
15, 12, 9, 16, 7, 3, 19, 5, 4, 1)
n <- length(x)
B <- 5000 # number of bootstrap replicates
# Bootstrap distribution of the mean
boot_means <- replicate(B, {
x_star <- sample(x, n, replace = TRUE)
mean(x_star)
})
cat('Original mean:', mean(x), '\n')
cat('Bootstrap mean:', mean(boot_means), '\n')
cat('Bootstrap SE:', sd(boot_means), '\n')
cat('Theoretical SE:', sd(x)/sqrt(n))适用于任意统计量的自助法
自助法适用于任何统计量——中位数、相关性、回归系数、截尾均值。只需在重复循环中将 mean() 替换为您感兴趣的统计量即可。
set.seed(42)
x <- c(2, 4, 6, 8, 10, 12, 100, 1, 3, 5)
# Notice the outlier 100 - makes the mean unreliable
B <- 5000; n <- length(x)
# Bootstrap distribution of the median
boot_medians <- replicate(B, {
median(sample(x, n, replace = TRUE))
})
# Bootstrap distribution of the trimmed mean (10%)
boot_tmeans <- replicate(B, {
mean(sample(x, n, replace = TRUE), trim = 0.1)
})
cat('Sample median:', median(x), '\n')
cat('Bootstrap SE of median:', sd(boot_medians), '\n')
cat('Bootstrap SE of trimmed mean:', sd(boot_tmeans))自助法置信区间：百分位数法
最简单的自助法置信区间是：取自助法分布的分位数。95% 置信区间使用第 2.5 和第 97.5 百分位数。该方法快速直观，但对于偏斜统计量可能存在覆盖率问题。
set.seed(42)
x <- rnorm(30, mean = 5, sd = 2)
B <- 5000; n <- length(x)
# Bootstrap CI for the mean
boot_means <- replicate(B, mean(sample(x, n, replace = TRUE)))
# Percentile CI
ci_percentile <- quantile(boot_means, c(0.025, 0.975))
cat('Percentile CI:', round(ci_percentile, 3), '\n')
# For comparison: t-interval (parametric)
t_ci <- mean(x) + qt(c(0.025, 0.975), df = n-1) * sd(x)/sqrt(n)
cat('t-interval CI:', round(t_ci, 3), '\n')
cat('True mean: 5')自助法标准误
自助法标准误就是 sd(boot_statistics)。它估计从总体中反复抽样时统计量会有多大变化，而无需作出分布假设。
set.seed(42)
# Bootstrap SE for the ratio of two means
group_a <- rnorm(20, mean = 10, sd = 2)
group_b <- rnorm(20, mean = 8, sd = 2)
all_data <- data.frame(
value = c(group_a, group_b),
group = rep(c('A', 'B'), each = 20)
)
B <- 4000; n <- nrow(all_data)
boot_ratio <- replicate(B, {
idx <- sample(1:n, n, replace = TRUE)
d <- all_data[idx, ]
mean(d$value[d$group == 'A']) /
mean(d$value[d$group == 'B'])
})
cat('Observed ratio:', mean(group_a)/mean(group_b), '\n')
cat('Bootstrap SE:', round(sd(boot_ratio), 4), '\n')
cat('95% CI:', round(quantile(boot_ratio, c(.025,.975)), 3))BCa 置信区间
BCa（偏差校正与加速）置信区间比百分位数置信区间更准确，尤其适用于偏斜或有偏的统计量。它会校正自助法分布中的偏差和偏斜。
# BCa CI implementation
bca_ci <- function(data, stat_fn, B = 2000, alpha = 0.05) {
n <- length(data)
theta_hat <- stat_fn(data)
# Bootstrap replicates
boot_vals <- replicate(B, stat_fn(sample(data, n, replace = TRUE)))
# Bias correction z0
z0 <- qnorm(mean(boot_vals < theta_hat))
# Acceleration a (jackknife)
jk <- sapply(seq_len(n), function(i) stat_fn(data[-i]))
num <- sum((mean(jk) - jk)^3)
den <- 6 * (sum((mean(jk) - jk)^2))^(3/2)
a <- num / den
# Adjusted quantiles
z_alpha <- qnorm(c(alpha/2, 1 - alpha/2))
p_adj <- pnorm(z0 + (z0 + z_alpha) / (1 - a * (z0 + z_alpha)))
quantile(boot_vals, p_adj)
}
set.seed(42)
x <- rexp(25, rate = 0.5) # skewed distribution
bca_ci(x, median)用于回归的自助法
通过对数据框的行进行重抽样来构造回归置信区间。这种方法可以处理非正态误差，并且无需正态性假设即可进行有效推断。
set.seed(42)
n <- 40
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)
df <- data.frame(x = x, y = y)
B <- 3000
# Bootstrap the slope coefficient
boot_slopes <- replicate(B, {
idx <- sample(1:n, n, replace = TRUE)
d_boot <- df[idx, ]
coef(lm(y ~ x, data = d_boot))[2]
})
cat('OLS slope:', round(coef(lm(y~x, data=df))[2], 3), '\n')
cat('Bootstrap SE of slope:', round(sd(boot_slopes), 4), '\n')
cat('Bootstrap 95% CI:', round(quantile(boot_slopes, c(.025,.975)), 3))自助法重复次数
您需要多少次自助法重复？对于标准误：B=200-500。对于百分位数置信区间：B=1000-2000。对于 BCa 置信区间：B=2000-5000。增加次数通常更好，但超过 10,000 次后收益会逐渐降低。
set.seed(42)
x <- rnorm(50, mean = 3, sd = 1)
B_values <- c(100, 500, 1000, 2000, 5000, 10000)
# See how CI width varies with B
results <- sapply(B_values, function(B) {
boot_m <- replicate(B, mean(sample(x, length(x), replace = TRUE)))
q <- quantile(boot_m, c(0.025, 0.975))
q[2] - q[1] # CI width
})
result_df <- data.frame(B = B_values, ci_width = round(results, 4))
print(result_df)
# CI width stabilizes as B increases配对自助法
对于双样本问题，请将成对数据 (x, y) 一起重抽样，而不是分别对 x 和 y 抽样。这样可以保留配对内的依赖结构，为相关性和配对差异提供有效的置信区间。
set.seed(42)
# Paired data: before/after treatment
before <- c(10, 12, 8, 15, 11, 9, 14, 13)
after <- c(12, 14, 9, 16, 11, 10, 15, 12)
n <- length(before)
B <- 4000
# Bootstrap paired mean difference
boot_diff <- replicate(B, {
idx <- sample(1:n, n, replace = TRUE)
mean(after[idx] - before[idx])
})
observed_diff <- mean(after - before)
cat('Observed mean diff:', observed_diff, '\n')
cat('Bootstrap SE:', round(sd(boot_diff), 4), '\n')
cat('95% CI:', round(quantile(boot_diff, c(0.025, 0.975)), 3))自助法假设检验
自助法也可以用于检验假设。要检验 H₀: θ=θ₀，请通过平移或重新居中自助法样本来生成原假设分布，然后将超过观测统计量的比例作为 p 值。
set.seed(42)
# Test H0: mean = 0 vs H1: mean != 0
x <- rnorm(30, mean = 0.5, sd = 2)
theta0 <- 0 # null hypothesis
observed_t <- (mean(x) - theta0) / (sd(x) / sqrt(length(x)))
B <- 5000; n <- length(x)
# Shift x to have mean = theta0 under H0
x_centered <- x - mean(x) + theta0
boot_t <- replicate(B, {
x_star <- sample(x_centered, n, replace = TRUE)
(mean(x_star) - theta0) / (sd(x_star) / sqrt(n))
})
p_value <- mean(abs(boot_t) >= abs(observed_t))
cat('Observed t:', round(observed_t, 3), '\n')
cat('Bootstrap p-value:', round(p_value, 4))Boot 程序包概览
boot 程序包提供 boot(data, statistic, R) 及多种高级功能：通过 boot.ci() 计算 BCa 置信区间、并行计算和分层抽样。statistic 函数接收数据和索引向量。
# library(boot) # uncomment to use
# The boot() statistic function signature:
# statistic(data, indices) -> scalar or vector
# Example (conceptual - requires boot package):
# mean_stat <- function(data, indices) {
# mean(data[indices])
# }
# results <- boot(data = x, statistic = mean_stat, R = 5000)
# boot.ci(results, type = c('perc', 'bca')) # CIs
# plot(results) # histogram of bootstrap distribution
# Manual equivalent:
set.seed(42)
x <- rnorm(30, mean = 5)
B <- 2000
boot_vals <- replicate(B, mean(sample(x, length(x), replace = TRUE)))
quantile(boot_vals, c(0.025, 0.975))快速检查
测试您对 R 中自助法重抽样的理解。
回顾：自助法重抽样
要点：自助法使用 sample(x, n, replace=TRUE) 对数据进行有放回重抽样。使用 replicate(B, ...) 生成 B 个自助法统计量。sd(boot_stats) 给出自助法标准误。百分位数置信区间：quantile(boot_stats, c(0.025, 0.975))。计算置信区间时请使用 B≥2000。boot 程序包增加了 BCa 置信区间和并行计算功能。自助法无需分布假设即可适用于任何统计量。
set.seed(42)
x <- c(3, 5, 7, 2, 9, 4, 6, 8, 1, 10)
n <- length(x); B <- 3000
# Generic bootstrap template:
boot_stat <- replicate(B, {
x_star <- sample(x, n, replace = TRUE)
median(x_star) # replace with any statistic
})
# Summary
cat('Original median:', median(x), '\n')
cat('Bootstrap SE:', round(sd(boot_stat), 3), '\n')
ci <- quantile(boot_stat, c(0.025, 0.975))
cat('95% Percentile CI: [', ci[1], ',', ci[2], ']')常见问题解答
「R 中的 Bootstrap 重抽样」课时是免费的吗？
是的 — 「R 中的 Bootstrap 重抽样」的完整文本可在网页上免费阅读。要进行交互式练习（内置代码编辑器和全天候 AI 导师）并解锁 R Academy 课程的其余内容，请升级到 CoddyKit PRO。 R Academy 课程共包含 4 节课。
「R 中的 Bootstrap 重抽样」这节课中我会学到什么？
使用带放回的 sample() 实现 Bootstrap 置信区间 你通过在浏览器中直接运行的动手代码来练习 R Academy，全天候 AI 导师会在你学习这节课的过程中回答你的问题。
学习 R Academy 需要有经验吗？
无需任何先前经验。CoddyKit 上的 R Academy 课程适合初学者到高级学习者，你可以从这里开始或从头开始，按照自己的节奏学习。 这是第 4 节课，共 4 节。
「R 中的 Bootstrap 重抽样」课时需要多长时间？
大多数 CoddyKit 课程大约需要 5–10 分钟。每节课都很精短且互动，所以你能稳步进步，并在网页和应用中从离开的地方继续。
我能在这节 R Academy 课中编写并运行代码吗？
能。每节 R Academy 课都包含内置代码编辑器，你可以在浏览器中直接编写并运行真实代码，并获得即时 AI 反馈 — 无需本地设置。
此课程中的所有课时
- set.seed() 与可复现性
- 生成随机分布
- 蒙特卡洛模拟基础
- R 中的 Bootstrap 重抽样