Эффективный расчет матрицы совокупного стандартного отклонения в r

Question

Я недавно разместил этот вопрос в списке рассылки r-help, но не получил ответов, поэтому я решил опубликовать его и здесь и посмотреть, есть ли какие-либо предложения.

Я пытаюсь вычислить совокупное стандартное отклонение матрицы. Я хочу функцию, которая принимает матрицу и возвращает матрицу того же размера, где выходная ячейка (i, j) установлена ​​на стандартное отклонение входного столбца j между строками 1 и i. NA следует игнорировать, если только ячейка (i, j) входной матрицы не является NA, и в этом случае ячейка (i, j) выходной матрицы также должна быть NA.

Мне не удалось найти встроенную функцию, поэтому я реализовал следующий код. К сожалению, для этого используется цикл, который оказывается несколько медленным для больших матриц. Есть ли более быстрая встроенная функция или кто-то может предложить лучший подход?

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
}

Спасибо.

Answer 1

Вы можете использовать cumsum для вычисления необходимых сумм от прямых формул для дисперсии / sd до векторизованных операций над матрицей:

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
}

просто для сравнения:

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

А по поводу времени:

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

Еще одна попытка (Марек быстрее)

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