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