[R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
David Winsemius
dwinsemius at comcast.net
Sun Nov 27 22:30:46 CET 2011
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
> <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.
>>
David Winsemius, MD
West Hartford, CT
More information about the R-help
mailing list