Is it possible to vectorise code like the following?
length(x) <- 100;
x[1] <- 1;
y <- rnorm(100);
for(i in 2:100) {
x[i] <- 2 * y[i] * x[i-1];
}
I appreciate that this is a trivial exam开发者_如何学JAVAple, but it serves to illustrate the idea.
I often need to write code where the i-th value in a vector depends on the (i-1)-th value and if possible, I'd like to write this without needing a for
loop, as profiling suggests the functions with this type of operation are the major bottlenecks in my code.
Is this operation vectorizable so I do not need to use a for()
loop in the calculation?
In general, if you want a vectorised solution you need to solve the recurrence relation.
In the example you have you could work out the formula for x[i] and see if it can be vectorized. In this case I think cumprod might work.
x <- c(1, cumprod(2*y)[1:99])
For some cases case you can also use the filter
command in convolution or recursive mode. See ?filter
However if it is isn't possible to work out a formula for the n'th value that fits one of the molds above, you could try using a package like inline
or Rcpp
to write this in loop in C/C++.
The interior of this plot command is equivalent. Rather interesting to repeatedly run it:
plot(c(1, 2^(2:length(x)-1)*cumprod(rnorm(99) )) )
You can write the non-vertorized code in 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")
Then call this function from R:
y <- rnorm(100);
x <- myfun(y);
I don't have full details on this yet, but it looks the function filter()
is going to be useful to do what I need.
精彩评论