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

William Dunlap wdunlap at tibco.com
Sun Nov 27 20:12:46 CET 2011 

The following quickly does the recursion suggested
by Florent D.:

  loopRec5 <- function(x, alpha) {
      as.vector(filter(x*alpha, alpha, method="recursive"))
  }

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Florent D.
> Sent: Sunday, November 27, 2011 10:01 AM
> To: David Winsemius
> Cc: r-help at r-project.org
> Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
> 
> ... 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.
> >
> 
> ______________________________________________
> 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