ssa.exact {adaptivetau}  R Documentation 
Exact stochastic simulation algorithm
Description
Simulates the trajectory of a continuoustime Markov process (aka. the Gillespie simulation algorithm).
Usage
ssa.exact(init.values, transitions, rateFunc, params, tf,
reportTransitions = FALSE)
Arguments
init.values 
Vector of initial values for all variables. This should be a simple onedimensional vector of real numbers. IMPORTANT: variables must never take negative values. 
transitions 
Twodimensional matrix of integers specifying how each state variable (rows) should be changed for a given transition (columns). Generally this will be a sparse matrix of primarily 1s and 1s. 
rateFunc 
Rfunction that returns instantaneous transition rates for each
transition in the form a realvalued onedimensional vector with
length equal to the number of transitions. The order of these
rates must match the order in which transitions were specified in
the

params 
any R variable to be passed asis to rateFunc, presumably specifying useful parameters. 
tf 
final time of simulation. Due to the approximate nature of the tauleaping algorithm, the simulation may overshoot this time somewhat. 
reportTransitions 
default FALSE; whether to include a matrix of executed transitions per time interval in the output (alongside with the states at each time point). 
Details
This function is supplied as a bonus with the adaptivetau package, since the C++ function that underlies this (exact) stochastic simulation algorithm is used in the (approximate) adaptive tau stochastic simulation as well.
The initial values, transition matrix, and transition rates must all be designed such that variables are always nonnegative. The algorithm relies on this invariant.
Consider calling enableJIT(1)
before running
ssa.exact. In most circumstances, this should yield some speedup by
bytecode compiling the rate function.
Value
Twodimensional matrix with rows for every timepoint at which the rateFunc was evaluated and columns giving the value for every state variable at that time. The first column specifies the time.
If the ‘reportTransitions’ option is used, then a list is returned with two elements. The first element is the ‘dynamics’ matrix described above. The second is a twodimensional matrix called ‘transitions’ with a row for each timepoint and columns giving the number of times each transition was executed between the current time and the previous time.
Author(s)
Philip Johnson
See Also
This function is a bonus the comes along with the approximate (but
faster) ssa.adaptivetau
. For systems with sparse
transition matrices, see helper function ssa.maketrans
.
Examples
## LotkaVolterra example
lvrates < function(x, params, t) {
with(params, {
return(c(preygrowth*x["prey"], ## prey growth rate
x["prey"]*x["pred"]*eat, ## prey death / predator growth rate
x["pred"]*preddeath)) ## predator death rate
})
}
params=list(preygrowth=10, eat=0.01, preddeath=10);
r=ssa.exact(c(prey = 1000, pred = 500),
matrix(c(1,0, 2,1, 0,1), nrow=2), lvrates, params, tf=2)
matplot(r[,"time"], r[,c("prey","pred")], type='l', xlab='Time', ylab='Counts')
legend("topleft", legend=c("prey", "predator"), lty=1:2, col=1:2)