r - 如何在 R 中优化以下程序以提高性能？ (涉及计算密集型置换测试的蒙特卡罗模拟)

Question

您好，我正在尝试针对 Monte Carlo Power 研究优化以下程序。我首先使用一种取自 Efron 和 Tibshirani (1993) 的算法来找到均值相等的置换检验(上尾)的 p 值。我编写了一个名为 perm_test() 的函数，其输出是单个 p 值。然后，我在另一个程序 power_p() 中调用此函数，该程序模拟 1000 个排列测试(返回 1000 个 p 值)。我的功效估计是这 1000 个 p 值在统计上显着的比例，即 < 0.05。整个过程大约需要我运行 8 分钟(2020 MacBook pro)。我想知道是否有人在优化此过程以使其运行得更快方面有任何建议。很多thnx。

perm_test <- function(delta) {
  
  # Permutation test as described in Efron and Tibshirani (1993), (Algorithm 15.1, p208)
  
  # Draw random samples from a normal distribution
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  # observed diff in means, denoted as D_obs 
  D_obs <- mean(x) - mean(y)
  
  # Create a data frame "N" with  n_x + n_y obs (20 rows in our case)
  N <- data.frame("v" = c(x, y))
  # create a group variable "g" indicating which group each observation belongs to
  N$g <- as.factor(c(rep("x", 10), rep("y", 10)))
  # arrange column "v" in ascending order 
  # notice that column "g" is also re-ordered
  N <- arrange(N, v)
  
  ###############################################################################################
  # There are 20 choose 10 (184756) possibilities for the ordering of "g"                       #       
  # corresponding to all possible ways of partitioning 20 elements into two subsets of size 10  #
  # we take only a random sample of 5000 from those 184756 possibilities                        #
  ###############################################################################################
  
  # Initialize variables
  B <- 5000
  x_mean <- 0
  y_mean <- 0
  D_perm <- rep(0, 5000)
  
  # Loop to randomly generate 5000 different orderings of "g"
  for (i in 1:B) {
    
    # Shuffle the ordering of "g"
    N$g <- sample(N$g)
    # Permuted means of x and y
    x_mean <- tapply(N$v, N$g, mean)[1]
    y_mean <- tapply(N$v, N$g, mean)[2]
    # Find permuted diff in means, denoted as D_perm
    D_perm[i] <- x_mean - y_mean 
    }
  
  # Find p-value 
  P_perm <- sum(D_perm >= D_obs)/ B
  
  # Output
  return(round(P_perm, digits = 5))
  
}

这是模拟 1000 次排列测试的程序:

power_p <- function(numTrial, delta) {
  
  # Initilize variables
  P_p <- rep(0, numTrial) 
  pwr_p <- 0
  
  # Simulation
  P_p <- replicate(n = numTrial, expr = perm_test(delta))
  
  # Power estimates are the proportions of p-values that are significant (i.e. less than 0.05)
  pwr_p <- sum(P_p < 0.05) / numTrial
  
  # Output 
  return(round(pwr_p, digits = 5))

}

Answer 1

perm_test2 <- function(delta) {
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  D_obs <- mean(x) - mean(y)
  v <- c(x, y)
  g <- rep(1:2, each = 10)
  B <- 5000
  y_mean <- x_mean <- 0
  D_perm <- rep(0, B)
  for (i in 1:B) {
    ii <- sample(g) == 1L
    D_perm[i] <- (sum(v[ii]) - sum(v[!ii]))/10 
  }
  P_perm <- sum(D_perm >= D_obs)/ B
  return(round(P_perm, digits = 5))
}

与更复杂的方法相比，我提出了上述方法，其中我对您现有的代码进行了一些简单的改进。

1. 我们不需要使用data.frame，每次只需要花费不必要的时间来子集向量
2. 我们可以使用简单的整数向量来代替因子组向量
3. 不需要数据排序
4. 我们可以减少均值计算。在您的代码中，tapply(N$v, N$g, mean) 的 2 次调用是最慢的部分。
5. mean(x) 比 sum(x)/n 慢，它会做额外的检查等，所以在这种情况下我们可以使用更快的方法。因为内部循环将执行 1000x5000 次(sims x B)。

---

bench::mark(perm_test_org(0.5), perm_test2(0.5), iterations = 5, check = F,
            relative = T)[, 1:8]
# A tibble: 2 x 8
#   expression           min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#   <bch:expr>         <dbl>  <dbl>     <dbl>     <dbl>    <dbl> <int> <dbl>
# 1 perm_test_org(0.5)  15.5   15.3       1        1.00     1        5    16
# 2 perm_test2(0.5)      1      1        13.5      1        1.69     5     2

在我的系统上大约快 15 倍。 1000 次迭代耗时 33.89 秒。

更新 2。

我们可以通过以下方式进一步提高速度:

1. 将 sample 替换为 sample.int
2. 然后我们看到，根本不需要 g 向量来随机选择两个组
3. 在循环中，我们不需要对向量 v 的两个部分求和，因为我们可以在循环之前 sum(v)，因此我们可以在循环内少做一个求和最后计算结果值。

---

perm_test3 <- function(delta) {
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  D_obs <- mean(x) - mean(y)
  v <- c(x, y)
  B <- 5000
  s <- sum(v)
  D_perm2 <- rep(0, B)
  for (i in 1:B) {
    D_perm2[i] <- sum(v[sample.int(10) < 6])
  }
  D_perm <- D_perm2 - (s - D_perm2)
  P_perm <- sum(D_perm/10 >= D_obs) / B
  return(round(P_perm, digits = 5))
}

在 +/- 20 秒内运行 1000 次迭代。现在最慢的部分是重复调用 sample.int。您可以查看更快的功能: https://www.r-bloggers.com/2019/04/fast-sampling-support-in-dqrng/