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

Question

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

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 有自己的滚动窗口。

2. 使用函数stats::filter 有效，但答案不正确，因为我不确定 ratio 参数的正确值。 Fast R implementation of an Exponentially Weighted Moving Average? 这是一个自定义函数。

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

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

以下是 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

谢谢！

Answer 1

这给出了预期的结果：

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

编辑：正如我在 cmets 中建议的那样，将 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])])
}

Answer 2

查看C source code of EMA 显示第一个值是平均窗口的平均值：

    /* 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