r中矩阵累积标准差的有效计算

Question

我最近在r-help邮件列表上发布了这个问题,但没有得到答案,所以我想我也会在这里发布,看看是否有任何建议.

我试图计算矩阵的累积标准差.我想要一个接受矩阵并返回相同大小的矩阵的函数,其中输出单元格(i,j)被设置为行1和i之间的输入列j的标准偏差.应忽略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