Simulates the SIR continuous-time Markov chain exactly. Two reaction
channels act on integer counts (S, I, R):
a1 = beta S I / N taking (S, I) -> (S - 1, I + 1), and
a2 = gamma I taking (I, R) -> (I - 1, R + 1), with total
propensity a0 = a1 + a2. Given two independent uniforms u1, u2 on
(0, 1) the direct method draws the waiting time and the channel as
tau = -log(u1) / a0 and channel 1 if and only if
u2 a0 < a1, which is exact – no time discretisation is involved.
The chain stops when I = 0 (the absorbing state) or when the clock
passes T.
Value
List with estimate (final size), S, I,
R, N, final_size, attack_rate,
peak_I, peak_time, t_end, extinction_time,
n_events, n_infections, R0, beta,
gamma, T, seed, method.
Details
The uniform stream is the Lehmer minstd generator
s <- 48271 s mod (2^31 - 1) shared with every other arm of this
package, so a given seed reproduces the same trajectory in R and in
Python bit for bit.
References
Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry 81(25), 2340-2361. doi:10.1021/j100540a008
Examples
Sirstn(S0 = 990, I0 = 10, beta = 0.3, gamma = 0.1, T = 100)
#> $estimate
#> [1] 922
#>
#> $S
#> [1] 75
#>
#> $I
#> [1] 3
#>
#> $R
#> [1] 922
#>
#> $N
#> [1] 1000
#>
#> $final_size
#> [1] 922
#>
#> $attack_rate
#> [1] 0.922
#>
#> $peak_I
#> [1] 304
#>
#> $peak_time
#> [1] 30.69388
#>
#> $t_end
#> [1] 100
#>
#> $extinction_time
#> [1] NA
#>
#> $n_events
#> [1] 1837
#>
#> $n_infections
#> [1] 915
#>
#> $R0
#> [1] 3
#>
#> $beta
#> [1] 0.3
#>
#> $gamma
#> [1] 0.1
#>
#> $T
#> [1] 100
#>
#> $seed
#> [1] 1
#>
#> $method
#> [1] "Stochastic SIR, Gillespie direct method (Gillespie 1977)"
#>
#> attr(,"class")
#> [1] "morie_rich_result" "list"