R高效递归

Question

让我们假设我们每月观察 1 年的数据 (k = 12)。让我们还假设以下递归结构（下面的示例）。

我想知道，考虑到以下几个方面（效率是绝对必要的），如何最有效地对该结构进行编程：

1年的时间是为了简单起见，一般不固定为这个值。因此，每次都要重新定义这个值，
必须存储每月的计算结果（我想到了一个列表？），因为我需要具体的值来进行进一步的计算。

MWE：

beta = c(0.95, 0.89, 0.23, 0.96, 0.98, 0.79, 0.99, 0.85, 0.97) 
dim(beta) <- c(3, 3)
phi = c(0.45, 0.12, 0.98, 0.66, 0.79, 0.24, 0.33, 0.19, 0.27)
dim(phi) <- c(3, 3)

（一般来说，参数beta和phi是无法观察到的，需要估计，所以这只是一个简单的假设。）

在k-1中我们计算出如下形式的矩阵：

K1 <- beta %*% t(beta)

在 k-2 中我们有类似的表达式，但现在它也取决于以前的值乘以 phi:

K2 <- beta %*% t(beta) + phi %*% K1

k-3期与k-2期相同，只是K1变为K2。所以从这里开始它是递归的并重复自己。所以对于最后一次观察 (k-11) 它与 K10 的公式相同。

Answer 1

首先你总是需要术语beta %*% t(beta)，所以存储它而不是在每次迭代中计算它：

btb <- beta %*% t(beta)

现在您可以使用 Reduce 功能了。

Reduce(f = function(kPrevious,someOther) {btb + phi %*% kPrevious},
       x = 1:10,
       init = btb,
       accumulate = TRUE
)

参数accumulate=TRUE保存所有中间结果

另一种方法是使用简单的 for 循环。但这不如 Reduce 高效。（尽管它“相同”）：

  history <- list(btb)
  for(k in 1:10) {
    history <- c(history,list(btb + phi %*% history[[length(history)]]))
  }
  history

为了计算一些基准，让我们根据迭代次数将所有内容放入一个函数中：

ff0 <- function(n=10) {
  Reduce(f = function(kPrevious,someOther) {btb + phi %*% kPrevious},
         x = 1:n,
         init = btb,
         accumulate = TRUE
  )
}

ff <- function(n=10) {
  history <- list(btb)
  for(k in 1:n) {
    history <- c(history,list(btb + phi %*% history[[length(history)]]))
  }
  history
}

现在，让我们看看这两个函数是否给出相同的输出：

all(mapply(FUN = function(x,y) {all(x==y)},ff0(),ff()))

要获得基准，让我们使用 microbenchmark:

microbenchmark::microbenchmark(ff0(),ff())
# Unit: microseconds
#  expr  min    lq   mean median    uq  max neval
# ff0() 20.6 22.80 25.388   24.1 25.20 79.2   100
#  ff() 11.4 12.25 13.922   13.0 13.65 80.0   100

microbenchmark::microbenchmark(ff0(100),ff(100))
#Unit: microseconds
#     expr   min     lq    mean median     uq   max neval
# ff0(100) 163.1 172.50 191.193 186.35 198.00 320.6   100
#  ff(100) 143.0 151.45 169.125 164.75 174.05 296.1   100

microbenchmark::microbenchmark(ff0(1000),ff(1000))
#Unit: milliseconds
#      expr    min      lq     mean  median      uq     max neval
# ff0(1000) 1.5680 1.62505 1.698665 1.67135 1.74185  2.6058   100
#  ff(1000) 4.6626 4.77065 5.293329 4.87640 5.10260 11.2255   100

microbenchmark::microbenchmark(ff0(10000),ff(10000))
#Unit: milliseconds
#       expr      min       lq      mean   median        uq      max neval
# ff0(10000)  15.9035  16.4428  17.28437  17.0618  17.71085  22.1506   100
#  ff(10000) 468.8924 484.7448 494.79808 489.4221 496.60075 581.0389   100

正如我们所见，for-loop 实现似乎更适合迭代次数相对较少的计算。但是随着迭代次数的增加，Reduce 变得更快。我认为原因是 Reduce 比简单的 foor-loop 产生更多的开销，但是这个开销在创建后并不重要。因此，您将使用哪种实现取决于您。

R高效递归

R highly efficient recursion

recursion

r

memory-efficient