Удивительно медленное стандартное отклонение в R

Question

Я рассчитываю стандартные отклонения в расширяющемся окне, где в каждой точке я пересчитываю стандартное отклонение.Это кажется довольно простой вещью, которая должна быть относительно быстрой.Однако это занимает намного больше времени, чем вы думаете (~ 45 секунд).Я что-то здесь упускаю?В Matlab это довольно быстро.

t0 <- proc.time()[[3]]
z <- rep(0, 7000)
x <- rnorm(8000)
for(i in 1000:8000){
##    print(i)
    z[i] <- sd(x[1:i])
}
print(proc.time()[[3]]- t0)

Answer 1

Вы также можете попробовать алгоритм, который обновляет стандартное отклонение (ну, на самом деле, сумму квадратов различий от среднего), когда вы идете. В моей системе это сокращает время с ~ 0,8 до ~ 0,002 с.

n <- length(x)
m <- cumsum(x)/(1:n)
m1 <- c(NA,m[1:(n-1)])
ssd <- (x-m)*(x-m1)
v <- c(0,cumsum(ssd[-1])/(1:(n-1)))
z <- sqrt(v)

Подробнее см. http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance.

Также см. Ответы на этот вопрос: Эффективный расчет кумулятивного стандартного отклонения матрицы в r

Answer 2

Отредактировано, чтобы исправить некоторые опечатки, извините.

На моем компьютере это займет ~ 1,3 секунды:

t0 <- proc.time()[[3]]
x <- rnorm(8000)
z <- sapply(1000:8000,function(y){sd(x[seq_len(y)])})
print(proc.time()[[3]]- t0)

и я бы поспорил, что есть еще более быстрые способы сделать это. Избегайте явных for петель!

Answer 3

Когда несколько дней назад несколько схожий вопрос о кумулятивной дисперсии и кумулятивной операции по куртозу возник в результате некоторых проблем, вот что я предложил:

daily <- rnorm(1000000)
mbar <- mean(daily)
cumvar <-  cumsum( (daily-cumsum(daily)/1:length(daily) )^2)
cumskew <- cumsum( (daily-cumsum(daily)/1:length(daily))^3)/cumvar^(3/2)

Это, конечно, быстрее, чем метод Саппи, но можетбыть сравнимым с Аароном.

 system.time( cumvar <-  cumsum( (daily-cumsum(daily)/1:length(daily) )^2) )
   user  system elapsed 
  0.037   0.026   0.061 
 system.time(cumsd <- sqrt(cumvar) )
   user  system elapsed 
  0.009   0.005   0.013