如何在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()或许使用？

或者是否有更好的方法？

编辑:更好的循环

下面David的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

这对你有什么帮助吗？