在R语言中高效计算矩阵累积标准差

Question

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

我正在尝试计算矩阵的累积标准差。我想要一个接受矩阵并返回相同大小的矩阵的函数，其中输出单元格(i,j)设置为在第1行和第i行之间输入列j的标准差。除非输入矩阵本身的单元格(i,j)是NA，否则应忽略NAs，在这种情况下，输出矩阵的单元格(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
}

感谢您。

- Abiel

Answer 1

你可以使用cumsum函数，通过直接公式计算方差/标准差的必要和来执行矩阵化操作：

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 

- Marek

---

非常好的Marek，这使得我的分析更加高效。顺便说一下，在函数中似乎没有使用变量n <- nrow(mat)。 - Abiel

这是早期版本的残留物 ;)。 - Marek

注意使用这个算法；@Marek 提出的想法很好，但是当标准差与平均值相比较小时，使用这个方差公式可能会得到有趣的结果。维基百科有更好的算法；还可以看看我在这里的回答。 - Aaron left Stack Overflow

@Aaron 这是非常好的观点。在空闲时间里，我会更新我的答案。谢谢。 - Marek

Answer 2

再试一次（Marek的更快）

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

- gd047