Michael Kao
2011-Nov-27 00:20 UTC
[R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
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.
Enrico Schumann
2011-Nov-27 09:43 UTC
[R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
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.
>
--
Enrico Schumann
Lucerne, Switzerland
http://nmof.net/