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

[Michael Kao]

Hi Florent,

That is great, I was working on solving the recurrence equation and this 
was part of that equation. Now I know how to put everything together, 
thanks for all the help e everybody!

Cheers,

On 27/11/2011 2:05 p.m., Florent D. wrote:
> 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<mkao006rmail at 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 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.
>>>>
>> ______________________________________________
>> 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.