7

最重要的是，我正在寻找一种快速（更好）的方法来多次对矩阵进行子集/索引：

for (i in 1:99000) {
  subset.data <- data[index[, i], ]
}

背景：
我正在实施一个涉及 R 中引导程序的顺序测试过程。想要复制一些模拟结果，我遇到了需要进行大量索引的瓶颈。为了实现块引导，我创建了一个索引矩阵，我用它对原始数据矩阵进行子集化以绘制数据的重新采样。

# The basic setup

B <- 1000 # no. of bootstrap replications
n <- 250  # no. of observations
m <- 100  # no. of models/data series

# Create index matrix with B columns and n rows.
# Each column represents a resampling of the data.
# (actually block resamples, but doesn't matter here).

boot.index <- matrix(sample(1:n, n * B, replace=T), nrow=n, ncol=B)

# Make matrix with m data series of length n.

sample.data <- matrix(rnorm(n * m), nrow=n, ncol=m)

subsetMatrix <- function(data, index) { # fn definition for timing
  subset.data <- data[index, ]
  return(subset.data)
}

# check how long it takes.

Rprof("subsetMatrix.out")
for (i in 1:(m - 1)) { 
  for (b in 1:B) {  # B * (m - 1) = 1000 * 99 = 99000
    boot.data <- subsetMatrix(sample.data, boot.index[, b])
    # do some other stuff
  }
  # do some more stuff
}
Rprof()
summaryRprof("subsetMatrix.out")

# > summaryRprof("subsetMatrix.out")
# $by.self
#              self.time self.pct total.time total.pct
# subsetMatrix      9.96      100       9.96       100

# In the actual application:
#########
# > summaryRprof("seq_testing.out")
# $by.self
#              self.time self.pct total.time total.pct
# subsetMatrix       6.78    53.98       6.78     53.98
# colMeans           1.98    15.76       2.20     17.52
# makeIndex          1.08     8.60       2.12     16.88
# makeStats          0.66     5.25       9.66     76.91
# runif              0.60     4.78       0.72      5.73
# apply              0.30     2.39       0.42      3.34
# is.data.frame      0.22     1.75       0.22      1.75
# ceiling            0.18     1.43       0.18      1.43
# aperm.default      0.14     1.11       0.14      1.11
# array              0.12     0.96       0.12      0.96
# estimateMCS        0.10     0.80      12.56    100.00
# as.vector          0.10     0.80       0.10      0.80
# matrix             0.08     0.64       0.08      0.64
# lapply             0.06     0.48       0.06      0.48
# /                  0.04     0.32       0.04      0.32
# :                  0.04     0.32       0.04      0.32
# rowSums            0.04     0.32       0.04      0.32
# -                  0.02     0.16       0.02      0.16
# >                  0.02     0.16       0.02      0.16
#
# $by.total
#              total.time total.pct self.time self.pct
# estimateMCS        12.56    100.00      0.10     0.80
# makeStats           9.66     76.91      0.66     5.25
# subsetMatrix        6.78     53.98      6.78    53.98
# colMeans            2.20     17.52      1.98    15.76
# makeIndex           2.12     16.88      1.08     8.60
# runif               0.72      5.73      0.60     4.78
# doTest              0.68      5.41      0.00     0.00
# apply               0.42      3.34      0.30     2.39
# aperm               0.26      2.07      0.00     0.00
# is.data.frame       0.22      1.75      0.22     1.75
# sweep               0.20      1.59      0.00     0.00
# ceiling             0.18      1.43      0.18     1.43
# aperm.default       0.14      1.11      0.14     1.11
# array               0.12      0.96      0.12     0.96
# as.vector           0.10      0.80      0.10     0.80
# matrix              0.08      0.64      0.08     0.64
# lapply              0.06      0.48      0.06     0.48
# unlist              0.06      0.48      0.00     0.00
# /                   0.04      0.32      0.04     0.32
# :                   0.04      0.32      0.04     0.32
# rowSums             0.04      0.32      0.04     0.32
# -                   0.02      0.16      0.02     0.16
# >                   0.02      0.16      0.02     0.16
# mean                0.02      0.16      0.00     0.00
#
# $sample.interval
# [1] 0.02
#
# $sampling.time
# [1] 12.56'

执行一次顺序测试过程大约需要 10 秒。在具有 2500 次复制和几个参数星座的模拟中使用它，大约需要 40 天。使用并行处理和更好的 CPU 能力可以做得更快，但仍然不是很令人愉快：/

  • 有没有更好的方法来重新采样数据/摆脱循环？
  • 可以在任何地方应用、矢量化、复制等吗？
  • 在 C 中实现子集是否有意义（例如操作一些指针）？

尽管 R 已经完成了每一步都非常快，但它还不够快。
对于任何形式的回应/帮助/建议，我确实非常高兴！

相关问题：
-通过“[”进行快速矩阵子集化：按行、按列还是无关紧要？
-用于在 R 中以矩阵形式生成引导样本的快速函数
-随机抽样 - 矩阵

从那里

mapply(function(row) return(sample.data[row,]), row = boot.index)
replicate(B, apply(sample.data, 2, sample, replace = TRUE))

并没有真正为我做。

4

1 回答 1

3

我重写了makeStats，makeIndex因为它们是两个最大的瓶颈：

makeStats <- function(data, index) {

  data.mean <- colMeans(data)
  m <- nrow(data)
  n <- ncol(index)
  tabs <- lapply(1L:n, function(j)tabulate(index[, j], nbins = m))
  weights <- matrix(unlist(tabs), m, n) * (1 / nrow(index))
  boot.data.mean <- t(data) %*% weights - data.mean

  return(list(data.mean = data.mean,
              boot.data.mean = boot.data.mean))
}

makeIndex <- function(B, blocks){

  n <- ncol(blocks)
  l <- nrow(blocks)
  z <- ceiling(n/l)
  start.points <- sample.int(n, z * B, replace = TRUE)
  index <- blocks[, start.points]
  keep <- c(rep(TRUE, n), rep(FALSE, z*l - n))
  boot.index <- matrix(as.vector(index)[keep],
                       nrow = n, ncol = B)

  return(boot.index)
}

这将我机器上的计算时间从 28 秒缩短到 6 秒。我敢打赌，代码的其他部分可以改进（包括我上面使用的 lapply/tabulate。）

于 2013-12-08T23:56:14.730 回答