如何在 data.table 中编写累积计算
How to write a cumulative calculation in data.table
顺序累积计算
我需要做一个时间序列计算，每行计算的值取决于上一行计算的结果。我希望使用 data.table 的便利。实际问题是一个水文模型——累积水量平衡计算，在每个时间步增加降雨量，并减去径流和蒸发量作为当前水量的函数。该数据集包括不同的流域和场景（组）。在这里我将使用更简单的问题来说明。
计算的简化示例如下所示，对于每个时间步长（行）i：
v[i] <- a[i] + b[i] * v[i-1]
a和b是参数值的向量，v是结果向量。对于第一行 (i == 1)，v 的初始值取为 v0 = 0.
第一次尝试
我的第一个想法是在 data.table 中使用 shift()。一个最小的例子，包括期望的结果 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 循环 and/or 函数定义？也许使用 set()？
或者是否有更好的方法？
编辑：更好的循环
下面 David 的 Rcpp 解决方案启发了我从 for 循环中删除 ifelse()：
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() 比 vcalc().
快 60%
它可能不是您要查找的 100%，因为它不使用 "data.table-way" 并且仍然使用 for 循环。但是，这种方法应该更快（我假设您想使用 data.table 和 data.table 方式来加速您的代码）。我利用 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
这对你有什么帮助吗？
我认为 Reduce 和 accumulate = TRUE 是此类计算的常用技术（参见 recursively using the output as an input for a function）。它不一定比编写良好的循环*更快，而且我不知道你认为它是如何 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 应用于行号 除了 第一行 (x = .I[-1]) 的整数向量。而是将 a[1] 添加到 x (init = a[1]) 的开头。
Reduce 然后 "successively applies f to the elements [...] from left to right"。
连续的 reduce 组合是 "accumulated" (accumulate = TRUE).
*参见例如here, where you also can read more about Reduce in this section。
顺序累积计算
我需要做一个时间序列计算，每行计算的值取决于上一行计算的结果。我希望使用 data.table 的便利。实际问题是一个水文模型——累积水量平衡计算，在每个时间步增加降雨量，并减去径流和蒸发量作为当前水量的函数。该数据集包括不同的流域和场景（组）。在这里我将使用更简单的问题来说明。
计算的简化示例如下所示，对于每个时间步长（行）i：
v[i] <- a[i] + b[i] * v[i-1]
a和b是参数值的向量，v是结果向量。对于第一行 (i == 1)，v 的初始值取为 v0 = 0.
第一次尝试
我的第一个想法是在 data.table 中使用 shift()。一个最小的例子，包括期望的结果 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 循环 and/or 函数定义？也许使用 set()？
或者是否有更好的方法？
编辑：更好的循环
下面 David 的 Rcpp 解决方案启发了我从 for 循环中删除 ifelse()：
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() 比 vcalc().
它可能不是您要查找的 100%，因为它不使用 "data.table-way" 并且仍然使用 for 循环。但是，这种方法应该更快（我假设您想使用 data.table 和 data.table 方式来加速您的代码）。我利用 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
这对你有什么帮助吗？
我认为 Reduce 和 accumulate = TRUE 是此类计算的常用技术（参见 recursively using the output as an input for a function）。它不一定比编写良好的循环*更快，而且我不知道你认为它是如何 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 应用于行号 除了 第一行 (x = .I[-1]) 的整数向量。而是将 a[1] 添加到 x (init = a[1]) 的开头。
Reduce 然后 "successively applies f to the elements [...] from left to right"。
连续的 reduce 组合是 "accumulated" (accumulate = TRUE).
*参见例如here, where you also can read more about Reduce in this section。