Эффективное вычисление матричного кумулятивного стандартного отклонения в r
Я недавно отправил этот вопрос в списке рассылки r-справки, но не получил ответов, таким образом, я думал, что буду отправлять его здесь также и видеть, были ли какие-либо предложения.
Я пытаюсь вычислить кумулятивное стандартное отклонение матрицы. Я хочу функцию, которая принимает матрицу и возвращает матрицу того же размера, где произведенная ячейка (я, j) установлена на стандартное отклонение входного столбца j между строками 1 и я. NAS должен быть проигнорирован, если ячейка (я, j) самой входной матрицы не является NA, в этом случае ячейка (я, 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
}
Спасибо.
2 ответа
Вы можете использовать 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
Еще одна попытка (Марека быстрее)
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)))
})
}