r中矩阵累积标准差的高效计算 [英] Efficient calculation of matrix cumulative standard deviation in r
问题描述
我最近在 r-help 邮件列表上发布了这个问题,但没有得到答案,所以我想我也将它发布在这里,看看是否有任何建议.
I recently posted this question on the r-help mailing list but got no answers, so I thought I would post it here as well and see if there were any suggestions.
我正在尝试计算矩阵的累积标准偏差.我想要一个函数,它接受一个矩阵并返回一个相同大小的矩阵,其中输出单元格 (i,j) 设置为第 1 行和第 i 行之间输入列 j 的标准偏差.NA 应该被忽略,除非输入矩阵的单元格 (i,j) 本身就是 NA,在这种情况下,输出矩阵的单元格 (i,j) 也应该是 NA.
I am trying to calculate the cumulative standard deviation of a matrix. I want a function that accepts a matrix and returns a matrix of the same size where output cell (i,j) is set to the standard deviation of input column j between rows 1 and i. NAs should be ignored, unless cell (i,j) of the input matrix itself is NA, in which case cell (i,j) of the output matrix should also be NA.
我找不到内置函数,所以我实现了以下代码.不幸的是,这使用了一个循环,最终对于大型矩阵来说有点慢.有没有更快的内置函数,或者有人可以提出更好的方法吗?
I could not find a built-in function, so I implemented the following code. Unfortunately, this uses a loop that ends up being somewhat slow for large matrices. Is there a faster built-in function or can someone suggest a better approach?
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
}
谢谢.
推荐答案
您可以使用 cumsum
来计算从方差/sd 的直接公式到矩阵矢量化运算所需的总和:
You could use cumsum
to compute necessary sums from direct formulas for variance/sd to vectorized operations on matrix:
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
这篇关于r中矩阵累积标准差的高效计算的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!