使用滞后计算连续值（矢量化而不是 for 循环）答案

【问题标题】：Using lags to calculate consecutive values (vectorization instead of for-loop)使用滞后计算连续值（矢量化而不是 for 循环）
【发布时间】：2023-03-26 09:06:01
【问题描述】：

如果有人能帮助我在 R 和 data.table 中从 for 循环转换为矢量化方法，我将不胜感激。

我有一个简单的 for 循环，它计算两个变量：库存（之前的库存 + 之前购买的数量 - 消耗量）和消耗量（基于恒定和之前库存的函数，计算一周消耗的数量）

w 周的当前库存取决于 w-1 周之前的库存和 w-1 周之前的消耗量，以及 w-1 周购买的数量。

起始示例如下所示：

> test2
    user       week amount base_consumption inventory consumption seq
 1:    1 2016-07-18  12.00        1.5865385         0           0   1
 2:    1 2016-07-25   0.00        1.5865385         0           0   2
 3:    1 2016-08-01   0.00        1.5865385         0           0   3
 4:    1 2016-08-08   0.00        1.5865385         0           0   4
 5:    1 2016-08-15   0.00        1.5865385         0           0   5
 6:    1 2016-08-22   0.00        1.5865385         0           0   6
 7:    1 2016-08-29  11.25        1.5865385         0           0   7
 8:    1 2016-09-05   0.00        1.5865385         0           0   8
 9:    1 2016-09-12   0.00        1.5865385         0           0   9
10:    1 2016-09-19   0.00        1.5865385         0           0  10
11:    1 2016-09-26   0.00        1.5865385         0           0  11
12:    2 2016-07-18   0.00        0.6923077         0           0   1
13:    2 2016-07-25   0.00        0.6923077         0           0   2
14:    2 2016-08-01   0.00        0.6923077         0           0   3
15:    2 2016-08-08   9.00        0.6923077         0           0   4
16:    2 2016-08-15   0.00        0.6923077         0           0   5
17:    2 2016-08-22   0.00        0.6923077         0           0   6

如果使用以下 for 循环来计算所需的值：

for(i in 1:nrow(test2)){
  if(test2[i,seq] > 1){
    inventory_new <- test2[i-1,inventory+amount-consumption]
    consumption_new <- test2[i-1,inventory_new*(base_consumption/(base_consumption+inventory_new))]
    test2[i,inventory:=inventory_new]
    test2[i,consumption:=consumption_new]
  }
}

这给了我：

> test2
    user       week amount base_consumption inventory consumption seq
 1:    1 2016-07-18  12.00        1.5865385  0.000000   0.0000000   1
 2:    1 2016-07-25   0.00        1.5865385 12.000000   1.4012739   2
 3:    1 2016-08-01   0.00        1.5865385 10.598726   1.3799689   3
 4:    1 2016-08-08   0.00        1.5865385  9.218757   1.3535875   4
 5:    1 2016-08-15   0.00        1.5865385  7.865170   1.3202264   5
 6:    1 2016-08-22   0.00        1.5865385  6.544943   1.2769880   6
 7:    1 2016-08-29  11.25        1.5865385  5.267955   1.2193189   7
 8:    1 2016-09-05   0.00        1.5865385 15.298636   1.4374666   8
 9:    1 2016-09-12   0.00        1.5865385 13.861170   1.4235949   9
10:    1 2016-09-19   0.00        1.5865385 12.437575   1.4070544  10
11:    1 2016-09-26   0.00        1.5865385 11.030521   1.3870384  11
12:    2 2016-07-18   0.00        0.6923077  0.000000   0.0000000   1
13:    2 2016-07-25   0.00        0.6923077  0.000000   0.0000000   2
14:    2 2016-08-01   0.00        0.6923077  0.000000   0.0000000   3
15:    2 2016-08-08   9.00        0.6923077  0.000000   0.0000000   4
16:    2 2016-08-15   0.00        0.6923077  9.000000   0.6428571   5
17:    2 2016-08-22   0.00        0.6923077  8.357143   0.6393443   6

您可以想象，它仅适用于 17 个示例行，但当我尝试将其应用于数千行时，它需要很长时间。我尝试过使用 data.table 的 shift()、rollapply() 和 cumsum()，但没有得到有效的“向量”解决方案。谁能帮我从这里开始？

这是一个输入，如果有人想复制：

> dput(test2)
structure(list(user = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 
2, 2, 2, 2), week = structure(c(17000, 17007, 17014, 17021, 17028, 
17035, 17042, 17049, 17056, 17063, 17070, 17000, 17007, 17014, 
17021, 17028, 17035), class = "Date"), amount = c(12, 0, 0, 0, 
0, 0, 11.25, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0), base_consumption = c(1.58653846153846, 
1.58653846153846, 1.58653846153846, 1.58653846153846, 1.58653846153846, 
1.58653846153846, 1.58653846153846, 1.58653846153846, 1.58653846153846, 
1.58653846153846, 1.58653846153846, 0.692307692307692, 0.692307692307692, 
0.692307692307692, 0.692307692307692, 0.692307692307692, 0.692307692307692
), inventory = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0), consumption = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0), seq = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 1L, 2L, 3L, 4L, 5L, 6L)), row.names = c(NA, -17L), class = c("data.table", 
"data.frame"))

