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

Florent D. flodel at gmail.com
Sun Nov 27 19:00:32 CET 2011 

... but the original request was to build a series, not approximate
its limit by building a different (but "faster") series.

What I suggested is linear in terms of the size of x so you won't find
a faster algorithm. What will make a difference is the implementation,
e.g. C loop versus R loop.



On Sun, Nov 27, 2011 at 11:46 AM, David Winsemius
<dwinsemius at comcast.net> wrote:
>
> On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
>
>> You also might look at EMA() in the TTR package for a C implementation. (I
>> think it matches your problem but can't promise)
>>
>
> It's pretty close for EMA(1:1000, n=1,  ratio=0.5) and 7 times faster.
>
>> EMA(1:10, n=1,  ratio=0.5)
>  [1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812
> 8.003906
> [10] 9.001953
>> loopRec3(1:10, 0.5)
>  [1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906
> 8.001953
> [10] 9.000977
>
> This series converges very quickly to n-1 when the ratio = 0.5 and to n-3
> when ratio = 0.25. It's already within 1% at the fifth iteration. There may
> be a simple analytical approximation that makes the process even more
> efficient. The asymptotic result appears to be:  n-(1-ratio)/ratio
>
>> tail(EMA(1:1000, n=1,  ratio=0.5))
> [1] 994 995 996 997 998 999
>> tail(loopRec3(1:1000, 0.5))
> [1] 994 995 996 997 998 999
>
>> EMA(1:1000, n=1,  ratio=0.1)[491:500]
>  [1] 482 483 484 485 486 487 488 489 490 491
> --
> David.
>
>
>> 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.
>>
>> ______________________________________________
>> 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.
>
> David Winsemius, MD
> West Hartford, CT
>
> ______________________________________________
> 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