Epidemic Model (SIR) for Julia and Turing.jl

This repository contains code for an SIR model with a random walk over the $\beta$ parameter. The code simulates an epidemic using a sample random walk and then uses turing.jl to infer the chosen model parameters from the simulated data.

Model definitions

Note that all distributions (except where specified) are parameterised as their default parameterisations in Distributions.jl.

SIR Model

To begin the SIR model chosen is defined as the following:

$$ \begin{align*} \frac{dS}{dt} &= - \frac{\beta(t) S I }{N}, \\ \frac{dI}{dt} &= \frac{\beta(t) S I }{N} - \gamma I ,\\ \frac{dR}{dt} &= \gamma I. \end{align*} $$

In the case of this model $\beta(t)$ is allowed to vary according to a random walk over $\log(\beta)$. This is written as the following (for some sequence of $\beta_k$ where $k$ indexes the sequence):

$$ \begin{equation} \log(\beta_k) \sim \begin{cases} \text{Normal}(\log(0.18), 0.15)\qquad &k = 0 \\ \text{Normal}(\log(\beta_{k-1}), 0.15)\qquad &k > 0 \\ \end{cases} \end{equation} $$

Simulated Data

To begin this model is simulated with data generated by:

$$ \begin{align*} N &= 1,000,000, \\ I_{t_0} &= 100, \\ \gamma &= 0.125 \\ t_0 &= 0, \\ t_{\text{end}} &= 120, \\ \delta t &= 1, \\ \delta_\beta &= 40 \\ \end{align*} $$

Here, $\delta_t$ is the timesteps when the results are saved rather than the timesteps used as part of the ODE solver. Additionally $\delta_\beta$ is the time between the iterations of the random walk. This selected value leads to two $\beta$ values one for the first $40$ days of the model and another for the final $40$ days. The $\beta$ values for the data are drawn from the random walk described in equation~\ref{eq:beta_rw}.

To perform inference observations were drawn from the model, where we draw a sample from a proportion of the number of daily recovered individuals:

$$ \begin{align*} \Delta_R(t) &= R(t) - R(t-1) \\ O(t) & \sim \text{NegativeBinomial}(\Delta_R(t) \times 0.3 + \epsilon, 10) \end{align*} $$

Here $\epsilon$ is some small constant (currently $1e-3$) used for stability. In particular this is added as due to the numerical computational methods it is possible for $\Delta_R(t)$ to be negative, (usually these estimates will be automatically rejected), but this allows the distribution to produce samples in the case the numerical error is negligible. The Negative Binomial distribution is parameterised by:

$$ Y \sim \text{NegativeBinomial}(\mu,\phi), $$

where $\mathbb{E}(Y) = \mu$ and $\text{Var}(Y)=\mu (1-\phi)$.

Priors

The parameters to be estimated in this simple model are $\beta$ and $I_0$ (the initial number of infected individuals. The $\beta_k$ terms use the same priors as their generative distribution. This can be written:

$$ \begin{align*} \log(\beta_{0}) &\sim \text{Normal}(\log(0.18), 0.15), \\ \log(\beta_{k}) &\sim \text{Normal}(\log(\beta_{k-1}), 0.15), \\ \log(I_0) &\sim \text{Normal}(-9.0,0.2). \end{align*} $$

These are written on the log space as the values of the parameters should be strictly positive.