Как записать накопительный расчет в data.table

Question

Последовательный, совокупный расчет

Мне нужно сделать расчет временных рядов, где значение, рассчитанное в каждой строке, зависит от результата, вычисленного в предыдущей строке. Я надеюсь использовать удобство data.table, Актуальной проблемой является гидрологическая модель - вычисление совокупного водного баланса, добавление осадков на каждом временном шаге и вычитание стока и испарения в зависимости от текущего объема воды. Набор данных включает в себя различные бассейны и сценарии (группы). Здесь я буду использовать более простую иллюстрацию проблемы.

Упрощенный пример расчета выглядит так, для каждого временного шага (строки) i:

 v[i] <- a[i] + b[i] * v[i-1]

a а также b являются векторами значений параметров, и v является вектором результата. Для первого ряда (i == 1) начальная стоимость v принимается как v0 = 0,

Первая попытка

Моей первой мыслью было использовать shift() в data.table, Минимальный пример, включая желаемый результат v.ans, является

library(data.table)        # version 1.9.7
DT <- data.table(a = 1:4, 
                 b = 0.1,
                 v.ans = c(1, 2.1, 3.21, 4.321) )
DT
#    a   b v.ans
# 1: 1 0.1 1.000
# 2: 2 0.1 2.100
# 3: 3 0.1 3.210
# 4: 4 0.1 4.321

DT[, v := NA]   # initialize v
DT[, v := a + b * ifelse(is.na(shift(v)), 0, shift(v))][]
#    a   b v.ans v
# 1: 1 0.1 1.000 1
# 2: 2 0.1 2.100 2
# 3: 3 0.1 3.210 3
# 4: 4 0.1 4.321 4

Это не работает, потому что shift(v) дает копию оригинального столбца v, сдвинут на 1 ряд. Это не зависит от назначения v,

Я также рассмотрел построение уравнения с использованием cumsum() и cumprod(), но это тоже не сработает.

Метод грубой силы

Поэтому для удобства я прибегаю к циклу for внутри функции:

vcalc <- function(a, b, v0 = 0) {
  v <- rep(NA, length(a))      # initialize v
  for (i in 1:length(a)) {
    v[i] <- a[i] + b[i] * ifelse(i==1, v0, v[i-1])
  }
  return(v)
}

Эта накопительная функция отлично работает с data.table:

DT[, v := vcalc(a, b, 0)][]
#    a   b v.ans     v
# 1: 1 0.1 1.000 1.000
# 2: 2 0.1 2.100 2.100
# 3: 3 0.1 3.210 3.210
# 4: 4 0.1 4.321 4.321
identical(DT$v, DT$v.ans)
# [1] TRUE

Мой вопрос

Мой вопрос, могу ли я написать этот расчет более кратким и эффективным data.table Кстати, без необходимости использовать цикл for и / или определение функции? С помощью set() возможно?

Или есть лучший подход все вместе?

Изменить: лучший цикл

Решение Дэвида Rcpp ниже вдохновило меня удалить ifelse() от for цикл:

vcalc2 <- function(a, b, v0 = 0) {
  v <- rep(NA, length(a))
  for (i in 1:length(a)) {
    v0 <- v[i] <- a[i] + b[i] * v0
  }
  return(v)
}

vcalc2() на 60% быстрее, чем vcalc(),

Answer 1

Это может быть не на 100% то, что вы ищете, так как он не использует "data.table-way" и все еще использует цикл for. Однако этот подход должен быть быстрее (я предполагаю, что вы хотите использовать data.table и data.table-way для ускорения вашего кода). Я использую Rcpp, чтобы написать короткую функцию под названием HydroFun, который может быть использован в R, как и любая другая функция (сначала нужно просто вызвать функцию). Мои интуитивные ощущения говорят мне, что способ data.table (если он существует) довольно сложен, потому что вы не можете вычислить решение в закрытой форме (но я могу ошибаться в этом вопросе...).

Мой подход выглядит так:

Функция Rcpp выглядит следующим образом (в файле: hydrofun.cpp):

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector HydroFun(NumericVector a, NumericVector b, double v0 = 0.0) {
  // get the size of the vectors
  int vecSize = a.length();

  // initialize a numeric vector "v" (for the result)
  NumericVector v(vecSize);

   // compute v_0
  v[0] = a[0] + b[0] * v0;

  // loop through the vector and compute the new value
  for (int i = 1; i < vecSize; ++i) {
    v[i] = a[i] + b[i] * v[i - 1];
  }
  return v;
}

Для получения и использования функции в R вы можете сделать:

Rcpp::sourceCpp("hydrofun.cpp")

library(data.table)
DT <- data.table(a = 1:4, 
                 b = 0.1,
                 v.ans = c(1, 2.1, 3.21, 4.321))

DT[, v_ans2 := HydroFun(a, b, 0)]
DT
# a   b v.ans v_ans2
# 1: 1 0.1 1.000  1.000
# 2: 2 0.1 2.100  2.100
# 3: 3 0.1 3.210  3.210
# 4: 4 0.1 4.321  4.321

Что дает результат, который вы ищете (по крайней мере, с точки зрения ценности).

Сравнение скоростей показывает ускорение примерно в 65 раз.

library(microbenchmark)
n <- 10000
dt <- data.table(a = 1:n,
                 b = rnorm(n))

microbenchmark(dt[, v1 := vcalc(a, b, 0)],
               dt[, v2 := HydroFun(a, b, 0)])
# Unit: microseconds
# expr                                min        lq       mean    median         uq       max neval
# dt[, `:=`(v1, vcalc(a, b, 0))]    28369.672 30203.398 31883.9872 31651.566 32646.8780 68727.433   100
# dt[, `:=`(v2, HydroFun(a, b, 0))]   381.307   421.697   512.2957   512.717   560.8585  1496.297   100

identical(dt$v1, dt$v2)
# [1] TRUE

Это тебе как-нибудь помогает?

Answer 2

Я думаю Reduce вместе с accumulate = TRUE является широко используемой техникой для этих типов вычислений (см., например, рекурсивное использование выходных данных в качестве входных данных для функции). Это не обязательно быстрее, чем хорошо написанный цикл *, и я не знаю, как data.table- Если вы верите, то все же я хочу предложить это для вашего набора инструментов.

DT[ , v := 0][
  , v := Reduce(f = function(v, i) a[i] + b[i] * v, x = .I[-1], init = a[1], accumulate = TRUE)]

DT
#    a   b v.ans     v
# 1: 1 0.1 1.000 1.000
# 2: 2 0.1 2.100 2.100
# 3: 3 0.1 3.210 3.210
# 4: 4 0.1 4.321 4.321

Объяснение:

Установите начальное значение v в 0 (v := 0). использование Reduce применить функцию f на целом векторе номеров строк, кроме первой строки (x = .I[-1]). Вместо этого добавьте a[1] к началу x (init = a[1]).Reduce затем "последовательно применяет f к элементам [...] слева направо". Последовательные комбинации уменьшения "накапливаются" (accumulate = TRUE).

---

* Смотрите, например, здесь, где вы также можете прочитать больше о Reduce в этом разделе.