Bayesian Inference of ODE
For this tutorial, we will show how to do Bayesian inference to infer the parameters of the Lotka-Volterra equations using each of the three backends:
- Turing.jl
- Stan.jl
- DynamicHMC.jl
Setup
First, let's set up our ODE and the data. For the data, we will simply solve the ODE and take that solution at some known parameters as the dataset. This looks like the following:
using DiffEqBayes, ParameterizedFunctions, OrdinaryDiffEq, RecursiveArrayTools,
Distributions
f1 = @ode_def LotkaVolterra begin
dx = a * x - x * y
dy = -3 * y + x * y
end a
p = [1.5]
u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
prob1 = ODEProblem(f1, u0, tspan, p)
σ = 0.01 # noise, fixed for now
t = collect(1.0:10.0) # observation times
sol = solve(prob1, Tsit5())
priors = [Normal(1.5, 1)]
randomized = VectorOfArray([(sol(t[i]) + σ * randn(2)) for i in 1:length(t)])
data = convert(Array, randomized)2×10 Matrix{Float64}:
2.77422 6.76906 0.973745 1.8835 … 4.36446 3.2481 1.0308
0.256485 2.02699 1.91762 0.312821 0.310398 4.56793 0.89817Inference Methods
Stan
using StanSample #required for using the Stan backend
bayesian_result_stan = stan_inference(prob1, :rk45, t, data, priors)Turing
bayesian_result_turing = turing_inference(prob1, Tsit5(), t, data, priors)╭─FlexiChain (1000 iterations, 1 chain) ───────────────────────────────────────╮
│ ↓ iter = 501:1500 │
│ → chain = 1:1 │
│ │
│ Parameters (2) ── AbstractPPL.VarName │
│ Float64 theta[1] │
│ Vector{Float64} σ (1,) │
│ │
│ Extras (14) │
│ Int64 n_steps, tree_depth │
│ Bool is_accept, numerical_error │
│ Float64 acceptance_rate, log_density, hamiltonian_energy, │
│ hamiltonian_energy_error, max_hamiltonian_energy_error, step_size, │
│ nom_step_size, logprior, loglikelihood, logjoint │
╰──────────────────────────────────────────────────────────────────────────────╯DynamicHMC
We can use DynamicHMC.jl as the backend for sampling with the dynamic_inference function. It is similarly used as follows:
bayesian_result_hmc = dynamichmc_inference(prob1, Tsit5(), t, data, priors)(posterior = [(parameters = [1.500538725384478], σ = [0.010290625715337665, 0.009449236042843118]), (parameters = [1.5006179195620397], σ = [0.01345420199633925, 0.006611136524768659]), (parameters = [1.5005968000625143], σ = [0.012384916809780362, 0.006952427937315721]), (parameters = [1.5004535728861081], σ = [0.008989835250784879, 0.009702824987678407]), (parameters = [1.5005289333848868], σ = [0.013169632555017293, 0.008307492043658717]), (parameters = [1.5003877976593176], σ = [0.016620018942906097, 0.013274547040344446]), (parameters = [1.5006363010336727], σ = [0.008794173120711764, 0.006217069786305236]), (parameters = [1.5006363010336727], σ = [0.008794173120711764, 0.006217069786305236]), (parameters = [1.5004721085649948], σ = [0.006960014209017536, 0.0069594288764821535]), (parameters = [1.5003730827194441], σ = [0.010463739492590484, 0.011335381428046028]) … (parameters = [1.5002026681697274], σ = [0.008371623171470067, 0.008093890742880607]), (parameters = [1.5004520196566504], σ = [0.01187240713607122, 0.007989323548247437]), (parameters = [1.5000550379582782], σ = [0.017333367976434602, 0.009109869837087864]), (parameters = [1.5006323616779167], σ = [0.010015042020824726, 0.010009671069269517]), (parameters = [1.500513622997267], σ = [0.00922867182933649, 0.00800539271029394]), (parameters = [1.5005176674232479], σ = [0.011971108334869128, 0.009056955150961828]), (parameters = [1.500390201612082], σ = [0.009336510032078112, 0.006945673221154631]), (parameters = [1.5006592334937712], σ = [0.009677262967608265, 0.01079647524471664]), (parameters = [1.5007373239576234], σ = [0.013116462306001374, 0.009783964040732963]), (parameters = [1.5005090202436895], σ = [0.010928415387650764, 0.009686686491726665])], posterior_matrix = [0.4058241938854677 0.40587696965624487 … 0.40595653664246184 0.40580439737217816; -4.576521922886018 -4.308463805641329 … -4.333887173136819 -4.516388965561074; -4.661821382771562 -5.0189996996915065 … -4.6270105559174945 -4.637002862876344], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(49.02373890134044, 3, turning at positions -4:3, 0.9311732468523832, 7, DynamicHMC.Directions(0x0d6ffd2b)), DynamicHMC.TreeStatisticsNUTS(49.86396661634861, 3, turning at positions -6:1, 0.9354444483519787, 7, DynamicHMC.Directions(0x9e388a79)), DynamicHMC.TreeStatisticsNUTS(50.34743300709363, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xd8a1498e)), DynamicHMC.TreeStatisticsNUTS(48.83207359077603, 3, turning at positions -5:2, 0.7710103786668887, 7, DynamicHMC.Directions(0xf869876a)), DynamicHMC.TreeStatisticsNUTS(50.13742199218824, 3, turning at positions -3:4, 0.9774363519563051, 7, DynamicHMC.Directions(0x7d485ed4)), DynamicHMC.TreeStatisticsNUTS(48.09354661955847, 2, turning at positions -3:0, 0.7610474157989539, 3, DynamicHMC.Directions(0xb5b570e8)), DynamicHMC.TreeStatisticsNUTS(47.37899225828778, 3, turning at positions -7:0, 0.9946515107905383, 7, DynamicHMC.Directions(0x32af2580)), DynamicHMC.TreeStatisticsNUTS(47.492984800779034, 1, turning at positions -1:0, 0.3204699641611959, 1, DynamicHMC.Directions(0x61e50320)), DynamicHMC.TreeStatisticsNUTS(48.026535233051995, 2, turning at positions -1:2, 0.9415495720162892, 3, DynamicHMC.Directions(0x3e73154a)), DynamicHMC.TreeStatisticsNUTS(46.475801328613784, 2, turning at positions -1:2, 0.695051098858281, 3, DynamicHMC.Directions(0xb210656a)) … DynamicHMC.TreeStatisticsNUTS(48.63568849821959, 2, turning at positions 0:3, 0.9843842422261359, 3, DynamicHMC.Directions(0x5c6ce4c7)), DynamicHMC.TreeStatisticsNUTS(49.46635720164691, 2, turning at positions -3:0, 0.9999999999999999, 3, DynamicHMC.Directions(0x21eb8df4)), DynamicHMC.TreeStatisticsNUTS(48.173276774167874, 2, turning at positions -3:0, 0.7335361514756996, 3, DynamicHMC.Directions(0xa5b16dcc)), DynamicHMC.TreeStatisticsNUTS(48.1165942588604, 2, turning at positions -4:-7, 0.9248650235136304, 7, DynamicHMC.Directions(0xe957c2e0)), DynamicHMC.TreeStatisticsNUTS(50.00370812565485, 3, turning at positions -6:1, 0.9012949812155282, 7, DynamicHMC.Directions(0x06186011)), DynamicHMC.TreeStatisticsNUTS(50.86410960769806, 3, turning at positions -2:5, 0.9929860352334748, 7, DynamicHMC.Directions(0xd24a55f5)), DynamicHMC.TreeStatisticsNUTS(49.26083133620436, 2, turning at positions -3:0, 0.9034100605546723, 3, DynamicHMC.Directions(0x746f522c)), DynamicHMC.TreeStatisticsNUTS(49.63434863177632, 2, turning at positions -1:2, 0.8425697816356861, 3, DynamicHMC.Directions(0xf61c9ce2)), DynamicHMC.TreeStatisticsNUTS(49.259801114373026, 2, turning at positions -1:2, 0.9790975134542009, 3, DynamicHMC.Directions(0x160a4a2e)), DynamicHMC.TreeStatisticsNUTS(49.395207840592676, 2, turning at positions 0:3, 0.9513020500431727, 3, DynamicHMC.Directions(0x376835bb))], logdensities = [51.03689796238601, 50.302356161924735, 50.742370547750696, 50.7825947504678, 50.87022714428036, 48.253526413737205, 49.769406820174126, 49.769406820174126, 48.905787216716504, 50.42558406565631 … 49.59658104502092, 51.1068913640723, 48.28561401703544, 50.607215585509735, 51.05640093019543, 51.03103995897771, 50.68057575002828, 50.14961000351767, 50.19593443708217, 51.008352331340525], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.00015632392408389015, 0.2514390340426794, 0.2543465015426426], ϵ = 0.7700584663466797)More Information
For a better idea of the summary statistics and plotting, you can take a look at the benchmarks.