用代码优化加速R中的仿真

Question

我试图做的一般版本是进行一项模拟研究，其中我操作了几个变量，看看它是如何影响结果的。我对R有一些速度问题，最近的模拟只进行了几次迭代(每次实验10次)。然而，当我移动到一个大规模的版本(每个实验10K)时，模拟已经运行了14个小时(并且仍然在运行)。

下面是我正在运行的代码(带注释)。作为一个新手与R，并正在努力优化的模拟是有效的。我希望从这里提供的评论和建议中学习优化代码，并将这些注释用于未来的仿真研究。

让我就这段代码应该做什么说几句话。我正在操作两个变量:效果大小和样本大小。每个组合运行10k次(即每个条件下的10k实验)。我初始化一个数据框架来存储我的结果(称为结果)。I循环三个变量:效果大小、样本大小和迭代(10k)。

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

我从控制条件(DVcontrol)的标准正常值(DVexperiment)中生成随机数据，并使用我的影响大小作为实验条件(DVexperiment)的平均值。我取值之间的差值，将结果放入R(配对样本t检验)中的t检验函数中。我将结果存储在一个名为Trials的列表中，并将其重新绑定到结果数据框架。这个过程重复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)
  }
}

谢谢你的意见和建议，乔什

Answer 1

优化的一般方法是运行一个轮廓仪，以确定解释器花在哪一部分代码中的时间最多，然后优化该部分。假设您的代码驻留在一个名为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

从这里，我们观察到rnorm和t.test是您最昂贵的操作(不应该是一个惊喜，因为这些都在您的内部-最循环中)。

一旦确定了昂贵的函数调用在哪里，实际的优化由两个步骤组成：

1. 优化函数和/或
2. 优化函数调用的次数。

由于t.test和rnorm是内置的R函数，所以以上步骤1的唯一选项是寻找替代包，这些包可以更快地实现来自正态分布的采样和/或运行多个t测试。第二步实际上是以一种不会多次重新计算相同事物的方式来重构您的代码。例如，以下代码行不依赖于i

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

将这些数据移出循环是否有意义，还是真的需要为每个不同的i值提供测试数据的新示例？