r - 使用滞后计算连续值(矢量化而不是 for 循环)

Question

如果有人能帮助我在 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，如果有人想复制:

> 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"))

Answer 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