We can easily solve an ordinary differential equation given the parameters and boundary conditions. For example, given the SIR model as described below

function sir!(du, u, p, t)
    S, I, R = u
    β, γ = p

    N = S + I + R
    a = β * S * I / N
    b = γ * I

    du[1] = -a
    du[2] = a - b
    du[3] = b

    return nothing
end

sir! (generic function with 1 method)

and its parameters, timeframe and boundary condition,

tmax = 50.0
tspan = (0.0,tmax)
obstimes = 1.0:1.0:tmax
u0 = [990.0,10.0,0.0] # S,I.R
p = [0.5,0.25]; # β,γ

we can solve it and get the usual curves.

prob_ode = ODEProblem(sir!,u0,tspan,p);
sol_ode = solve(prob_ode,
            Tsit5(),
            saveat = 1.0);

plot(sol_ode; label=["S(t)" "I(t)" "R(t)"], legend=:outertopright)

Bayesian Estimation

But what if we only has the noisy measurement of the infected number. Can we recover the parameters used? Let's extract the number of infected and add some Poisson noise.

sir_newcases = abs.(diff(sol_ode(0:40; idxs=1).u))
noisy_newcases = rand.(Poisson.(sir_newcases));

bar(noisy_newcases; label="observed")
plot!(sir_newcases; label="model",marker=:circle)

We can construct a bayesian model using Turing.jl

@model function sir_model(t, y, ::Type{T}=Float64) where {T}
    # sample the parameters: uniform priors for `i₀`, `β` and `γ`
    i₀ ~ Uniform(0.0, 1.0)
    β ~ Uniform(0.0, 1.0)
    γ ~ Uniform(0.0, 1.0)

    # simulate the SIR model and save the number of infected
    # individuals at times `t`
    I₀ = 1_000 * i₀
    prob = ODEProblem{true}(sir!, T[1_000 - I₀, I₀, 0], (0.0, last(t)), (; β, γ))
    sol = solve(prob, Tsit5(); saveat=t, save_idxs=1)

    # ensure that the simulation was successful
    if sol.retcode !== :Success
        Turing.@addlogprob! -Inf
    else
        # compute new cases
        cases = map(abs, diff(Array(sol)))

        # noisy observations
        for i in 1:length(y)
            y[i] ~ Poisson(cases[i])
        end
    end

    return (; i₀, β, γ, y)
end;

And use MCMC to estimate the parameters

sir_nuts = sample(
            sir_model(0:40, noisy_newcases),
            NUTS(0.65),
            MCMCThreads(),
            1_000,
            3;
            thinning=10,
            discard_initial=100,
)

┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC /home/catethos/.julia/packages/AbstractMCMC/0eT8o/src/sample.jl:291
┌ Info: Found initial step size
│   ϵ = 0.05
└ @ Turing.Inference /home/catethos/.julia/packages/Turing/uMoX1/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/catethos/.julia/packages/AdvancedHMC/51xgc/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/catethos/.julia/packages/AdvancedHMC/51xgc/src/hamiltonian.jl:47
┌ Info: Found initial step size
│   ϵ = 0.2
└ @ Turing.Inference /home/catethos/.julia/packages/Turing/uMoX1/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/catethos/.julia/packages/AdvancedHMC/51xgc/src/hamiltonian.jl:47
┌ Info: Found initial step size
│   ϵ = 0.025
└ @ Turing.Inference /home/catethos/.julia/packages/Turing/uMoX1/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/catethos/.julia/packages/AdvancedHMC/51xgc/src/hamiltonian.jl:47
Sampling (1 threads): 100%|█████████████████████████████| Time: 0:00:07

Chains MCMC chain (1000×15×3 Array{Float64, 3}):

Iterations        = 101:10:10091
Number of chains  = 3
Samples per chain = 1000
Wall duration     = 16.11 seconds
Compute duration  = 15.92 seconds
parameters        = i₀, β, γ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯

          i₀    0.0092    0.0018     0.0000    0.0000   2936.9917    0.9996    ⋯
           β    0.5000    0.0346     0.0006    0.0007   2915.3128    0.9999    ⋯
           γ    0.2437    0.0297     0.0005    0.0006   2864.2599    1.0004    ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          i₀    0.0063    0.0080    0.0090    0.0103    0.0132
           β    0.4304    0.4770    0.5005    0.5233    0.5657
           γ    0.1843    0.2245    0.2436    0.2639    0.3010

And we can see that mean of the transimission rate is 0.5 and the recovery rate is 0.2437, which is very close to the values we used initially.

plot(sir_nuts)