Cálculo eficiente de la desviación estándar acumulada de la matriz en r

Question

Recientemente publiqué esta pregunta en la lista de correo de r-help pero no obtuve ninguna respuesta, así que pensé en publicarla aquí también y ver si había alguna sugerencia.

Estoy tratando de calcular la desviación estándar acumulativa de una matriz. Quiero una función que acepte una matriz y devuelva una matriz del mismo tamaño donde la celda de salida (i, j) se establece en la desviación estándar de la columna de entrada j entre las filas 1 ei. Las NA deben ignorarse, a menos que la celda (i, j) de la matriz de entrada en sí sea NA, en cuyo caso la celda (i, j) de la matriz de salida también debe ser NA.

No pude encontrar una función incorporada, así que implementé el siguiente código. Desafortunadamente, esto utiliza un bucle que termina siendo algo lento para matrices grandes. ¿Hay una función incorporada más rápida o alguien puede sugerir un mejor enfoque?

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 
} 

Gracias.

Answer 1

Usted podría utilizar cumsum para calcular las sumas necesarias de fórmulas directas para la varianza/SD para operaciones vectorizadas de matriz:

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 
} 

sólo para comparación:

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 

Y sobre el calendario:

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

Otro intento (Marek es más rápido)

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