Wald区间覆盖率图表异常及代码计算失效问题问询

原代码的核心错误点 #

• Wald区间函数仅计算了置信上限，未生成完整的上下限区间
• 硬编码z值（2.576），未根据置信水平动态计算
• probs序列步长设置错误（0.8远大于区间范围0.2-0.4），导致仅生成一个数据点
• 覆盖率判断逻辑错误：始终判断是否覆盖第一个真实概率值，而非当前循环的真实概率；且仅用单个值（而非区间）判断覆盖情况
• 覆盖率结果存储错误：始终覆盖coverage[1]，未按真实概率对应存储

修正后的完整代码 #

# 修正后的Wald区间计算函数：生成完整的上下限区间，动态计算z值
wald_ci <- function(x, n, conf.level = 0.99) {
  p_hat <- x / n
  se <- sqrt(p_hat * (1 - p_hat) / n)
  z <- qnorm((1 + conf.level) / 2)  # 根据置信水平动态获取分位数
  lower <- p_hat - z * se
  upper <- p_hat + z * se
  return(c(lower, upper))
}

# 参数设置
numSamples <- 10000  # 模拟样本数
numTrials <- 10      # 每个样本的试验次数
probs <- seq(0.05, 0.95, 0.05)  # 真实概率范围，步长0.05保证足够数据点绘图
coverage <- numeric(length(probs))  # 预分配存储空间

# 循环计算每个真实概率对应的覆盖率
for (i in seq_along(probs)) {
  true_p <- probs[i]
  # 生成二项分布样本
  x <- rbinom(n = numSamples, size = numTrials, prob = true_p)
  isCovered <- logical(numSamples)  # 预分配逻辑向量
  
  for (j in 1:numSamples) {
    ci <- wald_ci(x = x[j], n = numTrials)
    # 判断真实概率是否落在区间内
    isCovered[j] <- (ci[1] < true_p) && (true_p < ci[2])
  }
  # 计算覆盖率并存储
  coverage[i] <- mean(isCovered) * 100
}

# 绘图：显示Wald区间的实际覆盖率
plot(probs, coverage, type = "l", ylim = c(70, 100), col = "blue", lwd = 2,
     frame.plot = FALSE, yaxt = 'n', main = "Coverage of Wald Confidence Interval",
     xlab = "True Population Proportion", ylab = "Coverage (%) for 99% CI")
# 添加目标覆盖率水平线
abline(h = 99, lty = 3, col = "maroon", lwd = 2)
# 设置y轴刻度
axis(side = 2, at = seq(70, 100, 5))
# 添加网格线辅助观察
grid(nx = NULL, ny = seq(70, 100, 5), lty = 2, col = "gray80")

关键修正说明 #

1. Wald函数修正：
   • 同时计算置信区间的上下限，返回长度为2的向量
   • 使用qnorm()动态计算对应置信水平的z分位数，避免硬编码
2. 真实概率序列修正：
   • 调整seq()的步长为0.05，覆盖0.05到0.95的范围，确保绘图时有足够多的点形成连续曲线
3. 覆盖率计算逻辑修正：
   • 循环中使用当前真实概率true_p判断覆盖情况，而非固定第一个值
   • 正确使用区间的上下限判断真实概率是否被包含
   • 预分配isCovered和coverage的存储空间，提升运行效率
4. 图表修正：
   • 统一标题和坐标轴标签（将错误的95% CI改为99%）
   • 调整y轴范围为70-100，更清晰展示Wald区间的覆盖率波动
   • 添加网格线，便于对比实际覆盖率与目标覆盖率（99%）

内容的提问来源于stack exchange，提问作者Mariana Adabuk