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

[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.