R中的贴现累积和

Posted 2023-02-16

【中文标题】R中的贴现累积和

【英文标题】：Discounted Cumulative Sum in R

【发布时间】：2020-08-31 06:04:03

【问题描述】：

我正在尝试计算折扣累积总和，其中后面的值更有价值。

假设我有以下数据集：

 dt <- data.table( "year" = c(79,80,81,82,83), "value" = c(5,2,6,8,9))  

> dt
   year value
1:   79     5
2:   80     2
3:   81     6
4:   82     8
5:   83     9

我想要以下输出：

> dt2
year value     DCS    
1:   79     5  5.0000  
2:   80     2  6.5000 
3:   81     6 11.8500
4:   82     8 18.6650 
5:   83     9 25.7985 

贴现累积总和 (DCS) 是通过以 10% 的年贴现率贴现之前的值计算的。因此，对于第一行第二行，DCS 值由 2 + 5*(0.9)^1 给出。对于第三行，DCS 为 6 + (0.9)^1*2 + (0.9)^2*5，以此类推。

形式上，折现和公式由下式给出：

最后，如果可能，最好使用 data.table 解决方案。

【问题讨论】：

我很抱歉，我用 85% 的折扣系数重新计算了它。已编辑。

【参考方案1】：

不是一个正确的答案，而只是其他答案的时间安排。希望这将有助于确定选择哪个选项：

加载库

library(data.table)
library(Rcpp)

创建数据集

set.seed(0L)
dt <- data.table(value = rpois(1e4, 100))

创建必要的函数

app_3 <- function(dt) 
  m <- matrix(0, nrow = nrow(dt), ncol = nrow(dt))
  v <- 0.9**(seq(nrow(dt)) - 1)
  m[lower.tri(m, diag = TRUE)] <- unlist(sapply(rev(seq_along(v)), function(k) head(v, k)))

  dt[, DCS3 := m %*% value]


