Sunday, December 13, 2015

GillespieSSA test. Single species

GillespieSSA test. Single species
h qin
December 13, 2015
Example from Pineda-Krch08 J of Statistical Software.
rm(list=ls())
library(GillespieSSA)
parms = c(b=2, d=1, K=1000)
x0 = c(N=500)
nu = matrix(c(+1,-1), ncol=2)
a = c("b*N", "(b+(b-d)*N/K)*N")
tf = 10
method = "D"
simName = "Logistic growth"
parms
##    b    d    K 
##    2    1 1000
nu 
##      [,1] [,2]
## [1,]    1   -1
a
## [1] "b*N"             "(b+(b-d)*N/K)*N"
out = ssa(x0, a, nu, parms, tf, method, simName, verbose = TRUE, consoleInterval = 1)
## Running D method with console output every 1 time step
## Start wall time: 2015-12-13 11:55:55...
## t=0 : 500
## (0.204s) t=1.000299 : 370
## (0.343s) t=2.000827 : 256
## (0.439s) t=3.001757 : 214
## (0.539s) t=4.001413 : 194
## (0.624s) t=5.003206 : 164
## (0.694s) t=6.001344 : 128
## (0.762s) t=7.000501 : 138
## (0.825s) t=8.000797 : 120
## (0.884s) t=9.00031 : 107
## t=10.0004 : 83
## tf: 10.0004
## TerminationStatus: finalTime
## Duration: 0.95 seconds
## Method: D
## Nr of steps: 8381
## Mean step size: 0.001193222+/-0.001501032
## End wall time: 2015-12-13 11:55:56
## --------------------
str(out)
## List of 3
##  $ data : num [1:8383, 1:2] 0 0.000869 0.003032 0.003286 0.003868 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:8383] "timeSeries" "" "" "" ...
##   .. ..$ : chr [1:2] "" "N"
##  $ stats:List of 8
##   ..$ startWallime      : chr "2015-12-13 11:55:55"
##   ..$ endWallTime       : chr "2015-12-13 11:55:56"
##   ..$ elapsedWallTime   : Named num 0.95
##   .. ..- attr(*, "names")= chr "elapsed"
##   ..$ terminationStatus : chr "finalTime"
##   ..$ nSteps            : int 8381
##   ..$ meanStepSize      : num 0.00119
##   ..$ sdStepSize        : num 0.0015
##   ..$ nSuspendedTauLeaps: num 0
##  $ args :List of 18
##   ..$ x0                 : Named num 500
##   .. ..- attr(*, "names")= chr "N"
##   ..$ a                  : chr [1:2] "b*N" "(b+(b-d)*N/K)*N"
##   ..$ nu                 : num [1, 1:2] 1 -1
##   ..$ parms              : Named num [1:3] 2 1 1000
##   .. ..- attr(*, "names")= chr [1:3] "b" "d" "K"
##   ..$ tf                 : num 10
##   ..$ method             : chr "D"
##   ..$ tau                : num 0.3
##   ..$ f                  : num 10
##   ..$ epsilon            : num 0.03
##   ..$ nc                 : num 10
##   ..$ hor                : num NaN
##   ..$ dtf                : num 10
##   ..$ nd                 : num 100
##   ..$ ignoreNegativeState: logi TRUE
##   ..$ consoleInterval    : num 1
##   ..$ censusInterval     : num 0
##   ..$ verbose            : logi TRUE

##   ..$ simName            : chr "Logistic growth"

No comments:

Post a Comment