R矢量化数组数据操作

Question

我想会有更多的人对这个主题感兴趣.我有一些特定的任务要以最有效的方式完成.我的基础数据是: - 买入和卖出信号的时间指数 - 在时间指示的诊断上我有最近的买卖对之间的ROC(变化率):

r <- array(data = NA, 
           dim = c(5, 5), 
           dimnames = list(buy_idx = c(1,5,9,12,16), 
                           sell_idx = c(3,7,10,14,19)))
diag(r) <- c(1.04,0.97,1.07,1.21,1.1)

任务是在每个可能的窗口(买卖对)上生成移动复合ROC,以及我目前正在解决我的任务的方式:

for(i in 2:5){
  r[1:(i-1),i] <- r[1:(i-1),i-1] * r[i,i]
}

直到我没有在上面的某个地方循环,我的解决方案的时间是非常可接受的.有没有办法将此循环更改为矢量化解决方案？是否有任何良好的文档化教程来学习R中的矢量化思维类型？ - 它比一次性解决方案更有价值!

编辑20130709:

下一个任务与先前的任务/示例高度相关.对每笔交易应用税额(税率为%值).当前解决方案

diag(r[,]) <- diag(r[,]) * ((1-(tax/100))^2)
for(i in 2:dim(r)[2]){
  r[1:(i-1),i] <- r[1:(i-1),i] * ((1-(tax/100))^(2*(i:2)))
}

你知道更有效的方法吗？或者更正确,如果这不能解决所有问题.

Answer 1

如果d是你的对角线元素,然后到处都是j >= i,r[i,j]是prod(d[i:j]),它也可以写prod(d[1:j]) / prod(d[1:(i-1)]).因此这个技巧使用outer累积产品的比例:

d <- c(1.04,0.97,1.07,1.21,1.1)
n <- length(d)
p <- cumprod(c(1, d))
r <- t(outer(p, 1/p, "*"))[-n-1, -1]
r[lower.tri(r)] <- NA

一些基准测试显示,对于某些(并非所有)输入大小,它确实优于OP:

OP <- function(d) {
   r <- diag(d)
   for(i in 2:length(d)){
     r[1:(i-1),i] <- r[1:(i-1),i-1] * r[i,i]
   }
   r
}

flodel <- function(d) {
   n <- length(d)
   p <- cumprod(c(1, d))
   r <- t(outer(p, 1/p, "*"))[-n-1, -1]
   r[lower.tri(r)] <- NA
   r
}

d <- runif(10)
microbenchmark(OP(d), flodel(d))
# Unit: microseconds
#        expr     min       lq   median      uq     max
# 1 flodel(d)  83.028  85.6135  88.4575  90.153 144.111
# 2     OP(d) 115.993 122.0075 123.4730 126.826 206.892

d <- runif(100)
microbenchmark(OP(d), flodel(d))
# Unit: microseconds
#        expr      min       lq    median       uq      max
# 1 flodel(d)  490.819  545.528  549.6095  566.108  684.043
# 2     OP(d) 1227.235 1260.823 1282.9880 1313.264 3913.322

d <- runif(1000)
microbenchmark(OP(d), flodel(d))
# Unit: milliseconds
#        expr      min        lq    median        uq       max
# 1 flodel(d) 97.78687 106.39425 121.13807 133.99502 154.67168
# 2     OP(d) 53.49014  60.10124  72.56427  85.17864  91.89011

编辑回答20130709加法:

我假设tax是一个标量而且让z <- (1- tax/100)^2.您的最终结果r乘以z不同权力的筹码矩阵.你想要避免的是一遍又一遍地计算这些权力.这是我要做的:

pow <- 1L + col(r) - row(r)
pow[lower.tri(pow)] <- NA
tax.mult <- (z^(1:n))[pow]
r <- r * tax.mult

Answer 2

我采取了一种不同的方法,归结为使用Reduce.给出一个简单的Reduce递归计算示例可能对某些人来说是值得的:

OP的预期结果:

> r
       sell_idx
buy_idx    3      7       10       14       19
     1  1.04 1.0088 1.079416 1.306093 1.436703
     5    NA 0.9700 1.037900 1.255859 1.381445
     9    NA     NA 1.070000 1.294700 1.424170
     12   NA     NA       NA 1.210000 1.331000
     16   NA     NA       NA       NA 1.100000

使用对角线起始值的基本示例 Reduce

x <- c(1.04,0.97,1.07,1.21,1.1)
Reduce(prod, tail(x,-1), x[1], accumulate=TRUE)

## gives first row of the answer 
## 1.04 / (1.04*0.97) / 1.07 * (1.04*0.97) etc etc etc

[1] 1.040000 1.008800 1.079416 1.306093 1.436703

循环起始值的长度并添加一些NA会得到完整的结果:

t(
  sapply(1:length(x),
    function(y) c(rep(NA,y-1),Reduce(prod, tail(x,-y), x[y], accumulate=TRUE))
    )
)

完整的结果:

     [,1]   [,2]     [,3]     [,4]     [,5]
[1,] 1.04 1.0088 1.079416 1.306093 1.436703
[2,]   NA 0.9700 1.037900 1.255859 1.381445
[3,]   NA     NA 1.070000 1.294700 1.424170
[4,]   NA     NA       NA 1.210000 1.331000
[5,]   NA     NA       NA       NA 1.100000

编辑

而且由于上面的Reduce幻想只是等同于cumprod,另一个更简单的解决方案就是:

rbind(
  cumprod(x),
  t(sapply(1:(length(x)-1),function(y) c(rep(NA,y),cumprod(tail(x,-y)))))
)

Answer 3

从矢量化开始走向不同的方向,这是一种产生速度增益的方法(对于小型阵列非常大,对于大型阵列则达到2-3倍范围):

library(inline)
library(Rcpp)

solver_fn = cxxfunction(signature(x = "numeric"), '
  NumericVector diag(x);

  unsigned n = diag.size();
  std::vector<double> result(n*n);

  result[0] = diag[0];

  unsigned col_shift_old = 0, col_shift = 0;
  for (unsigned col = 1; col < n; ++col) {
    col_shift = col * n;
    for (unsigned row = 0; row <= col; ++row) {
      if (result[row + col_shift_old] == 0)
        result[row + col_shift] = diag[col];
      else
        result[row + col_shift] = result[row + col_shift_old] * diag[col];
    }
    col_shift_old = col_shift;
  }

  return NumericVector(result.begin(), result.end());
', plugin = "Rcpp")

compute_matrix = function(d) {
  matrix(solver_fn(d), ncol = length(d))
}

这里有一些基准:

op = function(d) {
  r = diag(d)
  for (i in 2:length(d)) {
    r[1:(i-1), i] <- r[1:(i-1), i-1] * r[i,i]
  }
  r
}

d = runif(1e4)
system.time(op(d))
# user  system elapsed
#3.456   1.006   4.462
system.time(compute_matrix(d))
# user  system elapsed
#1.001   0.657   1.660

d = runif(1e3)
system.time(op(d))
# user  system elapsed
# 0.04    0.00    0.04
system.time(compute_matrix(d))
# user  system elapsed
#0.008   0.000   0.009

d = runif(1e2)
system.time(for (i in 1:1000) {op(d)})
# user  system elapsed
#1.075   0.000   1.075
system.time(for (i in 1:1000) {compute_matrix(d)})
# user  system elapsed
#0.075   0.000   0.075

Re 20130709编辑:

只需传递tax给C++函数并在那里进行乘法运算.如果您了解上述工作原理,那么更改将是微不足道的.