You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec <- function(x, alpha){
n <- length(x)
y <- numeric(n)
if (n == 0) return(y)
y[1] <- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i] <- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao <mkao006rm...@gmail.com> wrote:
>
> 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@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.
>>>
>>
>
> ______________________________________________
> R-help@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.
______________________________________________
R-help@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.
