Example from Pineda-Krch08 J of Statistical Software.

```
rm(list=ls())
library(GillespieSSA)
parms = c(b=2, d=1, K=1000, alpha=0.007, w=0.0035, c=2,g=2)
x0 = c(N=1000, P=100)
nu = matrix(c(+1,-1,-1, 0,0,
0, 0, 0,+1,-1), nrow=2, byrow=TRUE)
a = c("b*N",
"(d+(b-d)*N/K)*N",
"alpha/(1+w*N)*N*P",
"c*alpha/(1+w*N)*N*P",
"g*P")
tf = 100
method = "D"
simName = "Predator-prey model"
```

`parms`

```
## b d K alpha w c g
## 2.0e+00 1.0e+00 1.0e+03 7.0e-03 3.5e-03 2.0e+00 2.0e+00
```

`nu `

```
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 -1 0 0
## [2,] 0 0 0 1 -1
```

`a`

```
## [1] "b*N" "(d+(b-d)*N/K)*N" "alpha/(1+w*N)*N*P"
## [4] "c*alpha/(1+w*N)*N*P" "g*P"
```

```
out = ssa(x0, a, nu, parms, tf, method, simName,
verbose = TRUE, consoleInterval = 10, maxWallTime = 30)
```

```
## Running D method with console output every 10 time step
## Start wall time: 2015-12-13 20:44:52...
## t=0 : 1000,100
## (4.082s) t=10 : 639,313
## (10.931s) t=20.00066 : 123,499
## (18.091s) t=30.00005 : 112,193
## (26.468s) t=40.00018 : 271,95
## t=43.15012 : 374,361
## tf: 43.15012
## TerminationStatus: maxWallTime
## Duration: 30.019 seconds
## Method: D
## Nr of steps: 96521
## Mean step size: 0.0004470542+/-0.0006694862
## End wall time: 2015-12-13 20:45:22
## --------------------
```

`summary(out$data)`

```
## V1 N P
## Min. : 0.000 Min. : 46.0 Min. : 9
## 1st Qu.: 9.575 1st Qu.: 294.0 1st Qu.: 92
## Median :19.063 Median : 496.0 Median :206
## Mean :20.357 Mean : 487.3 Mean :237
## 3rd Qu.:31.002 3rd Qu.: 677.0 3rd Qu.:361
## Max. :43.150 Max. :1014.0 Max. :622
```

`plot( out$data[,2] ~ out$data[,1], ylab = "N")`

`plot( out$data[,3] ~ out$data[,1], ylab = "Predator")`

## No comments:

## Post a Comment