Ответ 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

Ответ 2

Я использовал другой метод, который сводится к использованию Reduce. Например, простой пример Reduce для рекурсивных вычислений может быть полезен кому-то:

ожидаемый результат:

> 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 fanciness просто эквивалентно cumprod, альтернативное более простое решение было бы просто:

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

Ответ 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++ и сделайте там умножения. Если вы понимаете, как это работает, изменение будет тривиальным.

Ответ 4

Отказ от ответственности: я использовал это в другом ответе. Так что это будет бесстыдный плагин.

Чтобы ответить на то, что кажется вашим общим вопросом, а не на пример, который вы указали, - как преобразовать цикл for в векторное решение --- следующие могут быть несколько полезных указателей:

Рассмотрим структуру объекта, который вы повторяете. Могут быть разные типы, например:

a) Элементы вектора/матрицы. б) Строки/столбцы матрицы. c) Размерность многомерного массива. d) Элементы списка (который внутри себя может быть одним  указанных выше объектов). e) Соответствующие элементы множества списков/векторов.

В каждом случае используемая вами функция может немного отличаться, но стратегия использования одинакова. Кроме того, изучите применимую семью. Различные функции * pply основаны на аналогичной абстракции, но отличаются тем, что они принимают в качестве входных данных и что они выбрасывают в качестве вывода.

В приведенном выше списке случаев, например.

a) Элементы вектора: Ищите уже существующие векторизованные решения (как указано выше), которые являются основной силой в R. В дополнение к этому рассмотрим матричную алгебру. Большинство проблем, которые, как представляется, требуют циклов (или вложенных циклов), могут быть записаны как уравнения в матричной алгебре.

b) Строки/столбцы матрицы: Используйте. Используйте правильное значение для аргумента MARGIN. Similary для c) для массивов с большими размерами.

d) Используйте лапку. Если результат, который вы возвращаете, является "простой" структурой (скалярным или векторным), вы можете рассмотреть sapply, который просто упрощает2array (lapply (...)) и возвращает массив в соответствующих измерениях.

e) Используйте mapply. "M" может выступать за многомерное применение.

Как только вы поняли объект, который вы итерируете, и соответствующий инструмент, упростите свою проблему. Подумайте не об общем объекте, который вы повторяете, а о одном его экземпляре. Например, при итерации по строкам матрицы забудьте о матрице и запомните только строку.

Теперь напишите функцию (или лямбда), которая работает только с одним экземпляром (элементом) вашего iterand и просто "примените" его, используя правильный член семейства * pply.

Вот моя попытка проблемы с помощью cumprod. Это попадает в сладкое пятно размером около 1000 х 1000 матриц, но оно возвращает список, а не матрицу, как вы ожидаете. Однако я не предлагаю это как решение, так как я думаю, что ваше решение в базе R лучше всего следует @eddi в Rcpp. Это всего лишь пример процесса, о котором я говорил выше:

asb <- function (d) lapply(X=seq.int(from=length(d), to=1),
                           FUN=function (k) cumprod(d[seq_len(k)]))