# ssa: Gillespie Stochastic Simulation Algorithm In INSP-RH/ssar: A speedy implementation of Gillespie's Stochastic Simulation Algorithm

## Description

Implementation of Gillespie's Stochastic Simulation Algorithm for in-homogeneous continuous-time Markov Chains.

## Usage

 ```1 2 3 4``` ```ssa(xinit, pfun, v, params = c(), tmin = 0, tmax = Inf, nsim = 10, maxiter = Inf, maxtime = Inf, print.time = FALSE, plot.sim = TRUE, title = "", xlab = "", ylab = "", kthsave = 1, keep.file = FALSE, file.only = FALSE, fname = "Temporary_File_ssa.txt") ```

## Arguments

 `xinit` ( matrix ) Numeric named vector with initial variables. `pfun` ( vector ) Character vector of propensity functions where state variables correspond to the names of the elements in xinit. `v` ( matrix ) Propensity Matrix `params` ( list ) List of parameters including functions used by the problem that are not contained in base R Optional `tmin` ( double ) Initial .time for the simulation `tmax` ( double ) Final .time for the simulation `nsim` ( integer ) Number of simulations `maxiter` ( numeric ) Maximum number of iterations for model `maxtime` ( numeric ) Maximum .time (in seconds) before breaking the proces `print.time` ( boolean ) Value indicating whether to print current estimation .time `plot.sim` ( boolean ) Variable indicating whether to plot the results `title` ( string ) Title for plot `xlab` ( string ) X-axis label for plot `ylab` ( string ) Y-axis label for plot `kthsave` ( numeric ) Integer indicating iteration-lapse to save iterations (saves every kthsave iteration) `keep.file` ( boolean ) Instruction whether to keep temporary file created by the process `file.only` ( boolean ) Instruction whether to return the information as a txt file when `TRUE`. this makes the package run faster. `fname` ( string ) Name of file where information will be saved if `keep.file == T`

## Note

To include .time dependent functions please express them as matrix functions of t. For example `pfun <- function(t){cbind(sin(t),2*t)}`

Rodrigo Zepeda

Dalia Camacho

## References

INCLUDE REFERENCES HERE TO GILLESPIE SSA WITH TIME-DEPENDENT PARAMETERS

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143``` ```#EXAMPLE 1 #------------------------ #Logistic Growth set.seed(123) params <- c(b=2, d=1, k=1000) X <- matrix(c(N=500), nrow = 1) pfun <- function(t,X,params){ cbind(params[1] * X[,1], (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1]) } v <- matrix( c(+1, -1),ncol=2) tmin <- 0 tmax <- 1 nsim <- 5 simulation <- ssa(X, pfun, v, params, tmin, tmax, nsim, title = "Logistic Growth example", xlab = "Time", ylab = "N") #EXAMPLE 2 #------------------------ #Time dependent logistic growth set.seed(123) params <- c(b=2, d=1, k=1000) X <- matrix(c(N=500), nrow = 1) pfun <- function(t,X,params){ cbind(params[1] *(1 + sin(t))* X[,1], (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1]) } v <- matrix( c(+1, -1),ncol=2) tmin <- 0 tmax <- 1 nsim <- 5 simulation <- ssa(X, pfun, v, params, tmin, tmax, nsim, title = "Time dependent Logistic Growth example", xlab = "Time", ylab = "N") #EXAMPLE 3 #------------------------ #Lotka Volterra #Set seed set.seed(123) #Get initial parameters params <- c(a = 3, b = 0.01, c = 2) #Notice that X must be inputed as matrix X <- matrix(c(100, 100), ncol = 2) #Notice that pfun returns a matrix pfun <- function(t, X, params){ cbind(params[1]*t*X[,1] + 1, params[2]*X[,1]*X[,2], params[3]*X[,2]) } #Propensity matrix v <- matrix(c(+1,-1,0,0,+1,-1),nrow=2,byrow=TRUE) #Additional variables tmin <- 0 nsim <- 10 maxiter <- 300 #Run the program simulation <- ssa(X, pfun, v, params, tmin, nsim = nsim, maxiter = maxiter, print.time = TRUE, title = "Lotka-Volterra example", xlab = "Time", ylab = "Individuals") #EXAMPLE 4 #------------------------ #Time dependent SIS model set.seed(123) #Initial parameters k <- 24576.5529836797 delta <- 0.0591113454895868 + 0.208953907151055 gamma_ct <- 0.391237630231631 params <- c(k = k, delta = delta, gamma_ct = gamma_ct) X <- matrix(c(S = 1000, I = 40), ncol = 2) pfun <- function(t, X, params){ #Value to return matreturn <- matrix(NA, nrow = length(t), ncol = 6) #Create birth function lambda <- function(t){ return(4.328e-4 - (2.538e-7)*t - (3.189e-7)*sin(2 * t * pi/52) - (3.812e-7)*cos(2 * t * pi/52) ) } #Create death function mu <- function(t){ return(9.683e-5 + (1.828e-8)*t + (2.095e-6)*sin(2 * t * pi/52) - (8.749e-6)*cos(2 * t * pi/52))} #Create infectives function beta_fun <- function(t){ return( 0.479120824267286 + 0.423263042762498*sin(-2.82494252560096 + 2*t*pi/52) )} #Estimate values matreturn[,1] <- lambda(t)*(X[,1] + X[,2]) matreturn[,2] <- mu(t)*X[,1] matreturn[,3] <- beta_fun(t)*X[,1]*X[,2]/(1 + params[1]*X[,2]) matreturn[,4] <- mu(t)*X[,2] matreturn[,5] <- params[2]*X[,2] matreturn[,6] <- params[3]*X[,2] #Return return(matreturn) } v <- matrix(c(1,-1, -1, 0, 0, 1, 0, 0, 1, -1, -1, -1), nrow = 2, byrow = TRUE) tmin <- 0 tmax <- 0.5 nsim <- 100 #Simulate the values simulation <- ssa(X, pfun, v, params, tmin, tmax, nsim = nsim, print.time = FALSE, plot.sim = FALSE, kthsave = 1) #Plot using ggplot2 library ## Not run: library(ggplot2) ggplot(data = simulation, aes(x = Time, y = Var2, group=as.factor(Simulation))) + geom_line(aes(color = as.factor(Simulation))) + theme_bw() + theme(legend.position="none") + ggtitle(paste0("SIS example; Infected cases ", nsim, " simulations")) + xlab("Time") + ylab("Individuals") + geom_vline(xintercept = tmax, linetype = 2) ## End(Not run) #EXAMPLE 5 #------------------------------------ set.seed(123) #Start the parameters params <- c(NA) X <- matrix(c(N = 10), nrow = 1) #Notice the cbind forces to return matrix pfun <- function(t, X, params){ cbind(1.1 + sin(pi*t/0.01))*X[,1] } v <- matrix( c(+1), ncol=1) tmin <- 0 nsim <- 25 maxtime <- 1 #Maximum .time for model: 1 seconds and stop iterating simulation <- ssa(X, pfun, v, params, tmin, nsim = nsim, maxtime = maxtime, print.time = FALSE, title = "Start of an epidemic with time-dependent R0", xlab = "Time", ylab = "Individuals") ```

INSP-RH/ssar documentation built on May 7, 2019, 6:02 a.m.