提高数据帧模拟的性能

Question

在编写代码时，我经常只是编写我想到的代码。虽然我认为我从一开始就学习了高效的 R 编码（例如，避免 for ... if 循环），但我的解决方案并不总是由性能驱动。不幸的是，有时知道什么是最有效的代码是至关重要的——我想学习这一点！

目前我正在模拟多个数据帧组合成一个列表。模拟之后，我需要第二个数据框，其中包含整个列表中所有列的平均值和 SD。 （这里的“模拟”表示某些变量是从其他数据帧中模拟/重新采样的，其他变量只是具有特定 b_0 的随机正态或二项分布值。为简洁起见，我省略了第一部分在这里重新采样。）

我的代码（参见下面的示例）完美地产生了预期的结果，但它似乎首先有点慢（我说的是真实的几个小时），其次是高度消耗 RAM（为此我暂时减少了列表中模拟dfs的数量）。

对于模拟，我知道在函数中定义 data.frame 可能是一个问题，但我不知道如何做得更好。对于均值/SD 数据帧，我只能说它更慢。

如何提高代码的性能？谁能提供一些关于这种性能提升的基本规则（或相关信息来源）？

（我正在使用 R 3.x/64 和 Win 7/64 AMD FX(tm)-8350 八核处理器，4 GHz，16 GB 机器。CPU 在运行时保持相当凉爽，RAM在它的极限中呻吟。）

这里我给出一个用 cmets 测量系统时间的示例代码：

# definitions
r <- 1e5 # number of rows
n <- 1e3 # number of dfs

# simulation of the list  
library(dplyr)
system.time(list <- lapply(1:n, function(i){       # 59.05 sec
  data.frame(a = rbinom(r, 1, .375)) %>%
    mutate(
      b = rnorm(r, 0, 2),
      c = .42 * rnorm(r, 0, 6),
      d = rbinom(r, 11, c(1:11)/11),
      e = rbinom(r, 1, .1),
      f = .02 * rnorm(r, 0, 5))
}))

# df w/ means & sds
system.time(list.s <- data.frame(                  # 73.20 sec
  list.mean = round(rowMeans(sapply(list, colMeans)), 2),
  list.sd = round(sapply(do.call(rbind, list), sd), 2)))

Answer 1

扩展 Rolands 评论，您可以预先创建大量人口数据，然后为每个“样本”/迭代简单地对其进行子集化。示例：

## create large population data:

s <- 1e6 # probably big enough for this problem
set.seed(12)
d <- matrix(NA, nrow = s, ncol = 6) #..
# using matrix is more efficient than data.frame
d[,1] <- rbinom(s, 1, .375)
d[,2] <- rnorm(s, 0, 2)
d[,3] <- .42 * rnorm(s, 0, 6)
d[,4] <- rbinom(s, 11, c(1:11)/11)
d[,5] <- rbinom(s, 1, .1)
d[,6] <- .02 * rnorm(s, 0, 5)
head(d)
#      [,1]        [,2]      [,3] [,4] [,5]        [,6]
# [1,]    0  0.73853351  1.097805    1    0 -0.06233008
# [2,]    1 -0.05311206  4.447807    2    0 -0.01117972
# [3,]    1  1.71576276 -3.619708    6    0  0.02962562
# [4,]    0  1.92188205 -1.062585    2    0  0.03195146
# [5,]    0 -1.41097404  1.706067    2    0 -0.07751285
# [6,]    0  4.19130890  2.663374    8    0 -0.02316172


r <- 1e4 # number of rows
n <- 1e2 # number of dfs

si <- replicate(n, sample.int(s, r)) # get indexes for each sample 

# loop trougth samples and subset data:
nSamples <- lapply(1:n, function(x) {
  d[si[, x],]
  })

# and calculate colMeans:
list.mean2 = round(rowMeans(sapply(nSamples, colMeans)), 3)
list.mean2
# [1]  0.376  0.000 -0.003  5.999  0.100  0.000

与你的结果比较：

require(dplyr)
list1 <- lapply(1:n, function(i){
  data.frame(a = rbinom(r, 1, .375)) %>%
    mutate(
      b = rnorm(r, 0, 2),
      c = .42 * rnorm(r, 0, 6),
      d = rbinom(r, 11, c(1:11)/11),
      e = rbinom(r, 1, .1),
      f = .02 * rnorm(r, 0, 5))
})

list.mean1 = round(rowMeans(sapply(list1, colMeans)), 3)
list.mean1
# a      b      c      d      e      f 
# 0.375 -0.002  0.004  6.001  0.100  0.000 

我们可以看到，平均值的估计值与这个小的 n 值非常相似。

附：由于“列表”是基本 R 函数，因此您不应使用该名称命名变量！

让我们将这两种方法包装成函数来测试时间：

mySim <- function(s, r, n) {
  d <- matrix(NA, nrow = s, ncol = 6)
  d[,1] <- rbinom(s, 1, .375)
  d[,2] <- rnorm(s, 0, 2)
  d[,3] <- .42 * rnorm(s, 0, 6)
  d[,4] <- rbinom(s, 11, c(1:11)/11)
  d[,5] <- rbinom(s, 1, .1)
  d[,6] <- .02 * rnorm(s, 0, 5)
  si <- lapply(1:n, function(x) sample.int(s, r))
  nSamples <- lapply(si, function(x) {
    d[x,]
  })
  list.mean2 = rowMeans(sapply(nSamples, colMeans))
  list.mean2
}

yourSim <- function(r, n) {
  require(dplyr)
  list1 <- lapply(1:n, function(i){
    data.frame(a = rbinom(r, 1, .375)) %>%
      mutate(
        b = rnorm(r, 0, 2),
        c = .42 * rnorm(r, 0, 6),
        d = rbinom(r, 11, c(1:11)/11),
        e = rbinom(r, 1, .1),
        f = .02 * rnorm(r, 0, 5))
  })
  list.mean1 = rowMeans(sapply(list1, colMeans))
  list.mean1
}

system.time(mySim(1e6, 1e4, 1e2)) # ~ 0.6 sek
system.time(yourSim(1e4, 1e2)) # ~ 1.5 sek

# if s = 1e7 :
system.time(mySim(1e7, 1e4, 1e2)) # ~ 4.53 sek

我们可以看到，为较小的 n 和 r 值创建大量人口数据并没有提高速度。

让s 为 1e6（100 万），但你应该自己调查一下 足够了。

如果我们为较大的 'r' 和 'n' 值计算时间：

system.time(r1 <- mySim(1e6, 1e5, 1e3)) # ~ 20 sek
system.time(r2 <- yourSim(1e5, 1e3)) # ~ 60 sek

round(r1, 3)
# [1]  0.376 -0.003 -0.002  6.001  0.100  0.00
round(r2, 3)
# a     b     c     d     e     f 
# 0.375 0.000 0.000 6.000 0.100 0.000 

关于计算 SD： 也许您想使用包 'matrixStats' 中的 'rowSds()' 或 'colSds()'？

另外我建议你研究一下 Rcpp 包，它可能有助于进一步加快代码速度。