On Nov 27, 2011, at 1:00 PM, Florent D. wrote:
... but the original request was to build a series, not approximate its limit by building a different (but "faster") series.
Right. Your function was faster that the earlier ones. But if speed were the issue, it might make more sense to use mathematics.
What I suggested is linear in terms of the size of x so you won't find a faster algorithm.
Actually, Dunlap's was 13-fold faster than yours on my machine.
> loopRec5 <- function(x, alpha) {
+ as.vector(filter(x*alpha, alpha, method="recursive"))
+ }
>
> loopRec4 <- 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)
+ }
>
>
> ## benchmark the functions.
> benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
+ loopRec3(1:1000, 0.5),loopRec4(1:1000, 0.5),loopRec5(1:1000, 0.5),
+ replications = 50, order = "relative")
test replications elapsed relative user.self
sys.self user.child
5 loopRec5(1:1000, 0.5) 50 0.017 1.00000 0.016
0.000 0
4 loopRec4(1:1000, 0.5) 50 0.228 13.41176 0.227
0.000 0
3 loopRec3(1:1000, 0.5) 50 1.458 85.76471 1.461
0.005 0
2 loopRec2(1:1000, 0.5) 50 1.579 92.88235 1.579
0.006 0
1 loopRec(1:1000, 0.5) 50 3.006 176.82353 3.006
0.012 0
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 <dwinsem...@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. (Ithink 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.[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812EMA(1:10, n=1, ratio=0.5)8.003906 [10] 9.001953[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906loopRec3(1:10, 0.5)8.001953 [10] 9.000977This 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 maybe a simple analytical approximation that makes the process even more efficient. The asymptotic result appears to be: n-(1-ratio)/ratiotail(EMA(1:1000, n=1, ratio=0.5))[1] 994 995 996 997 998 999tail(loopRec3(1:1000, 0.5))[1] 994 995 996 997 998 999EMA(1:1000, n=1, ratio=0.1)[491:500][1] 482 483 484 485 486 487 488 489 490 491 -- David.MichaelOn 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, thanksfor 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 ourimplementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:Hi Michaelhere are two variations of your loop, and both seem faster than theoriginal 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 gettest replications elapsed relative user.self sys.self2 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 giventhe 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 arejusttotally ridiculous so feel free to have a laugh but would appreciateifsomeone can give me a hint. Otherwise I guess I'll have to give RCppa 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.htmland 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.htmland 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.David Winsemius, MD West Hartford, CT ______________________________________________ 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.
David Winsemius, MD West Hartford, CT ______________________________________________ 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.
