[R-sig-hpc] Avoiding looping on vector with stochastic dependency

Martin Morgan mtmorgan at fhcrc.org
Sun Dec 28 19:37:34 CET 2008


Hi Guillaume --

Guillaume Chapron wrote:
> 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



More information about the R-sig-hpc mailing list