18

顺序累积计算

我需要进行时间序列计算，其中每一行计算的值取决于前一行计算的结果。我希望使用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)给出了原始 column 的副本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()。

4

2 回答 2

7

它可能不是您正在寻找的 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

这对你有任何帮助吗？

于 2016-11-03T23:50:28.633 回答
2

我认为Reducewithaccumulate = TRUE是这些类型计算的常用技术（参见例如递归地使用输出作为函数的输入）。它不一定比编写良好的循环更快*，而且我不知道data.table您认为它是多么-esque，但我仍想为您的工具箱推荐它。

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 应用于元素 [...]”。连续的 reduce 组合是“累积的”（）。fx = .I[-1]a[1]xinit = a[1]Reduceaccumulate = TRUE

*参见例如此处，您还可以Reduce在此部分阅读更多信息。

于 2016-11-04T11:06:22.767 回答