如何实现R语言中百万次模拟循环的最快运行速度？

我正在开展一项指数分布的模拟研究，已知其分布函数为 $F(x) = 1 - e^{-λx}$，由此推导出生成样本的公式为 $x = (-ln[1-F(x)])/(λ)$，对应R语言代码为：

x = (-log(1-p))/(lambda)

其中$F(x)$对应代码中的p。该分布的对数似然函数为：
$L_X = n_xlog(lambda) - lambdasum(x)$

由于计划运行多达100万次迭代，当前的实现速度无法满足需求：基础循环方案耗时近10小时，并行方案约3.4小时，急需进一步提升速度与效率。我的核心目标是收集所有迭代生成的“模拟lambda”值。

现有两种实现方案 #

(1) 基础循环方案 #

iterations <- 100
lambda = 35 #par[1]

simulation <- function(i) {
  
  p <- runif(1000, 0, 1)
  x <- (-log(1-p))/(lambda)
  n_x <- length(x)
  
  loglikelihood = function(par, x, y){
    return( -(n_x*log(par[1]) - par[1]*sum(x)))
  }
  
  result <- optim(c(10), fn = loglikelihood, x = x , method = "Brent",
                  lower = 0, upper = 100)
  
  return(result$par)
}

results_list <- lapply(1:iterations, simulation)                          
results <- do.call(rbind, results_list)
colnames(results) <- c("lambda_sim")
results

(2) 并行方案 #

future::plan(future::multisession, workers = parallelly::availableCores() - 1)

iterations <- 100
lambda = 35 #par[1]

simulation <- function(i) {
  
  p <- runif(1000, 0, 1)
  x <- (-log(1-p))/(lambda)
  n_x <- length(x)
  
  loglikelihood = function(par, x, y){
    return( -(n_x*log(par[1]) - par[1]*sum(x)))
  }
  
  result <- optim(c(10), fn = loglikelihood, x = x , method = "Brent",
                  lower = 0, upper = 100)   

  return(result$par)
}

results_list <- future.apply::future_lapply(1:iterations, simulation, future.seed = TRUE) 
results <- do.call(rbind, results_list)
colnames(results) <- c("lambda_sim")

future::plan(future::sequential)

疑问与后续计划 #

• 不确定clusterApply()是否是当前R语言中性能最优的并行实现方法
• 关于optim()初始值的确定，我将在新帖中单独提问

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