Fast (эр) способ матрицы индексации в R

Question

Прежде всего, я ищу быстрый (вдовцы) способ Подменит/индексирования матрицы много, много раз:

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

фон:
Я 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 дней. Использование параллельной обработки и более мощный процессор, что это можно сделать быстрее, но до сих пор не очень приятен:/

• Есть ли лучший способ ресэмплировать данные/избавиться от цикла?
• Может применяться, векторизовать, копировать и т.д.
• Имеет ли смысл реализовать подмножество в 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 выше.)