Для l oop векторизованный

Question

Как я могу сделать что-то подобное, но оптимальным (векторизованным) способом в R?

N=10000

f <- 1.005
S0 <- 100
p <- 1/10

n <- seq(3,N)
S <- c(f*S0, f^2*S0, f^3*S0)
P <- c(0, 0, p*(f-1)*f^2*S0)
for(i in n){
  R <- tail(S,1)-tail(P,1)
  S <- c(S, f*R)
  P <- c(P, p*(f-1)*R)
}

конечный желаемый результат, конечно же, S и P (вплоть до ряд N + 1). Это вычисляет последовательный временной ряд строка за строкой (каждая строка является функцией значений предыдущей строки, над строкой 3).

Я пытался использовать lapply, но сложно заставить функцию возвращать два изменения в глобальной среде ... (и итоговая таблица также плохо отформатирована)

Answer 1

Самый простой шаг для ускорения вашего кода - это предварительно выделить векторы. Начните S и P с их окончательной длиной, а не «увеличивайте» их на каждой итерации l oop. Это приводит к более чем 100-кратному ускорению вашего кода:

N <- 10000
f <- 1.005
S0 <- 100
p <- 1/10

original = function(N, f, S0, p) {
  n <- seq(3,N)
  S <- c(f*S0, f^2*S0, f^3*S0)
  P <- c(0, 0, p*(f-1)*f^2*S0)
  for(i in n){
    R <- tail(S,1)-tail(P,1)
    S <- c(S, f*R)
    P <- c(P, p*(f-1)*R)
  }
  return(list(S, P))
}

pre_allocated = function(N, f, S0, p) {
  n <- seq(3,N)
  S <- c(f*S0, f^2*S0, f^3*S0, rep(NA, N - 3))
  P <- c(0, 0, p*(f-1)*f^2*S0, rep(NA, N - 3))
  for(i in n){
    R <- S[i] - P[i]
    S[i + 1] <- f*R
    P[i + 1] <- p*(f-1)*R
  }
  return(list(S, P))
}


## Check that we get the same result
identical(original(N, f, S0, p), pre_allocated(N, f, S0, p))
# [1] TRUE

## See how fast it is
microbenchmark::microbenchmark(original(N, f, S0, p), pre_allocated(N, f, S0, p), times = 10)
# Unit: milliseconds
#                        expr      min       lq      mean    median       uq      max neval
#       original(N, f, S0, p) 414.3610 419.9241 441.26030 426.01610 454.6002 538.0523    10
#  pre_allocated(N, f, S0, p)   2.3306   2.6478   2.92908   3.05785   3.1198   3.2885    10

Возможно, что векторизованное решение, возможно с использованием такой функции, как cumprod, будет еще быстрее, но я не вижу четкий способ сделать это. Если вы можете математически записать свой результат в виде совокупной суммы или произведения, это сделает его более ясным и, возможно, откроет решение.