Calcolo efficiente della deviazione standard cumulativa della matrice in

Question

Recentemente ho postato questa domanda sulla mailing list di r-help ma non ho ricevuto risposte, quindi ho pensato di postarla anche qui e vedere se c'erano suggerimenti.

Sto tentando di calcolare la deviazione standard cumulativa di una matrice. Voglio una funzione che accetta una matrice e restituisce una matrice della stessa dimensione in cui la cella di output (i, j) è impostata sulla deviazione standard della colonna di input j tra le righe 1 e i. Le NA dovrebbero essere ignorate, a meno che la cellula (i, j) della matrice di input stessa sia NA, nel qual caso anche la cella (i, j) della matrice di output dovrebbe essere NA.

Impossibile trovare una funzione incorporata, quindi ho implementato il seguente codice. Sfortunatamente, questo usa un ciclo che finisce per essere un po 'lento per le matrici di grandi dimensioni. C'è una funzione incorporata più veloce o qualcuno può suggerire un approccio migliore?

cumsd <- function(mat) 
{ 
    retval <- mat*NA 
    for (i in 2:nrow(mat)) retval[i,] <- sd(mat[1:i,], na.rm=T) 
    retval[is.na(mat)] <- NA 
    retval 
} 

Grazie.

Answer 1

Si potrebbe utilizzare cumsum per calcolare le somme necessarie da formule dirette per varianza/sd alle operazioni vectorized su matrice:

cumsd_mod <- function(mat) { 
    cum_var <- function(x) { 
     ind_na <- !is.na(x) 
     nn <- cumsum(ind_na) 
     x[!ind_na] <- 0 
     cumsum(x^2)/(nn-1) - (cumsum(x))^2/(nn-1)/nn 
    } 
    v <- sqrt(apply(mat,2,cum_var)) 
    v[is.na(mat) | is.infinite(v)] <- NA 
    v 
} 

solo per il confronto:

set.seed(2765374) 
X <- matrix(rnorm(1000),100,10) 
X[cbind(1:10,1:10)] <- NA # to have some NA's 

all.equal(cumsd(X),cumsd_mod(X)) 
# [1] TRUE 

E una questione di tempi:

X <- matrix(rnorm(100000),1000,100) 
system.time(cumsd(X)) 
# user system elapsed 
# 7.94 0.00 7.97 
system.time(cumsd_mod(X)) 
# user system elapsed 
# 0.03 0.00 0.03 

Answer 2

Un'altra prova (Marek è più veloce)

cumsd2 <- function(y) { 
n <- nrow(y) 
apply(y,2,function(i) { 
    Xmeans <- lapply(1:n,function(z) rep(sum(i[1:z])/z,z)) 
    Xs <- sapply(1:n, function(z) i[1:z]) 
    sapply(2:n,function(z) sqrt(sum((Xs[[z]]-Xmeans[[z]])^2,na.rm = T)/(z-1))) 
}) 
}