veloce (er) modo di matrice indicizzazione in R

Question

Foremost, Cerco un modo veloce (er) di sottoinsiemi/indicizzazione di una matrice di molte, molte volte:

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

Background:
I' m implementando una procedura di test sequenziale che coinvolge il bootstrap in R. Volendo replicare alcuni risultati di simulazione, mi sono imbattuto in questo collo di bottiglia in cui è necessario fare molto indicizzazione. Per l'implementazione del blocco-bootstrap ho creato una matrice di indici con cui ho suddiviso la matrice di dati originale per disegnare i campioni dei dati.

# 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' 

L'esecuzione della procedura di test sequenziale richiede circa 10 secondi. Usando questo in simulazioni con 2500 repliche e diverse costellazioni di parametri , ci vorrebbe qualcosa come 40 giorni. Utilizzando l'elaborazione parallela e una migliore potenza della CPU è possibile fare più veloce, ma ancora non molto gradevole:/

• Esiste un modo migliore per ricampionare i dati/sbarazzarsi del ciclo?
• È possibile applicare, vettorizzare, replicare ecc. Venire in qualsiasi luogo?
• Sarebbe logico implementare il sottotitolo in C (ad esempio, manipolare alcuni puntatori)?

Anche se ogni singolo passaggio è già stato eseguito in modo incredibilmente veloce da R, non è abbastanza veloce.
Sarei davvero molto felice per qualsiasi tipo di risposta/aiuto/consiglio!

Qs correlati:
- 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

da lì

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

in realtà non lo fanno per me.

Answer 1

ho riscritto makeStats e makeIndex come erano due dei più grandi colli di bottiglia:

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) 
} 

Questo ha portato giù i tempi di calcolo da 28 a 6 secondi sulla mia macchina. Scommetto che ci sono altre parti del codice che possono essere migliorate (incluso il mio uso di lapply/tabulate sopra.)