Avoiding looping on vector with stochastic dependency

Hi Guillaume --

Hello,

I need to perform a computation on a vector, where the value at index i 
is a stochastic function of the value at index i-1. A example would the 
following code:


w <- numeric(100)
w[1] <- 10

# set.seed(2)

for (i in 2:length(w)) {
        w[i] <- sum(rbinom(w[i-1],1,0.5)) + sum(rpois(w[i-1],0.5))
    }


In my program, this code would be executed many times, and I would like 
to have it running faster. I assume I could do this be removing the loop 
and using a vector operation. But, whatever the way I think, I cannot 
figure out how to remove the loop on the index. I would be grateful if 
you could tell me if this is possible?

This question really belongs on the R-help list; this list is for 
questions related to cluster and other parallel computing solutions.

The sum of binomial or Poisson random variables are themselves binomial 
or Poisson distributed, so that sum(rpois(...)) can be replaced by 
rpois(sum(...), ...), etc. This replaces several 'expensive' calls 
(rbinom, rpois) with just one.

The transitions you describe are time-invariant (independent of i), so 
that you can pre-compute, probably analytically but certainly 
numerically, a matrix T[i,j] describing the probability of transition to 
state i, given current state j. The 'simulation' of a single transition 
is to sample from the distribution described by column j.

The simulation describes a first-order Markov process. For this there 
are explicit solutions, e.g., for the stationary distribution (wikipedia 
Markov Chain gives an introduction). So for the problem you outline what 
you likely want to do is to calculate T, and solve for the properties 
you're interested in. This will most easily involve a combination of 
analytic (formulating the solution in an abstract way) and numeric 
(finding the solution for your particular parameters) approaches.

This is not my area of expertise, so better advice might well be available.

Martin

Thanks so much!!

Guillaume

_______________________________________________
R-sig-hpc mailing list
R-sig-hpc at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-hpc

Thanks for your reply. In fact my code example was simplified from  
something more general

for (i in 2:length(w[,1])) {
		w[i,] <- foo( w[i-1,] )
	}

and I was trying to write this without a loop. I will look at the  
method you describe, but  I'm not sure if it can still be applied if  
foo() is quite complex. (sorry for sending this to the wrong list, I  
thought this one was about everything that can make R running faster!)

Cheers

Guillaume

On Sun, Dec 28, 2008 at 12:48 PM, Guillaume Chapron

Thanks for your reply. In fact my code example was simplified from something
more general

for (i in 2:length(w[,1])) {
               w[i,] <- foo( w[i-1,] )
       }

and I was trying to write this without a loop. I will look at the method you
describe, but  I'm not sure if it can still be applied if foo() is quite
complex. (sorry for sending this to the wrong list, I thought this one was
about everything that can make R running faster!)

In general, it is not possible to write out such recurrence relations
explicitly in closed form - see (e.g.)
http://en.wikipedia.org/wiki/Recurrence_relation for more details.

Hadley

http://had.co.nz/