Implementing the stochastic simulation algorithm in R
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Dear blog,
no I have not forgot you or the rest of the world. Sequestered in my humble cubicle I have been hammering away day and night for the last several weeks – hammering out a manuscript. Finally after what feels like an eternity I am almost there, i.e. getting close to submitting it. I just send a draft out for an in-lab review, so hopefully it is just a matter of some final polishing before I can get it out of my hair. If curiosity got the best of you, here’s the (draft) abstract.
Abstract
The deterministic dynamics of populations in continuous time is traditionally described using coupled, first-order ordinary differential equations. While this approach is accurate for large systems it is often inadequate for small systems where key species may be present in small numbers or where key reactions occur at a low rate. The stochastic simulation algorithm (SSA) is a procedure for generating time-evolution trajectories of finite populations in continuous time and has become the standard algorithm for these type of models. This article presents a simple to use and flexible framework for implementing the SSA using the high-level statistical computing language R. Using three biological models as examples, logistic growth, Rosenzweig-MacArthur predator-prey model, and Kermack-McKendrick SIRS metapopulation model, I show how the deterministic model can be formulated as a finite-population stochastic model within the framework of SSA theory and I give examples of it’s implementation in R. The example models are run using different SSA Monte Carlo methods, one exact method (Gillespie’s Direct method), and three approximate methods (Explicit, Binomial, and Optimized tau-leap methods). Comparison of simulation results from the different methods confirm that the time-evolution trajectories are indistinguishable while the approximate methods outperform the exact method by up to four orders of magnitude.
… yes, four orders of magnitude. I could hardly believe it myself when I saw it. For those of you not familiar with theoretical mumbojumbo, it means it is on the order of 10,000 faster. Actually, even that is an under estimate. The record I have is one of the accelerated methods running 63,000 times faster than the exact method (yes, that’s for exactly the same stochastic model)…, while still giving results that are indistinguishable from the results of the exact method.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.