GillespieSSA 0.5-1 is released

Mario Pineda-Krch

14 years ago

[This article was first published on Mario's Entangled Bank » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I just uploaded GillespieSSA 0.5-1 to CRAN. Now it’s just a matter of days before it has propagated itself across all CRAN mirrors. This version consists primarily of revisions I made in response to the reviewer comments on the paper where the package is introduced (submitted to the Journal of Statistical Software). There are some minor changes in the functionality of the ssa.plot() function but otherwise the changes consists entirely of buglet fixes and improvements to the documentation.

One of the more interesting comments I got addressed the use of a character vector to pass the propensity functions to the wrapper function ssa(). For example, normally one would define the propensity functions for a logistic growth model as

a <- c("c1*Y1", "c2*Y1*Y2","c3*Y2")

This is the way it would be defined if the simulation is invoked using the higher level wrapper function ssa(). One can, however, also pass the propensity vector as a function by directly invoking the lower-level method function (ssa.d, ssa.btl, ssa.etl, ssa.otl, …). For the logistic growth model this could be done like so,

a = function(parms,x){
b <- parms[1]
d <- parms[2]
K <- parms[3]
N <- x[1]
return(c(b*N , N*b + (b-d)*N/K))
}
parms <- c(2,1,1000,500)
x <- 500
nu <- matrix(c(+1, -1),ncol=2)
t <- 0
for (i in seq(100)) {
out <- ssa.d(a(parms,x),nu)
x <- x + out$nu_j
t <- t + 1
cat("t:",t,", x:",x,"n")
}
)

The obvious advantage of this approach is that the propensity vector is simpler to define and maintain throughout the simulation. It is also likely that the simulation would run faster (or at least being simpler to optimize) without the extra over head imposed by the higher-level wrapper function. The disadvantage is that setting up the simulation is a tad more involved since one has to “manually” update the state vector, time variable and collect the output data (which, by the way, the above routine does not do).

Perhaps the most interesting consequence of directly invoking the lower-level method functions is that one now can vary the parameters of the model during the simulation allowing for temporal environmental heterogeneity, e.g. varying vital rates and carrying capacity over time.