通过代码优化加速 R 中的仿真

Speed up Simulation in R with Code Optimization

我正在尝试做的通用版本是进行模拟研究,我在其中操纵一些变量以查看它如何影响结果。我在使用 R 时遇到了一些速度问题。最新的模拟进行了几次迭代(每个实验 10 次)。但是,当我转到大规模(每个实验 10k)版本时,模拟已经 运行ning 了 14 个小时(现在仍然是 运行ning)。

下面是我 运行ning 的代码(带注释)。作为 R 的新手,我正在努力优化模拟以提高效率。我希望能从此处提供的意见和建议中学习,以优化此代码,并将这些意见用于未来的模拟研究。

让我谈谈这段代码应该做什么。我正在操纵两个变量:效果大小和样本大小。每个组合是 运行 10k 次(即每个条件 10k 次实验)。我初始化一个数据框来存储我的结果(称为结果)。我遍历三个变量:效果大小、样本大小和迭代次数 (10k)。

在循环中,我初始化了四个 NULL 组件:p.test、p.rep、d.test 和 d.rep。前两个捕获初始 t 检验的 p 值和复制的 p 值(在类似条件下复制)。后两者计算效果大小(Cohen's d)。

我从控制条件 (DVcontrol) 的标准正态生成随机数据,并使用我的效应量作为实验条件 (DVexperiment) 的平均值。我取值之间的差异并将结果放入 R 中的 t 检验函数(配对样本 t 检验)。我将结果存储在一个名为 Trials 的列表中,并将其绑定到 Results 数据框。这个过程重复 10k 次直到完成。

# Set Simulation Parameters
## Effect Sizes (ES is equal to mean difference when SD equals Variance equals 1)
effect_size_range <- seq(0, 2, .1) ## ES
## Sample Sizes
sample_size_range <- seq(10, 1000, 10) ## SS
## Iterations for each ES-SS Combination
iter <- 10000

# Initialize the Vector of Results
Results <- data.frame()

# Set Random Seed
set.seed(12)

# Loop over the Different ESs
for(ES in effect_size_range) {

  # Loop over the Different Sample Sizes
  for(SS in sample_size_range) {

    # Create p-value Vectors
    p.test <- NULL
    p.rep <- NULL
    d.test <- NULL
    d.rep <- NULL

    # Loop over the iterations
    for(i in 1:iter) {

      # Generate Test Data
      DVcontrol <- rnorm(SS, mean=0, sd=1)
      DVexperiment <- rnorm(SS, mean=ES, sd=1)
      DVdiff <- DVexperiment - DVcontrol
      p.test[i] <- t.test(DVdiff, alternative="greater")$p.value
      d.test[i] <- mean(DVdiff) / sd(DVdiff)

      # Generate Replication Data
      DVcontrol <- rnorm(iter, mean=0, sd=1)
      DVexperiment <- rnorm(iter, mean=ES, sd=1)
      DVdiff <- DVexperiment - DVcontrol
      p.rep[i] <- t.test(DVdiff, alternative="greater")$p.value
      d.rep[i] <- mean(DVdiff) / sd(DVdiff)
    }

    # Results
    Trial <- list(ES=ES, SS=SS,
                  d.test=mean(d.test), d.rep=mean(d.rep),
                  p.test=mean(p.test), p.rep=mean(p.rep),
                  r=cor(p.test, p.rep, method="kendall"),
                  r.log=cor(log2(p.test)*(-1), log2(p.rep)*(-1), method= "kendall"))
    Results <- rbind(Results, Trial)
  }
}

提前感谢您的意见和建议, 乔什

一般的优化方法是运行一个profiler来确定解释器花费最多时间的代码部分,然后优化该部分。假设您的代码位于名为 test.R 的文件中。在 R 中,您可以通过 运行 执行以下命令序列来分析它:

Rprof()              ## Start the profiler
source( "test.R" )   ## Run the code
Rprof( NULL )        ## Stop the profiler
summaryRprof()       ## Display the results

(请注意,这些命令将在您的 R 会话目录中生成一个文件 Rprof.out。)

如果我们 运行 您代码上的分析器(使用 iter <- 10,而不是 iter <- 10000),我们会得到以下分析:

# $by.self
#                         self.time self.pct total.time total.pct
# "rnorm"                      1.56    24.53       1.56     24.53
# "t.test.default"             0.66    10.38       2.74     43.08
# "stopifnot"                  0.32     5.03       0.86     13.52
# "rbind"                      0.32     5.03       0.52      8.18
# "pmatch"                     0.30     4.72       0.34      5.35
# "mean"                       0.26     4.09       0.42      6.60
# "var"                        0.24     3.77       1.38     21.70

从这里,我们观察到 rnormt.test 是您最昂贵的操作(真的不应该感到惊讶,因为它们在您的最内层循环中)。

一旦你弄清楚了昂贵的函数调用在哪里,实际的优化包括两个步骤:

  1. 优化功能,and/or
  2. 优化函数调用次数

由于 t.testrnorm 是内置的 R 函数,您对上述步骤 1 的唯一选择是寻找替代包,这些包可能具有更快的正态分布采样实现 and/or 运行宁多 t 检验。第 2 步实际上是关于以一种不会多次重新计算同一事物的方式重构您的代码。例如,以下代码行不依赖于 i:

# Generate Test Data
DVcontrol <- rnorm(SS, mean=0, sd=1)
DVexperiment <- rnorm(SS, mean=ES, sd=1)

将这些移出循环是否有意义,或者您真的需要为每个不同的 i 值获取新的测试数据样本吗?