使用R语言Bootstrap方法计算单值置信区间返回0的问题

原代码的核心问题 #

1. 代码结构错误：my_function中return(median(data[index]))语句之后的for循环永远不会执行，且循环逻辑应放在函数外部而非内部。
2. Bootstrap对象错误：直接对单个值p_hat[i]做Bootstrap抽样，单个值的Bootstrap样本只能重复该值本身，无法产生变异，导致所有Bootstrap样本的中位数固定为该值，后续boot.ci的结果无法正确提取，最终intervals始终为初始的0矩阵。
3. 逻辑误解：Bootstrap需基于生成p_hat的原始数据集抽样，每次抽样后重新计算完整的p_hat向量，再针对每个位置的Bootstrap结果计算置信区间，而非对单个统计量值抽样。

修正后的代码（基于原始数据集的正确实现） #

假设original_data是生成p_hat的真实原始数据集（比如p_hat是该数据集ECDF的19个分位点取值），以下是符合统计逻辑的实现：

library(boot)

set.seed(123) # 设置随机种子保证结果可复现
# 替换为你的真实原始数据，这里用正态分布数据做示例
original_data <- rnorm(1000) 

# 定义计算p_hat的函数：输入样本，输出19个ECDF值
compute_p_hat <- function(data) {
  quantiles <- seq(0, 1, length.out = 19) # 19个等概率分位点
  ecdf_obj <- ecdf(data)
  return(ecdf_obj(quantiles) * 100) # 转换为百分比，匹配你提供的p_hat格式
}

# 初始的p_hat（基于原始数据的估计）
p_hat <- compute_p_hat(original_data)

# Bootstrap抽样函数：输入数据和索引，返回重新计算的p_hat向量
boot_function <- function(data, index) {
  boot_sample <- data[index]
  return(compute_p_hat(boot_sample))
}

# 执行Bootstrap抽样
B <- 1000
boot_result <- boot(data = original_data, statistic = boot_function, R = B)

# 计算每个p_hat元素的BCA置信区间
alpha <- 0.9
intervals <- matrix(NA, ncol = 2, nrow = length(p_hat))
for (i in 1:length(p_hat)) {
  ci <- boot.ci(boot_result, type = "bca", conf = alpha, index = i)
  intervals[i, ] <- ci$bca[4:5] # 提取BCA区间的上下限
}

# 整理结果
results <- data.frame(Estimate = p_hat, Lower_CI = intervals[,1], Upper_CI = intervals[,2])
print(results)

针对现有p_hat的演示性修正（无统计意义） #

如果没有原始数据，仅为让代码运行出结果，可构造伪样本（但此方法不符合Bootstrap统计逻辑，仅作演示）：

library(boot)

alpha <- 0.9
B <- 1000
p_hat <- c(0, 0, 6.70881, 14.16335, 26.08988, 41.33073, 57.23204, 
           74.02023, 88.54585, 95.48473, 98.97599, 99.90797, 
           99.94741, 100, 100, 100, 100, 100, 100)
intervals <- matrix(0, ncol = 2, nrow = length(p_hat))

# 定义Bootstrap函数
my_function <- function(data, index) {
  return(median(data[index]))
}

for (i in 1:length(p_hat)) {
  pseudo_sample <- rep(p_hat[i], 100) # 构造含重复值的伪样本
  boot_samples <- boot(pseudo_sample, my_function, R = B)
  ci <- boot.ci(boot_samples, type = "bca", conf = alpha)
  intervals[i, ] <- ci$bca[4:5]
}

results <- data.frame(p_hat, Lower_CI = intervals[,1], Upper_CI = intervals[,2])
print(results)

注意：第二种方法仅为代码运行演示，不具备统计学有效性，建议优先使用基于原始数据集的第一种方法。

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