14

我想会有更多的人对这个主题感兴趣。我有一些具体的任务要以最有效的方式完成。我的基础数据是： - 买入和卖出信号的时间指数 - 在时间指数上，我在最接近的买卖对之间拥有 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)))
}

你知道更有效的方法吗？如果这不能处理所有事情，或者更正确。

4

4 回答 4

13

如果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是一个标量，而 let z <- (1- tax/100)^2。您的最终结果r乘以z不同幂的raised矩阵。您要避免的是一遍又一遍地计算这些幂。这是我要做的：

pow <- 1L + col(r) - row(r)
pow[lower.tri(pow)] <- NA
tax.mult <- (z^(1:n))[pow]
r <- r * tax.mult
于 2013-06-19T23:20:24.383 回答
9

我采用了一种不同的方法，归结为使用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)))))
)
于 2013-06-19T23:44:15.380 回答
6

与矢量化方向不同，这是一种产生速度增益的方法（对于小型阵列来说非常大，对于大型阵列来说可以达到 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

关于 20130709 编辑：

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

于 2013-06-20T15:36:44.193 回答
1

免责声明：我在另一个答案中使用了这个。所以这将是一个无耻的插件。

要回答似乎是您的通用问题而不是您引用的示例 --- 如何将 for 循环转换为矢量化解决方案 --- 以下可能是一些有用的指针：

考虑您正在迭代的对象的结构。可能有不同的类型，例如：

a) 向量/矩阵的元素。b) 矩阵的行/列。c) 更高维数组的维度。d) 列表的元素（其本身可能是上述对象之一）。e) 多个列表/向量的对应元素。

在每种情况下，您使用的功能可能略有不同，但使用的策略是相同的。此外，学习应用家庭。各种 *pply 函数基于类似的抽象，但它们作为输入的内容和作为输出的内容有所不同。

例如，在上面的案例列表中。

a) 向量的元素：寻找已经存在的向量化解决方案（如上所示），这是 R 的核心优势。除此之外，还要考虑矩阵代数。大多数似乎需要循环（或嵌套循环）的问题都可以写成矩阵代数中的方程。

b) 矩阵的行/列：使用 apply。为 MARGIN 参数使用正确的值。c) 与更高维数组类似。

d) 使用 lapply。如果您返回的输出是“简单”结构（标量或向量），您可以考虑 sapply，它只是简单的 simple2array(lapply(...)) 并返回适当维度的数组。

e) 使用映射。'm' 可以代表多变量应用。

一旦你了解了你正在迭代的对象和相应的工具，就可以简化你的问题。不要想着你正在迭代的整个对象，而是它的一个实例。例如，当迭代矩阵的行时，忘记矩阵并只记住行。

现在，编写一个函数（或 lambda），它只对你的 iterand 的一个实例（元素）进行操作，并使用 *pply 系列的正确成员简单地“应用”它。

这是我使用cumprod. 这达到了大约 1000 x 1000 矩阵的最佳位置，但它返回一个列表而不是您期望的矩阵。但是，我没有将此作为解决方案提供，因为我认为您在基础 R 中的解决方案最好遵循 Rcpp 中的@eddi。这只是我上面讨论的过程的一个示例：

asb <- function (d) lapply(X=seq.int(from=length(d), to=1),
                           FUN=function (k) cumprod(d[seq_len(k)]))
于 2013-07-10T00:12:52.973 回答