[R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})

[Michael Kao]

Dear Enrico,

Brilliant! Thank you for the improvements, not sure what I was thinking 
using rev. I will test it out to see whether it is fast enough for our 
implementation, but your help has been SIGNIFICANT!!!

Thanks heaps!
Michael

On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
> Hi Michael
>
> here are two variations of your loop, and both seem faster than the 
> original loop (on my computer).
>
>
> require("rbenchmark")
>
> ## your function
> loopRec <- function(x, alpha){
>     n <- length(x)
>     y <- double(n)
>     for(i in 1:n){
>         y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>     }
>     y
> }
> loopRec(c(1, 2, 3), 0.5)
>
> loopRec2 <- function(x, alpha){
>     n <- length(x)
>     y <- numeric(n)
>     for(i in seq_len(n)){
>         y[i] <- sum(cumprod(rep.int(alpha, i)) * x[i:1])
>     }
>     y
> }
> loopRec2(c(1, 2, 3), 0.5)
>
> loopRec3 <- function(x, alpha){
>     n <- length(x)
>     y <- numeric(n)
>     aa <- cumprod(rep.int(alpha, n))
>     for(i in seq_len(n)){
>         y[i] <- sum(aa[seq_len(i)] * x[i:1])
>     }
>     y
> }
> loopRec3(c(1, 2, 3), 0.5)
>
>
> ## Check whether value is correct
> all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
> all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
>
>
> ## benchmark the functions.
> benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
>  loopRec3(1:1000, 0.5),
>  replications = 50, order = "relative")
>
>
> ... I get
>                    test replications elapsed relative user.self sys.self
> 2 loopRec2(1:1000, 0.5)           50    0.77 1.000000      0.76     0.00
> 3 loopRec3(1:1000, 0.5)           50    0.86 1.116883      0.85     0.00
> 1  loopRec(1:1000, 0.5)           50    1.84 2.389610      1.79     0.01
>
>
> Regards,
> Enrico
>
>
> Am 27.11.2011 01:20, schrieb Michael Kao:
>> Dear R-help,
>>
>> I have been trying really hard to generate the following vector given
>> the data (x) and parameter (alpha) efficiently.
>>
>> Let y be the output list, the aim is to produce the the following
>> vector(y) with at least half the time used by the loop example below.
>>
>> y[1] = alpha * x[1]
>> y[2] = alpha^2 * x[1] + alpha * x[2]
>> y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
>> .....
>>
>> below are the methods I have tried and failed miserably, some are just
>> totally ridiculous so feel free to have a laugh but would appreciate if
>> someone can give me a hint. Otherwise I guess I'll have to give RCpp a
>> try.....
>>
>>
>> ## Bench mark the recursion functions
>> loopRec <- function(x, alpha){
>> n <- length(x)
>> y <- double(n)
>> for(i in 1:n){
>> y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>> }
>> y
>> }
>>
>> loopRec(c(1, 2, 3), 0.5)
>>
>> ## This is a crazy solution, but worth giving it a try.
>> charRec <- function(x, alpha){
>> n <- length(x)
>> exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>> up.mat <- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
>> ", 0:(n - 1), ")", sep = ""),
>> paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
>> collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
>> colSums(up.mat * exp.mat)
>> }
>> vecRec(c(1, 2, 3), 0.5)
>>
>> ## Sweep is slow, shouldn't use it.
>> matRec <- function(x, alpha){
>> n <- length(x)
>> exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>> up.mat <- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
>> byrow = TRUE), 1,
>> c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
>> up.mat[lower.tri(up.mat)] <- 0
>> colSums(up.mat * exp.mat)
>> }
>> matRec(c(1, 2, 3), 0.5)
>>
>> matRec2 <- function(x, alpha){
>> n <- length(x)
>> exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>> up.mat1 <- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
>> up.mat2 <- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
>> up.mat <- up.mat1 * up.mat2
>> up.mat[lower.tri(up.mat)] <- 0
>> colSums(up.mat * exp.mat)
>> }
>>
>> matRec2(c(1, 2, 3), 0.5)
>>
>> ## Check whether value is correct
>> all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
>> all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
>> all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
>>
>> ## benchmark the functions.
>> benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 
>> 0.5),
>> matRec2(1:1000, 0.5), replications = 50,
>> order = "relative")
>>
>> Thank you very much for your help.
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>