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'

一連のテスト手順を 1 回実行するのに約 10 秒かかります。これを 2500 回の複製といくつかのパラメーター コンスタレーションを使用したシミュレーションで使用すると、40 日ほどかかることになります。並列処理とより優れた CPU パワーを使用すると、より高速に実行できますが、それでもあまり満足のいくものではありません :/

  • データを再サンプリングする/ループを取り除くより良い方法はありますか?
  • 適用、ベクトル化、複製などはどこにでも入ることができますか?
  • C でサブセット化を実装する (たとえば、いくつかのポインターを操作する) ことは理にかなっていますか?

すべてのステップは R によってすでに信じられないほど高速に実行されていますが、十分な速さではありません。
どんな種類の応答/ヘルプ/アドバイスも本当にうれしいです!

関連する Q:
- '[': 行ごと、列ごと、または関係ないことによる高速なマトリックスのサブセット化?
- R の行列形式でブートストラップ サンプルを生成するための高速関数
-ランダム サンプリング - 行列

そこから

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

本当に私のためにそれをしませんでした。

4

1 に答える 1