如何在R中向量化这些Bootstrap循环？

我刚从VB转R，习惯用循环写代码，但知道R里向量化处理数据效率更高。现在有一段Bootstrap的三层循环代码，想做向量化优化，先说明我的需求流程：

核心流程 #

针对样本量n=3:N，执行以下步骤：

• 从规模为N的原始样本中无放回抽取规模为n的随机子样本
• 对该子样本做B次重采样，完成Bootstrap参数估计（比如均值、标准差）
• 重复上述两步X次
• 对X次的参数估计结果取平均值，后续还要检查收敛性（比如看估计值的标准差）

目前代码只完成了前3步，第4步打算用rowMeans()在循环外实现。测试时B和X设为100，最终会用到10000甚至更大的值。

原始循环代码 #

# 模拟N=30的观测数据
bdf <- data.frame(sample(8:13, 30, rep = TRUE))

# 获取观测数
N <- length(bdf)

# 设置Bootstrap重复次数
B <- 100

# 设置重复估计的次数
X <- 100

# 创建空结果容器
result_vec <- vector(length=B)

# 外层循环：重复X次估计
for (j in 1:X) {
  
  # 中层循环：遍历样本量n=3到N
  for (i in 3:N) {
    
    # 无放回抽取n规模子样本
    boot_samp <- bdf[sample(N, size=i, replace=FALSE)]
    
    # 内层循环：B次Bootstrap重采样
    for(b in 1:B) {
      # 有放回抽取Bootstrap样本
      bsamp <- sample(boot_samp, size=i, replace=TRUE)
      
      # 计算参数（这里用均值，也可以换成标准差）
      p <- mean(bsamp)
      #p <- sd(bsamp)
      
      # 保存结果
      result_vec[b] <- p
      
    }
    # 绑定结果到数据框
    if (i==3) {
      df_res <- data.frame(result_vec)
    }
    else {
      df_temp <- data.frame(result_vec)
      df_res <- cbind(df_res, df_temp)
    }
    
    # 重命名列
    names(df_res)[ncol(df_res)] <- paste("n = ",i)
  }
  
  # 计算每列的均值
  allmeans <- colMeans(df_res)
  
  # 保存X次的均值结果
  if (j==1) {
    df_means <- data.frame(allmeans)
  }
  else {
    df_temp <- data.frame(allmeans)
    df_means <- cbind(df_means, df_temp)
  }
  
  # 重命名列
  names(df_means)[ncol(df_means)] <- j
}

向量化优化方案 #

原始三层循环在B和X很大时会非常慢，主要问题是循环内反复做数据框绑定、逐次赋值。下面用向量化采样+批量计算的方式优化，利用R的矩阵运算和apply系列函数减少循环层级：

优化后代码 #

# 模拟数据（设置随机种子保证结果可复现）
set.seed(123)
bdf <- sample(8:13, 30, replace = TRUE) # 直接用向量，比数据框计算效率更高
N <- length(bdf)
B <- 100
X <- 100

# 定义Bootstrap参数计算函数：输入子样本，返回B次重采样的参数（均值）
boot_param <- function(sub_samp) {
  n <- length(sub_samp)
  # 一次性生成B次重采样的索引矩阵，每行对应一次重采样的位置
  boot_indices <- matrix(sample(seq_len(n), size = n*B, replace = TRUE), nrow = B)
  # 批量计算每行的均值
  apply(boot_indices, 1, function(idx) mean(sub_samp[idx]))
}

# 外层：处理X次重复，每次生成所有n=3:N的子样本的Bootstrap均值
df_means <- do.call(cbind, lapply(1:X, function(j) {
  # 中层：遍历每个n，计算B次Bootstrap结果的均值
  sapply(3:N, function(n) {
    # 无放回抽取n规模子样本
    sub_samp <- sample(bdf, size = n, replace = FALSE)
    # 计算B次Bootstrap结果的均值
    mean(boot_param(sub_samp))
  })
}))

# 重命名行和列，提升可读性
rownames(df_means) <- paste("n =", 3:N)
colnames(df_means) <- paste("X =", 1:X)

优化点说明 #

1. 用向量替代数据框：原始代码用data.frame存储样本，向量的索引和计算效率更高，减少不必要的结构开销。
2. 批量生成重采样索引：在boot_param函数中一次性生成B×n的索引矩阵，避免逐次循环采样，直接通过矩阵索引批量获取样本。
3. 减少循环层级：用lapply和sapply替代嵌套循环，R内置的apply系列函数底层由C实现，比纯R循环速度快很多。
4. 避免循环内数据框绑定：直接用do.call(cbind, ...)一次性合并所有X次的结果，避免循环里反复cbind的性能损耗。
5. 函数化封装：把Bootstrap核心计算逻辑封装成函数，代码更清晰，修改参数（比如换成标准差）只需调整函数内的计算逻辑。

超大规模场景优化（B/X=10000+） #

如果B和X达到10000级，可以用parallel包做并行计算，利用多核CPU加速：

library(parallel)
# 获取可用核心数（留1个给系统）
num_cores <- detectCores() - 1

# 并行处理X次重复任务
df_means_par <- do.call(cbind, mclapply(1:X, function(j) {
  sapply(3:N, function(n) {
    sub_samp <- sample(bdf, size = n, replace = FALSE)
    mean(boot_param(sub_samp))
  })
}, mc.cores = num_cores))

rownames(df_means_par) <- paste("n =", 3:N)
colnames(df_means_par) <- paste("X =", 1:X)

收敛性检查（第4步实现） #

现在可以直接对df_means的行计算统计量，查看不同n下估计值的收敛情况：

# 计算每个n对应的X次估计值的均值和标准差
convergence_check <- data.frame(
  n = 3:N,
  mean_estimate = rowMeans(df_means),
  sd_estimate = apply(df_means, 1, sd)
)

print(convergence_check)

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