Reputation: 2273
I'm using the GillespieSSA package for R, which allows you to call a single instance of a stochastic simulation easily and quickly. But a single realization isn't super-useful, I'd like to be able to look at the variability that arises from chance, etc. Consider a very basic model:
library(GillespieSSA)
x0 <- c(S=499, I=1, R=0)
a <- c("0.001*{S}*{I}","0.1*{I}")
nu <- matrix(c(-1,0,
+1,-1,
0, +1),nrow=3,byrow=T)
out <- ssa(x0, a, nu, tf=100)
That puts out a profoundly complicated list, the interesting bits of which are in out$data
.
My question is, I can grab out$data for a single call of the instance, tag it with a variable indicating what call of the function it is, and then append that data to the old data to come out with one big set at the end. So in crude pseudo-R, something like:
nruns <- 10
for (i in 1:nruns){
out <- ssa(x0, a, nu, tf=100)
data <- out$data
run <- rep(i,times=length[data[,2]))
data <- cbind(data,run)
But where data isn't being overwritten at each time step. I feel like I'm close, but having hopped between several languages this week, my R loop fu, weak as it already is, is failing.
Upvotes: 1
Views: 1014
Reputation: 132576
I am not sure if I understood you correctly. Do you want to do something like the following?
out <- lapply(X=1:10,FUN=function(x) ssa(x0, a, nu, tf=100)$data)
This would do 10 runs and put the resulting data lists in a list. You could then access,e.g., the data from the second run with out[[2]]
.
Upvotes: 2