【问题讨论】：

我不认为你可以向量化它，因为inventory 依赖于consumption 的前一个值，而consumption 依赖于inventory 的前一个值。因此，必须先计算每一行，然后才能计算下一行。
是每个用户的base_consumption 常量吗？
是的，每个用户的 base_consumption 都是不变的

标签： r data.table

【解决方案1】：

这是一个使用data.table 的递归选项。您可能希望针对数千行的实际数据对其进行测试。

test2[, { 
        I <- inventory[1L]
        C <- consumption[1L]
        A <- amount[1L]
        b <- base_consumption[1L]
        ans <- .SD[-1L,
            {
                I <- I + A - C
                C <- I * b / (b + I)
                A <- copy(amount)
                b <- copy(base_consumption)
                .(I, C)   
            },
            seq]
        rbindlist(list(data.table(seq=1L, I=I, C=C), ans), 
            use.names=FALSE)
    }, user]

输出：

    user seq         I         C
 1:    1   1  0.000000 0.0000000
 2:    1   2 12.000000 1.4012739
 3:    1   3 10.598726 1.3799690
 4:    1   4  9.218757 1.3535875
 5:    1   5  7.865170 1.3202264
 6:    1   6  6.544943 1.2769880
 7:    1   7  5.267955 1.2193189
 8:    1   8 15.298636 1.4374666
 9:    1   9 13.861170 1.4235950
10:    1  10 12.437575 1.4070545
11:    1  11 11.030520 1.3870384
12:    2   1  0.000000 0.0000000
13:    2   2  0.000000 0.0000000
14:    2   3  0.000000 0.0000000
15:    2   4  0.000000 0.0000000
16:    2   5  9.000000 0.6428571
17:    2   6  8.357143 0.6393443

更快的选择是使用Rcpp，因为将 R 代码转换为 C++ 代码并不难：

library(Rcpp)
cppFunction('
List func(IntegerVector user, NumericVector amount, NumericVector base_consumption) {
    int len = amount.size();
    NumericVector inventory(len);
    NumericVector consumption(len);

    for (int i = 1; i < len; i++) {
        if (user[i-1]==user[i]) {
            inventory[i] = inventory[i-1] + amount[i-1] - consumption[i-1];
            consumption[i] = inventory[i] * (base_consumption[i-1]  / (base_consumption[i-1] + inventory[i]));
        }
    }

    return List::create(Named("inventory")=inventory,
        Named("consumption")=consumption);
}
')
test2[, c("inventory", "consumption") := func(user, amount, base_consumption)]

输出：

    user       week amount base_consumption inventory consumption seq
 1:    1 2016-07-18  12.00        1.5865385  0.000000   0.0000000   1
 2:    1 2016-07-25   0.00        1.5865385 12.000000   1.4012739   2
 3:    1 2016-08-01   0.00        1.5865385 10.598726   1.3799690   3
 4:    1 2016-08-08   0.00        1.5865385  9.218757   1.3535875   4
 5:    1 2016-08-15   0.00        1.5865385  7.865170   1.3202264   5
 6:    1 2016-08-22   0.00        1.5865385  6.544943   1.2769880   6
 7:    1 2016-08-29  11.25        1.5865385  5.267955   1.2193189   7
 8:    1 2016-09-05   0.00        1.5865385 15.298636   1.4374666   8
 9:    1 2016-09-12   0.00        1.5865385 13.861170   1.4235950   9
10:    1 2016-09-19   0.00        1.5865385 12.437575   1.4070545  10
11:    1 2016-09-26   0.00        1.5865385 11.030520   1.3870384  11
12:    2 2016-07-18   0.00        0.6923077  0.000000   0.0000000   1
13:    2 2016-07-25   0.00        0.6923077  0.000000   0.0000000   2
14:    2 2016-08-01   0.00        0.6923077  0.000000   0.0000000   3
15:    2 2016-08-08   9.00        0.6923077  0.000000   0.0000000   4
16:    2 2016-08-15   0.00        0.6923077  9.000000   0.6428571   5
17:    2 2016-08-22   0.00        0.6923077  8.357143   0.6393443   6

【讨论】：

非常感谢！两种解决方案都有效并显着提高了速度。除了性能之外，我还应该使用 C++ 解决方案而不是 data.table 解决方案吗？
好的，这是我的结论（以 1200 万行为基准）：for-loop（如上所述）：~ 18 小时，data.table 中的函数内的 for-loop（不包括在这篇文章中） : ~ 6 小时, chinsoon12 (recursive data.table) 建议的解决方案: