R语言中replicate()的高效替代方案：蒙特卡洛模拟加速

问题背景 #

现有一段蒙特卡洛模拟代码，针对含11628行的final.distrib数据框，逐行计算损失分布的95分位数。测试显示5行耗时50秒，预估总耗时超32小时，急需高效替代方案。

原代码核心问题 #

1. 函数冗余：6个建筑类型的模拟函数逻辑完全一致，仅参数不同，重复代码导致维护和执行效率低下
2. 重复随机数生成：每个函数每次调用都生成10万条对数正态分布数据再筛选，重复计算浪费资源
3. 低效循环：使用replicate逐次生成单建筑模拟结果，再用rowSums累加，循环开销大；apply逐行处理为单线程，未利用多核CPU
4. 种子设置不当：set.seed(1)放在每行计算的函数内，重复初始化种子，既影响效率也没必要

优化方案 #

1. 参数化统一函数 #

将6个重复函数合并为一个带参数的函数，传入各建筑类型的成本、损失分布参数，避免代码冗余。

2. 预生成复用随机数 #

一次性生成足够多的符合要求的面积数据，后续所有模拟直接从中采样，避免重复生成和筛选。

3. 向量化替代循环 #

直接生成批量模拟结果矩阵，用矩阵乘法替代replicate+rowSums，大幅减少循环开销。

4. 并行化处理多行数据 #

利用foreach和doParallel包实现多线程并行，同时处理多行数据，充分利用多核CPU资源。

优化后完整代码 #

第一步：预生成复用数据 #

# 加载所需包
library(partitions)
library(triangle)
library(foreach)
library(doParallel)

# 生成final.distrib数据框（原代码保留，此处省略重复部分）
# ...（原final.distrib生成代码不变）

# 预生成符合条件的面积数据，一次性生成足够多的样本
set.seed(1)
total_area_samples <- 1e6  # 生成100万条，足够所有模拟使用
random_area_all <- rlnorm(total_area_samples, meanlog = 3.61, sdlog = 0.82)
valid_area <- random_area_all[random_area_all >=5 & random_area_all <=200]

# 定义各建筑类型的参数表
building_params <- data.frame(
  cost_center = c(8407.29, 10328.43, 12055.57, 17935.86, 19363.57, 25182.71),
  loss_a = c(0.0890, 0.0423, 0.0453, 0.0645, 0.0263, 0.0583),
  loss_b = c(0.6655, 0.5572, 0.5665, 0.6361, 0.5009, 0.6078),
  loss_c = c(0.3212, 0.2038, 0.2136, 0.2791, 0.1537, 0.2519),
  row.names = c("LWA", "LWB", "SCA", "SCB", "RCA", "RCB")
)

第二步：参数化模拟函数 #

# 单建筑类型的批量模拟函数
simulate_building <- function(n_sim, cost_center, loss_a, loss_b, loss_c, valid_area) {
  # 生成成本、损失的随机数
  cost <- rtriangle(n_sim, a = 0.9*cost_center, b = 1.1*cost_center, c = cost_center)
  loss <- rtriangle(n_sim, a = loss_a, b = loss_b, c = loss_c)
  # 从预生成的有效面积中采样
  area <- sample(valid_area, n_sim, replace = TRUE)
  # 返回单建筑的模拟损失向量
  cost * area * loss
}

# 单一行的计算函数
calculate_row_percentile <- function(row, n_sim = 10000, building_params, valid_area) {
  total_loss <- numeric(n_sim)
  # 遍历每个建筑类型，累加对应数量的模拟结果
  for (bldg in names(row)) {
    count <- row[[bldg]]
    if (count == 0) next
    # 批量生成count个建筑的模拟结果，每行是一个建筑的n_sim次模拟
    bldg_sim <- replicate(count, simulate_building(n_sim, 
                                                   building_params[bldg, "cost_center"],
                                                   building_params[bldg, "loss_a"],
                                                   building_params[bldg, "loss_b"],
                                                   building_params[bldg, "loss_c"],
                                                   valid_area))
    # 按列求和（每个模拟场景下的该类型建筑总损失）
    total_loss <- total_loss + rowSums(bldg_sim)
  }
  # 返回95分位数
  quantile(total_loss, probs = 0.95)
}

第三步：并行处理所有行 #

# 设置并行计算，使用所有可用核心
cl <- makeCluster(detectCores())
registerDoParallel(cl)

# 并行遍历每一行计算
loss.PEISVII.95p <- foreach(i = 1:nrow(final.distrib), .combine = c) %dopar% {
  calculate_row_percentile(final.distrib[i, ], 
                           n_sim = 10000, 
                           building_params = building_params, 
                           valid_area = valid_area)
}

# 转换为数据框
loss.PEISVII.95p <- setNames(data.frame(loss.PEISVII.95p), "PEIS VII loss")

# 关闭并行集群
stopCluster(cl)

额外优化建议 #

• 调整预生成样本量：如果valid_area的数量足够（比如100万条），可以满足所有模拟的采样需求，无需重复生成
• 减少模拟次数：如果业务允许，可适当降低n_sim（比如从1万降到5千），能大幅缩短时间，同时对分位数结果影响有限
• 使用更快的随机数生成包：比如random包的随机数生成速度比基础包更快，可进一步提升效率

性能预期 #

经过上述优化，单批次计算效率可提升5-10倍，并行处理可再提升与CPU核心数相当的倍数。假设使用8核CPU，总耗时可从32小时缩短至约40分钟到1小时，完全可接受。

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