r - 对具有后向依赖性的 R 循环进行向量化

Question

我有一个随机向量vec并且想在不使用循环的情况下创建一个新向量L。 L 的新元素依赖于 L 和 vec 的旧元素。

set.seed(0)

vec <- rnorm(20,0)
i = 2; 
N <- length(vec) -1
L <- numeric(N-1)
constant <- 0.6

while (i < N){  

  L[i] = vec[i + 1] - vec[i] - constant * L[i - 1]

  i <- i + 1          
}

L
#  [1]  0.0000000  1.6560326 -1.0509895 -0.2271942 -1.8182750  1.7023480 -0.3875622  0.5214906  2.0975262 -2.8995756  0.1771427
# [12] -0.4549334  1.1311555 -0.6884468  0.3007724  0.4832709 -1.4341071  2.1880687

Answer 1

你想要

L[1] = 0
L[i] = -constant * L[i - 1] + (vec[i + 1] - vec[i]),  i = 2, 3, ..., 

让dv <- diff(vec) ，第二行变为

L[i] = -constant * L[i - 1] + dv[i],  i = 2, 3, ...

具有滞后 1 自相关性的 AR1 过程 -constant和创新dv[-1] 。 AR1进程可以通过filter高效生成使用“递归”方法。

dv <- diff(vec)
L <- c(0, filter(dv[-1], -constant, "recursive"))

# [1]  0.0000000  1.6560326 -1.0509895 -0.2271942 -1.8182750  1.7023480
# [7] -0.3875622  0.5214906  2.0975262 -2.8995756  0.1771427 -0.4549334
#[13]  1.1311555 -0.6884468  0.3007724  0.4832709 -1.4341071  2.1880687
#[19] -2.9860629

我猜你的意思是while (i <= N)在你的问题中。如果您确实想要i < N ，那么你必须去掉上面的最后一个元素。可以通过以下方式完成

dv <- diff(vec)
L <- c(0, filter(dv[2:(length(dv) - 1)], -constant, "recursive"))

---

几个小时后...

我被 Rui Barradas 引起注意的基准。简称vec ，任何方法都足够快。长久以来vec , filter肯定更快，但实际上会受到强制，如 filter期望并返回一个“ts”(时间序列)对象。最好直接调用其主力 C 例程:

AR1_FILTER <- function (x, filter, full = TRUE) {
  n <- length(x)
  AR1 <- .Call(stats:::C_rfilter, as.double(x), as.double(filter), double(n + 1L))
  if (!full) AR1 <- AR1[-1L]
  AR1
  }

dv <- diff(vec)
L <- AR1_FILTER(dv[-1], -constant)
#L <- AR1_FILTER(dv[2:(length(dv) - 1)], -constant)

我对比较 AR1_FILTER 不感兴趣与R级循环。我将它与 filter 进行比较.

library(microbenchmark)
v <- runif(100000)
microbenchmark("R" = c(0, filter(v, -0.6, "recursive")),
               "C" = AR1_FILTER(v, -0.6))

Unit: milliseconds
 expr      min       lq     mean   median       uq       max neval
    R 6.803945 7.987209 11.08361 8.074241 9.131967 54.672610   100
    C 2.586143 2.606998  2.76218 2.644068 2.660831  3.845041   100