0

Say I have a function func that takes two scalar numeric inputs and delivers a scalar numeric result, and I have the following code to calculate a result vector u, based on input numeric vector v and initial value u0 for the result vector:

u<-rep(u0,1+length(v))
for (k in 2:length(u)){
  u[k]<-func(u[k-1],v[k-1])
}

Note how a component of the result vector depends not only on the corresponding element of the input vector but also on the immediately prior element of the result vector. I can see no obvious way to vectorise this.

It is common to do this sort of thing in financial simulations, for instance when projecting forward company accounts, rolling them up with interest or inflation and adding in operational cash flows each year.

For some specific instances, it is possible to find a case-specific, non-iterative coding, but I would like to know if there's a general solution.

The problem can also be coded by recursion, as follows:

calc.u<-function(v,u0){
  if (length(v)<2){
    func(u0,v[1]) } 
  else {
    u.prior<-func(u0,v[-length(v),drop=FALSE])       
    c(u.prior,func(u.prior[length(u.prior)],v[length(v)]) )
  }
u<-calc.u(v,u0)

Is there an R tactic for doing this without using either iteration or recursion, ie for vectorising it?

Answered: Thank you @MrFlick for introducing me to the Reduce function, which does exactly what I was wanting. I see that

Reduce('+',v,0,accumulate=T)[-1]

gives me

cumsum(v)

and 

Reduce('*',v,0,accumulate=T)[-1]

gives me

cumprod(v)

as expected, where the [-1] is to discard the initial value.

Very nice indeed! Thanks again.

Andrew Kirk
  • 2,027
  • 2
  • 11
  • 16
  • I'm going to go out on a limb here and say, no, there isn't really a way to vectorise it generally with R code only. If you can use existing things like `cumsum` `cumprod` etc then I would do so. You could pass any specific functions off to more efficient looping code, maybe through Rcpp? – thelatemail Apr 12 '17 at 03:09
  • This would be easier to answer with a proper [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) because there's nothing we can really test with here. I think you should be able to use `Reduce(..., acc=TRUE)` to calculate such values. – MrFlick Apr 12 '17 at 03:10
  • 1
    I don't there is any way, fundamentally, to do it without iteration or recursion. As you said, the results depends on prior results, so how can you perform the computation in parallel (e.g., vectorized)? – thc Apr 12 '17 at 03:10
  • Are you hoping that not using loops will offer speed advantages? What do you think the benefits of vectorising will be? – Marius Apr 12 '17 at 03:20
  • @Marius. I have a vague expectation that it might be faster. But it's mostly just a feeling of unease about using loops in R, because it has been so thoroughly drummed into me that one should vectorise rather than loop wherever possible. So - kind of an irrational guilt thing, maybe? – Andrew Kirk Apr 12 '17 at 08:09

1 Answers1

2

If you have this example

u0 <- 5
v <- (1:5)*2
func <- function(u,v) {u/2+v}

u <- rep(u0,1+length(v))
for (k in 2:length(u)){
  u[k]<-func(u[k-1],v[k-1])
}

this is equivalent to 

w <- Reduce(func, v, u0, accumulate=TRUE)

And we can check that

all(u==w)
# [1] TRUE
MrFlick
  • 195,160
  • 17
  • 277
  • 295