R中向量的指数移动平均值

Exponential moving average of a vector in R

我有一个简单的向量如下:

x = c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

我正在尝试使用以下函数找到此向量的滚动 EMA -

library(TTR)
y = EMA(x, 5)

我得到的结果如下 -

 [1]     NA     NA     NA     NA 13.33400 13.22267 13.52844 14.44563 16.51042 16.88695

但是,我想要如下结果 -

 [1]     14.24 14.03 13.06 13.43 13.33400 13.22267 13.52844 14.44563 16.51042 16.88695
  1. 第一个值应与原始向量中的值相同
  2. 第二个值应该是第一个和第二个值的 EMA
  3. 第三个值应该是前三个值的EMA 矢量
  4. 第四个值应该是向量中前四个值的 EMA

其余计算由函数正确处理 EMA

我试过的解决方案-

  1. 运行 下面的命令—— zoo::rollapplyr(x, width = 5, FUN = EMA, partial = TRUE) 将给出错误,因为 EMA 有自己的滚动 window。

  2. 使用函数 stats::filter 有效,但答案不正确,因为我不确定比率参数的正确值。 这是一个自定义函数。

ema_2 <- function (k, width) {
  ratio <- 2/(width + 1)
  c(stats::filter(k * ratio, 1 - ratio, "convolution", init = k[1]))
}

理想的解决方案最多应该是 TTR 库的 EMA 函数所用计算时间的两倍。

以下是 Waldi 和 Andre 分享的 2 个解决方案的性能结果。

              expr     min       lq     mean   median       uq      max neval cld
    TTR::EMA(x, 5) 433.593 457.5815 500.9478 477.0535 530.7105  1128.49  1000   a
        EMA3(x, 5) 445.388 468.9585 515.2009 490.4345 546.5025  1843.46  1000   a
 rollmeanEMA(x, 5) 436.689 461.0885 535.7035 481.8815 538.3150 33122.75  1000   a

谢谢!

这给出了期望的结果:

require(TTR)

x <- c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

rollmeanEMA <- function(vec, len) {
  c(cumsum(vec[1:(len-1)]) / seq_along(vec[1:(len-1)]),
    EMA(vec, len)[len:length(vec)])
}

rollmeanEMA(x,5)
#[1] 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#[9] 16.51042 16.88695

编辑:正如我在评论中所建议的那样,将 NA 部分替换为 mean()。这提供了巨大的加速。另外,去掉了周围的条件。

y <- rnorm(1000000)

system.time( rollmeanEMA(y,10000) )
#   user  system elapsed
#  0.031   0.003   0.034

system.time( EMA(y,10000) )
#   user  system elapsed
#  0.018   0.002   0.019

添加了NA“处理”:

rollmeanEMA <- function(vec, len) {
  v_n <- !is.na(vec)
  c( vec[is.na(vec)],
     cumsum(vec[v_n][1:(len-1)]) / seq_along(vec[v_n][1:(len-1)]),
     EMA(vec[v_n], len)[len:length(vec[v_n])])
}

查看C source code of EMA表明第一个值是平均值window:

    /* Raw mean to start EMA */
    double seed = 0.0;
    for(i = first; i < first + i_n; i++) {
      d_result[i] = NA_REAL;
      seed += d_x[i] / i_n;
    }
    d_result[first + i_n - 1] = seed;

这可以很容易地计算出来以替换 NA:

x = c(14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)

EMA2 <- function(x,n) {
  y = TTR::EMA(x, n)
  noNA <- which.min(is.na(x))
  y[noNA:(noNA+n-2)] <- cumsum(x[noNA:(noNA+n-2)])/1:(n-1)
  y
}

EMA2(x,5)
#>  [1] 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#>  [9] 16.51042 16.88695

这也适用于领先的 NA:

x = c(NA, NA, 14.24, 13.82, 12.75, 12.92, 12.94, 13.00, 14.14, 16.28, 20.64, 17.64)
EMA2(x,5)
#> [1]       NA       NA 14.24000 14.03000 13.60333 13.43250 13.33400 13.22267 13.52844 14.44563
#> [11] 16.51042 16.88695

这个短向量的计算开销很小,在更长的向量上应该更好:

microbenchmark::microbenchmark(TTR::EMA(x,5),EMA2(x,5),times=1000)

#> Unit: microseconds
#>           expr   min    lq     mean median     uq   max neval cld
#> TTR::EMA(x, 5) 157.7 161.8 181.6156  164.0 180.55 593.5  1000  a 
#>     EMA2(x, 5) 164.2 167.5 193.0643  170.6 193.20 857.1  1000   b