Блог

Предположим, что мы ежемесячно наблюдаем данные за горизонтом в 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.

Прежде всего, вам всегда нужен термин 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 цикл, но эти накладные расходы не имеют большого веса после его создания.
Так что вам решать, какую реализацию вы будете использовать.