system.time(
cppFunction("
NumericVector dcs(NumericVector x, double disc) 
    int n = x.size();
    NumericVector res(n);
    res[0] = x[0];
    for (int i=1; i<n; i++) 
        res[i] += x[i] + res[i-1]*disc;
    
    return res;
"))
#   user  system elapsed 
#   0.03    0.16   20.03 

基准测试

res <- bench::mark(time_unit="s",
  app_1 = dt[, DCS1 := sapply(1:.N, function(k) sum(0.9**(k - 1:k)*head(value, k)))],
  app_2 = dt[, DCS2 := dt[, Reduce(function(x, y) 0.9 * x + y, as.list(value), accumulate = TRUE)]],
  app_3 = app_3(dt),

  dt_rcpp = dt[, DCS4 := dcs(value, 0.9)],
  dt_recursive = s <- 0
  dt[, DCS5 := 
    s <- value + s*0.9
    s
  , 1L:nrow(dt)]
  ,

  min_time = 1
)

res

时间安排：

# A tibble: 5 x 13
  expression                   min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result       memory      time   gc        
  <bch:expr>                 <dbl>    <dbl>     <dbl> <bch:byt>    <dbl> <int> <dbl>      <dbl> <list>       <list>      <list> <list>    
1 app_1                   6.34     6.34         0.158    1.12GB    0.315     1     2      6.34  <df[,7] [10~ <df[,3] [5~ <bch:~ <tibble [~
2 app_2                   0.0109   0.0123      71.3    612.34KB   21.8      72    22      1.01  <df[,7] [10~ <df[,3] [2~ <bch:~ <tibble [~
3 app_3                   3.93     3.93         0.255     4.1GB    0.764     1     3      3.93  <df[,7] [10~ <df[,3] [2~ <bch:~ <tibble [~
4 dt_rcpp                 0.000308 0.000337  2681.     195.46KB    6.01   2679     6      0.999 <df[,7] [10~ <df[,3] [2~ <bch:~ <tibble [~
5 dt_recursive            0.00939  0.00972     99.2    294.52KB    6.94    100     7      1.01  <df[,7] [10~ <df[,3] [3~ <bch:~ <tibble [~

另一个 1e6 行的计时：

# A tibble: 3 x 13
  expression                  min  median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result        memory       time   gc        
  <bch:expr>                <dbl>   <dbl>     <dbl> <bch:byt>    <dbl> <int> <dbl>      <dbl> <list>        <list>       <list> <list>    
1 app_2                   1.52    1.52        0.659    53.5MB    6.59      1    10       1.52 <df[,5] [1,0~ <df[,3] [27~ <bch:~ <tibble [~
2 dt_rcpp                 0.00731 0.00942    89.9      15.3MB    0.899   100     1       1.11 <df[,5] [1,0~ <df[,3] [20~ <bch:~ <tibble [~
3 dt_recursive            0.902   0.905       1.10     22.9MB    1.66      2     3       1.81 <df[,5] [1,0~ <df[,3] [4,~ <bch:~ <tibble [~

^{由reprex package (v0.3.0) 于 2020-05-15 创建}

【参考方案2】：

这里有另外 2 个选项。

1) 使用Rcpp，然后在data.table 中通过引用更新：

library(Rcpp)
cppFunction("
NumericVector dcs(NumericVector x, double disc) 
    int n = x.size();
    NumericVector res(n);
    res[0] = x[0];
    for (int i=1; i<n; i++) 
        res[i] += x[i] + res[i-1]*disc;
    
    return res;
")
dt[, DCS := dcs(value, 0.9)]

2) 或者在data.table中递归：

s <- 0
dt[, dcs2 := 
       s <- value + s*0.9
       s
    , 
    1L:nrow(dt)]

#or simply: s <- 0; dt[, dcs2 := s <- value + s*0.9, 1L:nrow(dt)]

输出：

   year value     DCS    dcs2
1:   79     5  5.0000  5.0000
2:   80     2  6.5000  6.5000
3:   81     6 11.8500 11.8500
4:   82     8 18.6650 18.6650
5:   83     9 25.7985 25.7985

---

编辑：回应关于分组的评论：

dt <- data.table(ID=c(1,1,2,2), value=1:4)
dt[, 
    n <- .N
    s <- 0;
    .SD[, 
        s <- value + s*0.9;
        s
      , 
      1L:n]
  ,  
  ID]

输出：

   ID n  V1
1:  1 1 1.0
2:  1 2 2.9
3:  2 1 3.0
4:  2 2 6.7

【讨论】：

你能比较一下函数的时间吗？

第二种方法很棒，我只是从另一个答案中尝试了 sapply，然后我要建议的一个实现以及具有 100 万行数据集的递归 data.table 解决方案，而前两个非常缓慢且消耗内存，data.table 中的递归方式非常神奇。

您能否分享在哪里可以找到有关第二种方法中 data.table 语法的更多信息？由于某种原因，我似乎无法在文档中甚至通过谷歌找到它：/

@ira 这是 Matt Dowle 在***.com/questions/38285789/… 的帖子。从那以后出现了回归，我在项目的 github 站点上发布了一个问题。你愿意分享时间吗？

在具有10k 观察值的数据集上，我在微基准测试的 30 次迭代中具有以下中值时间：sapply: 5945ms、recursively in data.table: 11ms。此外，此答案中提出的第二个选项似乎内存效率更高。不过我还没有尝试过 Rcpp 方法。

【参考方案3】：

也许你可以试试下面的代码。

---

方法 1

通过使用sum直接遵循公式

dt[,DCS:=sapply(1:.N,function(k) sum(0.9**(k-1:k)*head(value,k)))]

---

方法2

使用来自基础 R 的Reduce

dt[,Reduce(function(x,y) 0.9*x+y,as.list(value),accumulate = TRUE)]

---

方法 3

首先，您可以构造一个矩阵m，它给出类似卷积的系数

m <- matrix(0,nrow = nrow(dt),ncol = nrow(dt))
v <- 0.9**(seq(nrow(dt))-1)
m[lower.tri(m,diag = TRUE)] <- unlist(sapply(rev(seq_along(v)),function(k) head(v,k)))

或者使用shift获取矩阵m（感谢@chinsoon12）

x <- 0L:(nrow(dt)-1L); 
m <- t(do.call(cbind, shift(0.9^x, x, fill=0)))

然后就可以运行了

dt[,DCS:=m%*%value]

---

结果

> dt
   year value     DCS
1:   79     5  5.0000
2:   80     2  6.5000
3:   81     6 11.8500
4:   82     8 18.6650
5:   83     9 25.7985