R中滞后的向量化差异函数答案

【问题标题】：Vectorizing diff function on lags in RR中滞后的向量化差异函数
【发布时间】：2014-06-11 22:04:23
【问题描述】：

我想做一个函数来做以下事情：

c <- rnorm(100)
n <- 10
sum.diff<- integer(n)

for (k in 1:n) {
   sum.diff[k] <- sum(diff(c, lag=k))
}

通过矢量化而不是循环。具体来说，我想发送一个向量（这里是'c'）和一个滞后值向量（这里是'1：n'），并得到第k个条目中第k个差异的总和输出向量（此处为“sum.lags”）。

例如，c <- 1:100 和 n <- 10 应该产生：

99 196 291 ... 900

对应于：

sum(diff(1:100,lag=1)) sum(diff(1:100,lag=2)) sum(diff(1:100,lag=3)) ... sum(diff(1:100,lag=10)) 有任何想法吗？

【问题讨论】：

我可以想办法用sapply 来代替for 循环，但它并没有真正矢量化，也不是更快（事实上，有点慢）。将这个向量化的问题是每次迭代都要处理不同数量的数字。您是否希望矢量化，因为您的生产代码正在处理相当大的数据集？如果是这样，我会推荐 compiler::cmpfun，或者更好的是 Rcpp。
这是一个更大的数据集（因此非常慢），但它是针对一个特定的问题，可能不值得投入那么多时间。证明 rpp 的合理性有点困难（我有从未使用过它，但一直渴望学习它，虽然我确实有 C++ 经验），但我会检查两者。如果不出意外，它可以作为一种熟悉我应该知道的事情的方式。谢谢！
带有“sugar”和类似 R 语法的 Rcpp 并不难学。如果您已经有一些 c++ 经验，那么学习 Rcpp 可能会比您想象的要容易得多。（顺便说一句：编码“最佳实践”通常不鼓励使用 c 作为变量名......它有效，但是......）
谢谢，这令人鼓舞。是的，我只使用c 来回答这个问题。我的代码中的实际变量名是ln.pre.grid，但谢谢:)
无法理解“输出向量的第 k 个条目中的滞后总和”的含义。该代码显然是不正确的，并且没有真正的帮助。请指定一个小的可重现数据集和“正确”答案。

标签： r vector vectorization

【解决方案1】：

由于在 cmets 中提到了有关性能和 C/C++ 的内容，因此这里有一种使用 .Call 的方法，似乎有效：

library(inline)

ff = cfunction(sig = c(R_x = "numeric", R_lag = "integer"), body = '
   SEXP x, lag, ans;
   PROTECT(x = coerceVector(R_x, REALSXP));
   PROTECT(lag = coerceVector(R_lag, INTSXP));
   PROTECT(ans = allocVector(REALSXP, LENGTH(lag)));

   double *px = REAL(x), *pans = REAL(ans);
   memset(pans, 0, sizeof(double)*LENGTH(ans));
   R_len_t *plag = INTEGER(lag);

   for(int l = 0; l < LENGTH(lag); l++) 
       for(int i = 0; i < LENGTH(x) - plag[l]; i++) 
           pans[l] += px[i + plag[l]] - px[i];

   UNPROTECT(3);

   return(ans);
')

ff(1:100, 1:10)
#[1]  99 196 291 384 475 564 651 736 819 900

还有一些基准测试：

OPff = function(c, n) {
   sum.diff <- integer(n)
   for (k in 1:n) sum.diff[k] <- sum(diff(c, lag = k))
   sum.diff
}

ff2 = function(c, n) unlist(lapply(1:n, function(i) sum(diff(c, lag = i))))

xx = runif(1e4)
l = 1e3

identical(OPff(xx, l), ff(xx, 1:l))
#[1] TRUE
identical(OPff(xx, l), ff2(xx, l))
#[1] TRUE
library(microbenchmark)
microbenchmark(OPff(xx, l), ff(xx, 1:l), ff2(xx, l), times = 10)
#Unit: milliseconds
#        expr       min        lq    median        uq       max neval
# OPff(xx, l) 387.49171 390.43269 407.25796 427.09764 485.97181    10
# ff(xx, 1:l)  37.73505  38.27028  39.10201  41.33271  46.84648    10
#  ff2(xx, l) 384.35098 389.70397 401.51451 423.38513 436.85008    10

【讨论】：

【解决方案2】：

尝试以下方法：

 sum.diff <- function(c, n) sapply(n, function(k) sum(diff(c, lag = k)))

现在运行测试：

 sum.diff(1:100, 1:10)
 ## [1]  99 196 291 384 475 564 651 736 819 900

【讨论】：

这个向量化或与for 循环有何不同
@David Arenburg，它不涉及索引或显式循环。
它仍然不是向量化的 imo，只是以更好的语法呈现
最终所有操作都由最低级别的循环完成，无论方法如何。通常矢量化是指对整个对象进行操作而不使用索引。

【解决方案3】：

data.table 实现（应该比您在大数据集上的代码稍快）

a <- 1:100
b <- 1:10
library(data.table)
DT <- data.table(b)[, Res := sum(diff(a, b)), by = b]
DT

# b Res
# 1:  1  99
# 2:  2 196
# 3:  3 291
# 4:  4 384
# 5:  5 475
# 6:  6 564
# 7:  7 651
# 8:  8 736
# 9:  9 819
# 10: 10 900

【讨论】：