R高效递归

Question

让我们假设我们观察到的数据在一个为期一年的基础上，每月(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)  

(一般来说，参数β和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-function了。

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

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

另一种选择是使用简单的for-loop。但这并不像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产生了更多的开销，但是一旦创建它，这种开销就不会很大。因此，这取决于您要使用的实现。