Можно ли векторизовать последовательное обновление элементов вектора в R?

Question

Можно ли оцифровать код, как указано ниже?

length(x) <- 100; 
x[1]  <- 1; 
y   <- rnorm(100); 

for(i in 2:100) { 
    x[i] <- 2 * y[i] * x[i-1]; 
} 

Я ценю, что это тривиальный пример, но он служит для иллюстрации идеи.

Мне часто приходится писать код, где i-е значение в векторе зависит от (i-1) -го значения, и, если возможно, я бы хотел написать это без необходимости в цикле for, поскольку профилирование предлагает функции с этим типом операций являются основными узкими местами в моем коде.

Является ли эта операция векторируемой, поэтому мне не нужно использовать цикл for() в расчете?

Answer 1

В общем, если вы хотите получить векторное решение, вам необходимо решить recurrence relation.

Answer 2

В примере вы можете найти формулу для x [i] и посмотреть, можно ли ее векторизовать. В этом случае я думаю, что cumprod может работать.

x <- c(1, cumprod(2*y)[1:99]) 

В некоторых случаях случае Вы можете также использовать команду filter в свертке или рекурсивном режиме. См ?filter

Однако если это не представляется возможным выработать формулу для значения n'th, что соответствует одной из форм выше, вы можете попробовать использовать пакет как inline или Rcpp, чтобы написать это в цикле в C/C++.

Answer 3

Внутренняя часть этой сюжетной команды эквивалентна. Довольно интересно повторно запустить его:

участок (с (1, 2^(2: длина (х) -1) * cumprod (RNorm (99))))

Answer 4

У меня нет полной информации по-прежнему, но выглядит, что функция filter() будет полезна для того, чтобы делать то, что мне нужно.

Answer 5

Вы можете написать не-vertorized код в C++:

library(inline) 
myfun <- cxxfunction(signature(y="numeric"), body=' 
Rcpp::NumericVector yvec(y); 
int ysize = yvec.size(); 
Rcpp::NumericVector result(ysize); 
if (ysize > 0) { 
    result[0] = 1; 
    for (int i = 1; i < ysize; i++) { 
     result[i] = 2 * yvec[i] * result[i-1]; 
    } 
} 
return result; 
', plugin="Rcpp") 

Затем вызовите эту функцию из R:

y <- rnorm(100); 
x <- myfun(y);