arrays - 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 中的矢量化思维类型？ - 这将比一次性解决方案更有值(value)!

编辑 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