Hiệu quả tính toán độ lệch chuẩn tích lũy ma trận trong r

Question

Gần đây tôi đã đăng câu hỏi này trên danh sách gửi thư r-help nhưng không có câu trả lời, vì vậy tôi nghĩ tôi cũng sẽ đăng nó ở đây và xem có đề xuất nào không.

Tôi đang cố tính toán độ lệch chuẩn tích lũy của ma trận. Tôi muốn một hàm chấp nhận ma trận và trả về ma trận có cùng kích thước trong đó ô đầu ra (i, j) được đặt thành độ lệch chuẩn của cột đầu vào j giữa các hàng 1 và i. Các NA nên được bỏ qua, trừ khi ô (i, j) của ma trận đầu vào là NA, trong đó ô mẫu (i, j) của ma trận đầu ra cũng phải là NA.

Tôi không thể tìm thấy chức năng tích hợp, vì vậy tôi đã triển khai mã sau. Thật không may, điều này sử dụng một vòng lặp mà kết thúc là hơi chậm cho ma trận lớn. Có chức năng tích hợp nhanh hơn hay ai đó có thể đề xuất phương pháp tiếp cận tốt hơn?

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 
} 

Cảm ơn.

Answer 1

Bạn có thể sử dụng cumsum để tính toán số tiền cần thiết từ thức trực tiếp cho đúng/sd đến hoạt động vectorized trên ma trận:

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 
} 

chỉ để so sánh:

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 

Và khoảng thời gian:

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

Một lần thử khác (Marek's là nhanh hơn)

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