[R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
R. Michael Weylandt <michael.weylandt@gmail.com>
michael.weylandt at gmail.com
Sun Nov 27 15:25:40 CET 2011
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise)
Michael
On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
> 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.
>
> ______________________________________________
> 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.
More information about the R-help
mailing list