Hi Is there a way to link a frequently-used sub-routine (in binary form or preferably in ascii form) in a program, similar to the following in Gauss? Thanks.
#include c:\programs\mylib1.g; -- Steven T. Yen, Professor of Agricultural Economics The University of Tennessee http://web.utk.edu/~syen/ ________________________________________ From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of David Winsemius [dwinsem...@comcast.net] Sent: Sunday, November 27, 2011 11:46 AM To: R. Michael Weylandt <michael.weyla...@gmail.com> Cc: r-help@r-project.org Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1}) 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 <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. 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. ______________________________________________ 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.