r的矩陣累積標準偏差的有效計算

我最近在r-help郵件列表上發佈了這個問題，但沒有得到答案，所以我想我會在這裏發佈它，看看是否有任何建議。r的矩陣累積標準偏差的有效計算

我想計算一個矩陣的累積標準偏差。我想要一個接受矩陣的函數，並返回一個相同大小的矩陣，其中輸出單元格（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 
}

謝謝。

來源

2010-05-04 Abiel

你可以使用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

來源

2010-05-04 14:13:41 Marek

非常好的馬立克，這使得我的分析更有效率。僅供參考，它看起來不像你在函數中使用變量n < - nrow（mat）。 – Abiel 2010-05-04 15:07:36

這是來自早期版本之一的殘留物;）。 – Marek 2010-05-04 18:57:11

小心使用這種算法; @Marek有一個好主意，但是當sd相對於均值小時，使用這個方差來表示方差可以給出有趣的結果。維基百科有[更好的算法]（http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance）;也看到我的答案[這裏]（http://stackoverflow.com/questions/7474943/surprisingly-slow-standard-deviation-in-r/7475664#7475664）。 – Aaron 2011-09-19 19:05:01

另一個嘗試（馬立克更快）

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

來源

2010-05-04 15:19:34

r的矩陣累積標準偏差的有效計算

回答

相關問題