Bayesian Inference of Pendulum Parameters
In this tutorial, we will perform Bayesian parameter inference of the parameters of a pendulum.
Set up simple pendulum problem
using DiffEqBayes, OrdinaryDiffEq, RecursiveArrayTools, Distributions, Plots, StatsPlots,
BenchmarkTools, TransformVariables, StanSample, DynamicHMCLet's define our simple pendulum problem. Here, our pendulum has a drag term ω and a length L.
We get first order equations by defining the first term as the velocity and the second term as the position, getting:
function pendulum(du, u, p, t)
ω, L = p
x, y = u
du[1] = y
du[2] = -ω * y - (9.8 / L) * sin(x)
end
u0 = [1.0, 0.1]
tspan = (0.0, 10.0)
prob1 = ODEProblem(pendulum, u0, tspan, [1.0, 2.5])ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
1.0
0.1Solve the model and plot
To understand the model and generate data, let's solve and visualize the solution with the known parameters:
sol = solve(prob1, Tsit5())
plot(sol)
It's the pendulum, so you know what it looks like. It's periodic, but since we have not made a small angle assumption, it's not exactly sin or cos. Because the true dampening parameter ω is 1, the solution does not decay over time, nor does it increase. The length L determines the period.
Create some dummy data to use for estimation
We now generate some dummy data to use for estimation
t = collect(range(1, stop = 10, length = 10))
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(2)) for i in 1:length(t)])
data = convert(Array, randomized)2×10 Matrix{Float64}:
0.0604183 -0.381395 0.119573 0.0691436 … -0.00671097 0.0211086
-1.20825 0.357891 0.322074 -0.25523 0.0127449 -0.00606225Let's see what our data looks like on top of the real solution
scatter!(data')
This data captures the non-dampening effect and the true period, making it perfect for attempting a Bayesian inference.
Perform Bayesian Estimation
Now let's fit the pendulum to the data. Since we know our model is correct, this should give us back the parameters that we used to generate the data! Define priors on our parameters. In this case, let's assume we don't have much information, but have a prior belief that ω is between 0.1 and 3.0, while the length of the pendulum L is probably around 3.0:
priors = [
truncated(Normal(0.1, 1.0), lower = 0.0),
truncated(Normal(3.0, 1.0), lower = 0.0)
]2-element Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}}:
Truncated(Distributions.Normal{Float64}(μ=0.1, σ=1.0); lower=0.0)
Truncated(Distributions.Normal{Float64}(μ=3.0, σ=1.0); lower=0.0)Finally, let's run the estimation routine from DiffEqBayes.jl with the Turing.jl backend to check if we indeed recover the parameters!
bayesian_result = turing_inference(prob1, Tsit5(), t, data, priors;
syms = [:omega, :L], sample_args = (num_samples = 10_000,))╭─FlexiChain (10000 iterations, 1 chain) ──────────────────────────────────────╮
│ ↓ iter = 1001:11000 │
│ → chain = 1:1 │
│ │
│ Parameters (3) ── AbstractPPL.VarName │
│ Float64 omega, L │
│ Vector{Float64} σ │
│ │
│ 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 │
╰──────────────────────────────────────────────────────────────────────────────╯Notice that while our guesses had the wrong means, the learned parameters converged to the correct means, meaning that it learned good posterior distributions for the parameters. To look at these posterior distributions on the parameters, we can examine the chains:
plot(bayesian_result)
As a diagnostic, we will also check the parameter chains. The chain is the MCMC sampling process. The chain should explore parameter space and converge reasonably well, and we should be taking a lot of samples after it converges (it is these samples that form the posterior distribution!)
plot(bayesian_result, colordim = :parameter)
Notice that after a while these chains converge to a “fuzzy line”, meaning it found the area with the most likelihood and then starts to sample around there, which builds a posterior distribution around the true mean.
DiffEqBayes.jl allows the choice of using Stan.jl, Turing.jl and DynamicHMC.jl for MCMC, you can also use ApproxBayes.jl for Approximate Bayesian computation algorithms. Let's compare the timings across the different MCMC backends. We'll stick with the default arguments and 10,000 samples in each. However, there is a lot of room for micro-optimization specific to each package and algorithm combinations, you might want to do your own experiments for specific problems to get better understanding of the performance.
@btime bayesian_result = turing_inference(prob1, Tsit5(), t, data, priors;
syms = [:omega, :L], sample_args = (num_samples = 10_000,))╭─FlexiChain (10000 iterations, 1 chain) ──────────────────────────────────────╮
│ ↓ iter = 1001:11000 │
│ → chain = 1:1 │
│ │
│ Parameters (3) ── AbstractPPL.VarName │
│ Float64 omega, L │
│ Vector{Float64} σ │
│ │
│ 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 │
╰──────────────────────────────────────────────────────────────────────────────╯@btime bayesian_result = stan_inference(prob1, :rk45, t, data, priors;
sample_kwargs = Dict(:num_samples => 10_000), print_summary = false)10000×4 DataFrame
Row │ sigma1.1 sigma1.2 theta_1 theta_2
│ Float64 Float64 Float64 Float64
───────┼───────────────────────────────────────
1 │ 0.296876 0.208416 1.04929 2.36886
2 │ 0.409671 0.214435 0.976861 2.72706
3 │ 0.275503 0.412957 0.819203 3.14381
4 │ 0.297976 0.474928 0.794019 3.14398
5 │ 0.276293 0.468015 0.604142 3.26871
6 │ 0.186098 0.203846 1.16423 2.30113
7 │ 0.212553 0.309322 0.965825 2.66113
8 │ 0.255146 0.220842 1.25559 2.71028
⋮ │ ⋮ ⋮ ⋮ ⋮
9994 │ 0.30101 0.241782 0.924249 2.87271
9995 │ 0.216602 0.268113 0.984717 2.15207
9996 │ 0.244895 0.19948 1.03245 2.14294
9997 │ 0.188135 0.322937 0.869665 3.16394
9998 │ 0.234297 0.328753 0.925642 3.16433
9999 │ 0.318712 0.311578 0.953452 2.3913
10000 │ 0.355381 0.326101 0.765456 2.3878
9985 rows omitted@btime bayesian_result = dynamichmc_inference(prob1, Tsit5(), t, data, priors;
num_samples = 10_000)(posterior = [(parameters = [0.9976716126147783, 2.498723007826637], σ = [0.01666918662544548, 0.010104515550848727]), (parameters = [0.9829999453412314, 2.488019492735679], σ = [0.015532027341730321, 0.00981410894663699]), (parameters = [0.9866310709885024, 2.4961301949078885], σ = [0.017694759936525632, 0.010361815955728]), (parameters = [0.9923008797774152, 2.4967976892433263], σ = [0.0076089672428565685, 0.009273131934172731]), (parameters = [0.9926895355136842, 2.507321108155256], σ = [0.011680697836312409, 0.011631798217782924]), (parameters = [0.985255011580237, 2.5182540498832906], σ = [0.009773292763196854, 0.008968276951176633]), (parameters = [0.9974083871718615, 2.480207008603498], σ = [0.015392720079414439, 0.010236277279649847]), (parameters = [0.9820482764286853, 2.4898829467512282], σ = [0.00942955514155643, 0.017061802232494788]), (parameters = [1.0033546231356978, 2.483685745394725], σ = [0.01224143769761466, 0.013191514308900313]), (parameters = [0.9808494319178026, 2.5122429030846316], σ = [0.008362065862357649, 0.012892297807207888]) … (parameters = [0.9869346179670394, 2.5050758235541055], σ = [0.013064198180103584, 0.01516825134483162]), (parameters = [0.9941444287748558, 2.4969216482477883], σ = [0.013544290399559206, 0.00966810303878483]), (parameters = [1.0009239932058713, 2.489140325273548], σ = [0.014014106585047256, 0.013930688038595263]), (parameters = [1.014745673138074, 2.5006972940749366], σ = [0.011116910585491484, 0.014989280067267416]), (parameters = [0.9804678781565546, 2.456448744439375], σ = [0.011040317788933097, 0.01730402225698871]), (parameters = [0.9973926207898712, 2.462066162658121], σ = [0.010323963404114551, 0.01607661061841648]), (parameters = [0.9961912176875635, 2.4769591864676763], σ = [0.008551854061286716, 0.01430953217865128]), (parameters = [1.0174267662045853, 2.49089291174018], σ = [0.006399765016338183, 0.019015582113227774]), (parameters = [0.9755413692662005, 2.5023771055954516], σ = [0.01713188896608961, 0.007244574041941086]), (parameters = [0.9917769692191678, 2.5095618990409903], σ = [0.01347503477388078, 0.008202768546998154])], posterior_matrix = [-0.0023311022941881173 -0.017146214439009367 … -0.024762711549390623 -0.008257026391347864; 0.9157798045036474 0.9114870095470774 … 0.917241122348203 0.9201081956946825; -4.094193376124566 -4.164851106765932 … -4.066813700332807 -4.306916581568456; -4.5947728708048565 -4.623934240231047 … -4.927502498320763 -4.803283554012628], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(47.55725047116711, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xa090ef2a)), DynamicHMC.TreeStatisticsNUTS(46.304605342902605, 2, turning at positions -2:1, 0.7649093755396424, 3, DynamicHMC.Directions(0xed2b2439)), DynamicHMC.TreeStatisticsNUTS(46.39824957611053, 2, turning at positions -2:-5, 0.966329695191803, 7, DynamicHMC.Directions(0xd012be32)), DynamicHMC.TreeStatisticsNUTS(45.24439211383659, 3, turning at positions -7:0, 0.7913871990278395, 7, DynamicHMC.Directions(0x0fe000e8)), DynamicHMC.TreeStatisticsNUTS(45.948975538303124, 2, turning at positions -3:0, 0.7206981840737016, 3, DynamicHMC.Directions(0xa9759f68)), DynamicHMC.TreeStatisticsNUTS(45.69840553443062, 2, turning at positions 0:3, 0.8947006414695298, 3, DynamicHMC.Directions(0xe64bebef)), DynamicHMC.TreeStatisticsNUTS(46.54142721804882, 3, turning at positions -4:3, 0.9057976223308417, 7, DynamicHMC.Directions(0x83da230b)), DynamicHMC.TreeStatisticsNUTS(43.408944610724035, 3, turning at positions -6:1, 0.9999999999999999, 7, DynamicHMC.Directions(0x95b57ec9)), DynamicHMC.TreeStatisticsNUTS(46.10487167288595, 3, turning at positions -3:4, 0.998438977549923, 7, DynamicHMC.Directions(0xa30d76ac)), DynamicHMC.TreeStatisticsNUTS(44.46293898791441, 3, turning at positions -6:1, 0.8946721505114522, 7, DynamicHMC.Directions(0x66794ef9)) … DynamicHMC.TreeStatisticsNUTS(47.13341112221245, 3, turning at positions -3:4, 0.992146664396116, 7, DynamicHMC.Directions(0x79f3c954)), DynamicHMC.TreeStatisticsNUTS(45.87372787993138, 3, turning at positions -2:5, 0.8911523017465679, 7, DynamicHMC.Directions(0xa906bbfd)), DynamicHMC.TreeStatisticsNUTS(47.459406788450764, 2, turning at positions -2:-5, 0.9369450598328976, 7, DynamicHMC.Directions(0x3d11fa02)), DynamicHMC.TreeStatisticsNUTS(45.42324529840362, 3, turning at positions -5:2, 0.8232012449357279, 7, DynamicHMC.Directions(0x1ce6a732)), DynamicHMC.TreeStatisticsNUTS(40.3895880879215, 3, turning at positions -5:2, 0.9029312804428035, 7, DynamicHMC.Directions(0x71348082)), DynamicHMC.TreeStatisticsNUTS(41.80263798447655, 2, turning at positions -1:2, 0.9277062763262607, 3, DynamicHMC.Directions(0x794496de)), DynamicHMC.TreeStatisticsNUTS(44.845382594871516, 2, turning at positions -3:-6, 0.9702439686812336, 7, DynamicHMC.Directions(0x5b6230b9)), DynamicHMC.TreeStatisticsNUTS(35.652253955837615, 2, turning at positions -2:1, 0.3433104927439454, 3, DynamicHMC.Directions(0xd4cdfe35)), DynamicHMC.TreeStatisticsNUTS(38.55603431137531, 3, turning at positions -7:0, 0.9999999999999999, 7, DynamicHMC.Directions(0xc7c4b958)), DynamicHMC.TreeStatisticsNUTS(44.42662148594073, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xa4c1f80e))], logdensities = [48.29099824569924, 46.98659042212942, 47.609555309934784, 48.304539868810984, 48.898949824430765, 47.95196580702652, 47.44974319557871, 47.02241771921936, 48.03462629345883, 47.43457595039126 … 47.65894621619614, 49.28123283791366, 47.71227911725524, 45.79867188767801, 42.90460752519547, 45.45024760103899, 47.04879574321352, 38.547654204702496, 44.777201586398135, 48.83613993920133], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.011714459952281788, 0.005269837182639025, 0.22694681760100857, 0.25277587996736844], ϵ = 0.6367544715712073)