[R] aplpy recursive function on a list

[Berend Hasselman]

On 26-01-2012, at 19:10, Berend Hasselman wrote:

> 
> On 26-01-2012, at 17:58, Brian Diggs wrote:
> 
>> On 1/25/2012 10:09 AM, patzoul wrote:
>>> I have 2 series of data a,b and I would like to calculate a new series which
>>> is z[t] = z[t-1]*a[t] + b[t] , z[1] = b[1].
>>> How can I do that without using a loop?
>> 
>> ........
> 
> I don't think so.
> 
> a <- c(2,4,3,5)
> b <- c(1,3,5,7)
> 
> z <- rep(0,length(a))
> z[1] <- b[1]
> for( t in 2:length(a) ) { z[t] <- a[t] * z[t-1] + b[t] }
> z
> 
> gives
> 
> [1]   1   7  26 137
> 
> and this agrees with a manual calculation.
> 
> You get a vector of length 5 as result. It should be of length 4 with your data.
> If you change the Reduce expression to this
> 
> u <- Reduce(function(zm1, coef) {coef[1] * zm1 + coef[2]},
>           Map(c, a[-1], b[-1]),
>           init = b[1], accumulate = TRUE)
> 
> then you get the correct result.
> 
>> u
> [1]   1   7  26 137


And the loop especially if byte compiled with cmpfun  appears to be quite a bit quicker.

Nrep <- 1000

tfrml.loop <- function(a,b) {
    z <- rep(0,length(a))
    z[1] <- b[1]
    for( t in 2:length(a) ) {
        z[t] <- a[t] * z[t-1] + b[t]
    }

    z
}

tfrml.rdce <- function(a,b) {
    u <- Reduce(function(zm1, coef) {coef[1] * zm1 + coef[2]},
               Map(c, a[-1], b[-1]),
               init = b[1], accumulate = TRUE)
    u
}

library(compiler)
tfrml.loop.c <- cmpfun(tfrml.loop)
tfrml.rdce.c <- cmpfun(tfrml.rdce)

z.loop <- tfrml.loop(a,b)
z.rdce <- tfrml.rdce(a,b)
all.equal(z.loop, z.rdce)

library(rbenchmark)

N <- 500
set.seed(1)
a <- runif(N)
b <- runif(N)

benchmark(tfrml.loop(a,b), tfrml.rdce(a,b), tfrml.loop.c(a,b), tfrml.rdce.c(a,b),
                replications=Nrep, columns=c("test", "replications", "elapsed"))

                test replications elapsed
1   tfrml.loop(a, b)         1000   2.665
3 tfrml.loop.c(a, b)         1000   0.554
2   tfrml.rdce(a, b)         1000   4.082
4 tfrml.rdce.c(a, b)         1000   3.143

Berend

R-2.14.1 (32-bits), Mac OS X 10.6.8.