通過代碼優化加速R中的仿真

Question

我試圖做的通用版本是進行仿真研究，在該仿真研究中我操縱一些變量以查看結果如何影響結果。我遇到了一些與速度有關的問題。最新的仿真工作只需幾次迭代（每個實驗10次）。但是，當我進行大規模（每個實驗10k）版本時，仿真已經運行了14個小時（並且仍在運行）。

下面是我正在運行的代碼（帶有註釋）。作爲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的列表中，並將其與結果數據框綁定。這個過程重複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

以優化的一般方法是運行profiler，以確定哪些部分代碼的解釋花費的時間最多，然後優化該部分。假設您的代碼位於名爲test.R的文件中。在R，您可以通過運行下面的命令序列圖譜是：（請注意，這些命令將在您的R會話目錄下生成一個文件Rprof.out）

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

如果我們運行分析器在你的代碼（與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檢驗替代包。第2步實際上是以不重新計算同一事物的方式重構代碼。例如，下面的代碼行不依賴於i：

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

是否有意義移動這些外循環，或者你真的需要你的測試數據的新樣本的i每個不同的值？