[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