如何在二项比例置信区间模拟中实现递增样本量循环

看起来你已经搭好了不错的模拟基础！要实现从n=10到100的样本量测试，我们可以通过循环遍历每个样本量，在每次迭代中复用你已经写好的模拟逻辑。下面一步步来优化和实现：

第一步：优化置信区间计算函数 #

先把上下限的计算合并成一个函数，让代码更简洁，也方便后续调用：

# 合并计算置信区间上下限的函数
ci_binom <- function(x, cl = 0.95) {
  n <- length(x)
  p_est <- mean(x)
  z <- abs(qnorm((1 - cl)/2))
  se <- sqrt(p_est * (1 - p_est)/n)
  return(c(lower = p_est - z*se, upper = p_est + z*se))
}

第二步：编写样本量循环逻辑 #

这里提供两种实现方式，你可以根据习惯选择：

方式1：使用for循环（更直观易读） #

# 设置模拟核心参数
true_p <- 0.4
reps <- 200  # 每个样本量的重复模拟次数
sample_sizes <- 10:100  # 要测试的样本量范围
coverage_results <- numeric(length(sample_sizes))  # 存储每个样本量的覆盖率

# 循环遍历每个样本量
for (i in seq_along(sample_sizes)) {
  n <- sample_sizes[i]
  # 生成模拟数据：reps个独立样本，每个样本包含n个0/1观测
  dat <- rbinom(reps * n, 1, true_p)
  # 把数据转成矩阵：每行是一个观测，每列是一个完整样本
  x_matrix <- matrix(dat, nrow = n, ncol = reps)
  
  # 对每个样本计算置信区间
  ci_results <- apply(x_matrix, 2, ci_binom)
  ll_res <- ci_results["lower", ]
  ul_res <- ci_results["upper", ]
  
  # 计算覆盖率：真实比例p落在置信区间内的样本占比
  hits <- ll_res <= true_p & true_p <= ul_res
  coverage_results[i] <- mean(hits)
}

# 把结果整理成数据框，方便查看和后续分析
results_df <- data.frame(
  sample_size = sample_sizes,
  coverage = coverage_results
)

# 查看前5个样本量的结果
head(results_df)

方式2：使用sapply（更紧凑简洁） #

如果你偏好更凝练的代码，可以用sapply替代for循环：

true_p <- 0.4
reps <- 200
sample_sizes <- 10:100

# 用sapply遍历每个样本量，直接返回覆盖率结果
coverage_results <- sapply(sample_sizes, function(n) {
  dat <- rbinom(reps * n, 1, true_p)
  x_matrix <- matrix(dat, nrow = n, ncol = reps)
  ci_results <- apply(x_matrix, 2, ci_binom)
  hits <- ci_results["lower", ] <= true_p & true_p <= ci_results["upper", ]
  mean(hits)
})

results_df <- data.frame(sample_size = sample_sizes, coverage = coverage_results)

第三步：可视化结果（可选但推荐） #

为了直观看到样本量对覆盖率的影响，可以画个折线图：

plot(results_df$sample_size, results_df$coverage, 
     type = "l", lwd = 2, col = "blue",
     xlab = "样本量n", ylab = "置信区间覆盖率",
     main = "二项比例置信区间覆盖率随样本量变化")
# 添加目标置信水平的参考线
abline(h = 0.95, col = "red", lty = 2)

小提示 #

• 注意矩阵维度：之前的代码里matrix(dat, ncol=rep)会让每行是一个样本，其实更合理的是让每列对应一个样本（每行是一个观测），这样apply(x_matrix, 2, ...)就是对每个完整样本计算置信区间。
• 如果模拟次数reps太小，结果会有随机波动，你可以适当增大（比如1000次）来得到更稳定的覆盖率估计。

内容的提问来源于stack exchange，提问作者Bea Hs