```
using DifferentialEquations
using DiffEqSensitivity
using Random
using Distributions
using Turing
using DataFrames
using StatsPlots
```

## Exploring Julia and ODEs

This post is exploring the use of Julia, Turing, and ODEs and largely expands on the work from Simon Frost’s awesome epirecipes.

## Getting Started with Julia

Julia as a statistical programming language/environment has some really neat properties.

Speed, generally speed.

Turing also has some neat elements when it comes to probabilistic computing/programming. The key feature that interests me is the ability to switch between sampling approaches (e.g., Gibbs, Hamiltonian Monte Carlo, no u-turn sampling, sequential) while keeping the same syntax and model structure. While I absolutely love Stan and have spent many hours learning and the syntax (and it is actively developed, highly optimized, etc), I sometimes wish I could explore other sampling approaches due to strange posteriors (especially when it comes to times when SMC might be better).

## The Approach

This is a pretty vanilla post, in which we will first create some fake data from a known distribution and then fit said data. As always, we will bring in the required packages.

Now we can generate a simple three compartment model, with S (susceptibles), I, (infected), and R (recovered).

Additionally:

- On average 10 contacts are made per day
- The contact rate is 0.05 (probability of passing infection to any given contact)
- Those who are infected recover on average every 5 days
- Immunity lasts for a 180 days on average before people become susceptible again
- Population size of 6,000 individuals
- Standard mass-action approach

Something that looks like the following:

```
using Luxor
using MathTeXEngine
```

```
Drawing(200, 200, "sirs.png")
origin()
setline(10)
arrow(Point(-40,0),Point(0,0))
Luxor.arrow(Point(-0,0),Point(40,0))
Luxor.fontsize(15)
text(L"S",Point(-40,20), halign = :center)
Luxor.
text(L"I",Point(0,20), halign = :center)
Luxor.text(L"R",Point(40,20), halign = :center)
Luxor.fontsize(12)
text(L"β*c",Point(-20,10), halign = :center)
Luxor.text(L"γ",Point(20,10), halign = :center)
Luxor.text(L"δ",Point(0,-15), halign = :center)
Luxor.
= 30
loopx = 40
loopy = 6
adjx = -2
adjy arrow(
Luxor.Point(40,0) + Point(-adjx,adjy),
Point(40,0) + Point(-loopx,-loopy),
Point(-40,0)+ Point(loopx,-loopy),
Point(-40,0)+Point(adjx,adjy)
)finish()
```

```
function sir_ode!(du,u,p,t)
= u
(S,I,R,C) = p
(β,c,γ, δ) = S+I+R
N = β*c*I*S/N
infection = γ*I
recovery = δ * R
wane @inbounds begin
1] = -infection + wane
du[2] = infection - recovery
du[3] = recovery - wane
du[4] = infection
du[end
nothing
end;
= 365.0
tmax = (0.0,tmax)
tspan = 1.0:1.0:tmax
obstimes = [5990.0,10.0,0.0,0.0] # S,I.R,C
u0 = [0.05,10.0,0.20, 1.0/180]; # β, c, γ, δ
p
= ODEProblem(sir_ode!,u0,tspan,p);
prob_ode
= solve(prob_ode,
sol_ode Tsit5(),
= 1.0); saveat
```

We can then visualize our expected cases:

```
= Array(sol_ode)[4,:] # Cumulative cases
C = C[2:end] - C[1:(end-1)];
X
Random.seed!(1234)
= rand.(Poisson.(X));
Y
bar(obstimes,Y,legend=false)
plot!(obstimes,X,legend=false)
```

Note the second wave of cases occuring as immunity begins to wane.

## Solve Using Turing

Now we build our model in Turing, using our prior defined ODE. Note that I included a check on the ODE output to try and catch values less than one, otherwise the sampler would reject the solution (because Poisson distributions cannot have a rate parameters less than 0).

```
@model bayes_sir(y) = begin
# Calculate number of timepoints
= length(y)
l ~ Uniform(0.0,.2)
i₀ ~ Uniform(0, .1)
β ~ Uniform(90,365)
immunity
= i₀*6000.0
I =[6000.0-I,I,0.0,0.0]
u0=[β,10.0,0.2, 1.0/immunity]
p= (0.0,float(l))
tspan = ODEProblem(sir_ode!,
prob
u0,
tspan,
p)= solve(prob,
sol Tsit5(),
= 1.0)
saveat = Array(sol)[4,:] # Cumulative cases
sol_C = sol_C[2:end] - sol_C[1:(end-1)]
sol_X = length(y)
l for i in 1:l
~ Poisson(ifelse(sol_X[i] <0 , 1e-6, sol_X[i] ))
y[i] end
end;
= bayes_sir(Y[1:200]);
mod
= sample(mod, NUTS(0.65),10000); ode_nuts
```

`describe(ode_nuts)`

```
2-element Vector{ChainDataFrame}:
Summary Statistics (3 x 8)
Quantiles (3 x 6)
```

Now let’s see if we were able to recover our parameters:

```
plot(ode_nuts)
= DataFrame(ode_nuts);
posterior = mean(posterior[!,:β])
beta_hat = mean(posterior[!, :immunity])
immunity_hat println("Estimate for β: ", beta_hat)
println("Estimate for immunity duration: ", immunity_hat)
println("Estimate for initial proportion infected: ", mean(posterior[!, :i₀]))
```

```
Estimate for β: 0.04984659207646735
Estimate for immunity duration: 178.34561091700505
Estimate for initial proportion infected: 0.0016216537054861632
```

Not too bad! Now we can simulate the dynamics going forward:

```
function predict(y,chain)
# Length of data
= length(y)
l # Length of chain
= length(chain)
m # Choose random
= sample(1:m)
idx = chain[:i₀][idx]
i₀ = chain[:β][idx]
β = chain[:immunity][idx]
immunity = i₀*6000.0
I =[6000.0-I,I,0.0,0.0]
u0=[β,10.0,0.2, 1.0/immunity]
p= (0.0,float(l))
tspan = ODEProblem(sir_ode!,
prob
u0,
tspan,
p)= solve(prob,
sol Tsit5(),
= 1.0)
saveat = Array(sol)
out = [0.0; out[4,2:end] - out[4,1:(end-1)]]
sol_X hcat(sol_ode.t,out',sol_X)
end;
= []
Xp for i in 1:365
= predict(Y,ode_nuts)
pred push!(Xp,pred[2:end,6])
end
scatter(obstimes,Y,legend=false, ylims = [0,600])
plot!(obstimes,Xp,legend=false, color = "grey50", alpha = .2)
```

## Conclusions

This shows a very intuitive approach towards compartmental modeling in a Bayesian framework using Julia.

## Reuse

## Citation

```
@online{dewitt2022,
author = {Michael DeWitt},
title = {Turing and {ODEs}},
date = {2022-08-13},
url = {https://michaeldewittjr.com/programming/2022-08-13-turing-and-odes},
langid = {en}
}
```