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 <mkao006rm...@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<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. ______________________________________________ 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.