SIR model with uncertainty

ode_uncertainty.jmd

# Ordinary differential equation model with uncertainty quantification
Simon Frost (@sdwfrost), 2021-12-20

## Introduction

The classical ODE version of the SIR model is:

- Deterministic
- Continuous in time
- Continuous in state

In this notebook, we accommodate uncertainty in parameter values on the numerical solution of the ODE system in different ways.

## Libraries

```julia
using DifferentialEquations
using OrdinaryDiffEq
using DiffEqUncertainty
using Distributions
using Plots
```

## Transitions

The following function provides the derivatives of the model, which it changes in-place. State variables and parameters are unpacked from `u` and `p`; this incurs a slight performance hit, but makes the equations much easier to read.

A variable is included for the cumulative number of infections, $C$.

```julia
function sir_ode!(du,u,p,t)
    (S,I,R,C) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end;
```

## Time domain

We set the timespan for simulations, `tspan`, initial conditions, `u0`, and parameter values, `p` (which are unpacked above as `[β,γ]`).

```julia
δt = 1.0
tmax = 40.0
tspan = (0.0,tmax)
obstimes = 1.0:1.0:tmax;
```

## Initial conditions

```julia
u0 = [990.0,10.0,0.0,0.0]; # S,I.R,Y
```

## Parameter values with uncertainty


```julia
p_fixed = [0.05,10.0,0.25]; # β,c,γ
p_uncertain = [Uniform(0.01,0.09),10.0,0.25];
```

## Running the model

```julia
prob_ode = ODEProblem(sir_ode!,u0,tspan,p_fixed)
sol_ode = solve(prob_ode,Tsit5(),saveat=δt);
```

## Using DiffEqUncertainty.jl

We first define a function that returns the observables of interest - here, the cumulative number of infections over time, `C`.

```julia
g(sol) = hcat(sol(obstimes).u...)[4,:]
```

### Monte Carlo

```julia
C_mean_mc = expectation(g, prob_ode, u0, p_uncertain, MonteCarlo(), Tsit5(); nout=40, trajectories = 1000)
```

```julia
C_var_mc = centralmoment(2, g, prob_ode, u0, p_uncertain, MonteCarlo(), Tsit5(); nout=40, trajectories = 1000)
```

### Koopman

```julia
C_mean_koopman = expectation(g, prob_ode, u0, p_uncertain, Koopman(), Tsit5(); nout=40)
```

```julia
C_var_koopman = centralmoment(2, g, prob_ode, u0, p_uncertain, Koopman(), Tsit5(); nout=40)
```

## Plotting

```julia
plot(obstimes,C_mean_mc)
plot!(obstimes,C_mean_koopman)
```