# Stochastic SIR model with Julia

using RandomNumbers
using DataFrames

@inline @fastmath function randbn(n,p,rng)
q = 1.0 - p
s = p/q
a = (n+1)*s
r = exp(n*log(q))
x = 0
u = rand(rng)
while true
if (u < r)
return x
end
u -= r
x += 1
r *= (a/x)-s
end
end

randbn (generic function with 1 method)

@inline @fastmath function sir(u, parms, rng)
(S, I, R, Y) = u
(β, γ, ι, N, δt) = parms
λ = β * (I + ι) / N
ifrac = 1.0 - exp(-λ * δt)
rfrac = 1.0 - exp(-γ * δt)
infection = randbn(S, ifrac, rng)
recovery = randbn(I, rfrac, rng)
return (S - infection, I + infection - recovery, R + recovery, Y + infection)
end

sir (generic function with 1 method)

function simulate(r)
parms = (0.1, 0.05, 0.01, 1000.0, 0.1)
tf = 200
t = 0:0.1:tf
tl = length(t)
S = zeros(tl)
I = zeros(tl)
R = zeros(tl)
Y = zeros(tl)
u0 = (999, 1, 0, 0)
(S,I,R,Y) = u0
u = u0
for i in 2:tl
u = sir(u, parms, r)
(S[i],I[i],R[i],Y[i]) = u
end
return DataFrame(Time=t,S=S,I=I,R=R,Y=Y)
end

simulate (generic function with 1 method)

seed = 42
r = Xorshifts.Xorshift128Plus(seed);

sir_out = simulate(r);

head(sir_out)

TimeSIRY
Float64Float64Float64Float64Float64
10.0999.01.00.00.0
20.1999.01.00.00.0
30.2999.01.00.00.0
40.3999.01.00.00.0
50.4999.01.00.00.0
60.5999.01.00.00.0

#### Visualisation

using Plots, StatPlots

ArgumentError: Package StatPlots not found in current path:
- Run import Pkg; Pkg.add("StatPlots") to install the StatPlots package.




Stacktrace:

  require(::Module, ::Symbol) at ./loading.jl:820

  top-level scope at In:1

@df sir_out plot(:Time, [:S :I :R], colour = [:red :green :blue], xlabel="Time",ylabel="Number")

UndefVarError: @df not defined