快速（ER）的方式

Question

最重要的，我在尋找子集的快速（ER）的方式/索引矩陣很多很多次m實現了一個涉及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功率是可以做到更快，但 仍然不是很討好：/

• 有沒有更好的方式來重新採樣數據/擺脫循環？
• 可以申請，Vectorize，複製等進來嗎？
• 在C中實現子集（例如操作某些指針）有意義嗎？

即使每一步已經被R完成得非常快，但它不夠快。
對於任何形式的回覆/幫助/建議，我會非常高興！

相關Qs的：
- Fast matrix subsetting via '[': by rows, by columns or doesn't matter?
- fast function for generating bootstrap samples in matrix forms in R
- random sampling - matrix

從那裏

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

並沒有真正爲我做的。

Answer 1

我改寫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）。