指数加重移動平均の高速R実装ですか？

Question

私は指数関数的に加重移動平均を実行したいです（パラメータ化はhereと定義されています）。以下の最初の試みよりも優れた実装がありますか？

私の最初の試みだった：

ewma <- function(x, a) { 
    n <- length(x) 
    s <- rep(NA,n) 
    s[1] <- x[1] 
    if (n > 1) { 
    for (i in 2:n) { 
     s[i] <- a * x[i] + (1 - a) * s[i-1] 
    } 
    } 
    return(s) 
} 

y <- 1:1e7 
system.time(s <- ewma(y,0.5)) 
#user system elapsed 
# 2.48 0.00 2.50 

私の第二の試みで、私はベクトル化でもっとうまくできると思っていた：

ewma_vectorized <- function(x,a) { 
    a <- 0.1 
    n <- length(x) 
    w <- cumprod(c(1, rep(1-a, n-1))) 
    x1_contribution <- w * x[1] 
    w <- a * w 
    x <- x[-1] 
    s <- apply(as.array(1:(n-1)), 1, function(i,x,w){sum(w[i:1] * x[1:i])}, x=x, w=w) 
    s <- x1_contribution + c(0,s) 
    return(s) 
} 

system.time(s <- ewma_vectorized(y,0.5)) 
# I stopped the program after it continued to run for 4min 

私は結果によってあまりにも驚いていたべきではないと思います私の2回目の試みで。ベクトル化の試みはかなり醜い試みでした。しかし、のようなものがのようなものでなければならない。

UPDATE：

私はより良い実装hereを見つけたし、次のようにそれを適応：

ewma_vectorized_v2 <- function(x, a) { 
    s1 <- x[1] 
    sk <- s1 
    s <- vapply(x[-1], function(x) sk <<- (1 - a) * x + a * sk, 0) 
    s <- c(s1, s) 
    return(s) 
} 

system.time(s <- ewma_vectorized_v2(y,0.5)) 
# user system elapsed 
# 1.74 0.01 1.76 

Answer 1

あなたはstats::filterでこれを行うことができます：あなたは本当に速いようにしたい場合は、私はあなたがRcppを使用して、次の例に見ることができるようにあなたは、C/C++を使用すべきだと思う。とにかく

、

ewma.filter <- function (x, ratio) { 
    c(filter(x * ratio, 1 - ratio, "recursive", init = x[1])) 
} 
set.seed(21) 
x <- rnorm(1e4) 
all.equal(ewma.filter(x, 0.9), ewma(x, 0.9)) 
# [1] TRUE 

これはあなたの最初の試みのコンパイル済みのバージョンよりも少し速いです：

ewma <- compiler::cmpfun(function(x, a) { 
    n <- length(x) 
    s <- rep(NA,n) 
    s[1] <- x[1] 
    if (n > 1) { 
    for (i in 2:n) { 
     s[i] <- a * x[i] + (1 - a) * s[i-1] 
    } 
    } 
    return(s) 
}) 
microbenchmark(ewma.filter(x,0.9), ewma(x, 0.9)) 
Unit: microseconds 
       expr  min  lq median  uq  max neval 
ewma.filter(x, 0.9) 318.508 341.7395 368.737 473.254 1477.000 100 
     ewma(x, 0.9) 1364.997 1403.4015 1458.961 1503.876 2221.252 100 

Answer 2

私のマシン（R 3.3.2 Windowsの場合）で、あなた最初のループ〜16秒かかります。 jitコンパイルを有効にするには、関数定義の前に行compiler::enableJIT(2)を追加します。コードは〜1秒で実行されます。

library(Rcpp) 

sourceCpp(
    code = 
    " 
    #include <Rcpp.h> 
    // [[Rcpp::export]] 
    Rcpp::NumericVector ewmaRcpp(Rcpp::NumericVector x, double a){ 
     int n = x.length(); 
     Rcpp::NumericVector s(n); 
     s[0] = x[0]; 
     if (n > 1) { 
     for (int i = 1; i < n; i++) { 
      s[i] = a * x[i] + (1 - a) * s[i-1]; 
     } 
     } 
     return s; 
    } 

    ") 

y <- 1:1e7 
system.time(s2 <- ewmaRcpp(y,0.5)) 
# user system elapsed 
# 0.06 0.01 0.07 

Answer 3

@digEmAllはRcppバージョンでは非常に親切でしたが、TTRパッケージや著者のメモとして、10年前に（現在は存在していない）Rグラフギャラリーの記事で使用したstats::filter()のアプローチを使用することができます。

とにかく、素早いシュートアウト戦は、おそらく我々はパラメータ化間違って持って意味している...同じくらい速くRcppバージョンを示しています

R> sourceCpp("/tmp/ema.cpp") 

R> library(TTR) 

R> library(microbenchmark) 

R> y <- as.numeric(1:1e6) # else the sequence creates ints 

R> microbenchmark(ewmaRcpp(y,0.5), EMA(y, n=10)) 
Unit: milliseconds 
      expr  min  lq  mean median  uq  max neval cld 
ewmaRcpp(y, 0.5) 2.43666 2.45705 3.06699 2.47283 2.51439 9.76883 100 a 
    EMA(y, n = 10) 15.13208 15.37910 21.36930 15.59278 17.26318 76.45934 100 b 
R> 

実は、lambda=0.5が半分に相当するであろう非常に強い減衰でありますある日の生活、またはN=1。私がそれを使用すると、ギャップ はさらに広いです。

完全にするために、ファイル全体がどのちょうどRcpp::sourceCpp() -edことができます。

#include <Rcpp.h> 
// [[Rcpp::export]] 
Rcpp::NumericVector ewmaRcpp(Rcpp::NumericVector x, double a){ 
    int n = x.length(); 
    Rcpp::NumericVector s(n); 
    s[0] = x[0]; 
    if (n > 1) { 
    for (int i = 1; i < n; i++) { 
     s[i] = a * x[i] + (1 - a) * s[i-1]; 
    } 
    } 
    return s; 
} 

/*** R 
library(TTR) 
library(microbenchmark) 

y <- as.numeric(1:1e6) # else the sequence creates ints 
microbenchmark(ewmaRcpp(y,0.5), EMA(y, n=1)) 
*/