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, CmdStan, 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
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.0642063 -0.383084 0.12968 … -0.0276154 -0.00289514 0.0126156
-1.21817 0.35675 0.324259 -0.023583 0.0155866 0.00109738Let'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,))Chains MCMC chain (10000×15×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 18.68 seconds
Compute duration = 18.68 seconds
parameters = omega, L, σ[1]
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 mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
omega 1.0260 0.1733 0.0026 5211.9476 4585.3839 1.0000 ⋯
L 2.5137 0.1935 0.0025 6033.3051 4837.9827 1.0000 ⋯
σ[1] 0.1584 0.0378 0.0006 4618.5587 4599.0182 1.0002 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
omega 0.7510 0.9121 1.0037 1.1171 1.4348
L 2.1421 2.3916 2.5071 2.6316 2.9209
σ[1] 0.1008 0.1316 0.1528 0.1791 0.2490
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,))Chains MCMC chain (10000×15×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 17.57 seconds
Compute duration = 17.57 seconds
parameters = omega, L, σ[1]
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 mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
omega 3.1096 0.4396 0.0161 835.6966 375.9912 1.0006 ⋯
L 0.1309 0.0177 0.0005 817.4396 717.5835 1.0022 ⋯
σ[1] 0.2566 0.0564 0.0012 2328.1224 2408.0799 1.0005 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
omega 2.1438 2.8276 3.1440 3.4162 3.9003
L 0.1096 0.1177 0.1268 0.1396 0.1751
σ[1] 0.1724 0.2161 0.2489 0.2870 0.3896
@btime bayesian_result = stan_inference(prob1, :rk45, t, data, priors;
sample_kwargs = Dict(:num_samples => 10_000), print_summary = false)Chains MCMC chain (10000×4×1 Array{Float64, 3}):
Iterations = 1:1:10000
Number of chains = 1
Samples per chain = 10000
parameters = sigma1.1, sigma1.2, theta_1, theta_2
internals =
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
sigma1.1 0.2804 0.0811 0.0062 146.3898 967.0798 1.0668 ⋯
sigma1.2 0.3271 0.1188 0.0146 59.1124 323.6778 1.1854 ⋯
theta_1 1.7420 0.9334 0.1938 27.9240 186.1701 1.5547 ⋯
theta_2 1.5990 1.2200 0.2656 26.9401 39.1245 1.5694 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
sigma1.1 0.1582 0.2225 0.2700 0.3201 0.4779
sigma1.2 0.1670 0.2443 0.3061 0.3785 0.6252
theta_1 0.6822 0.9659 1.2985 2.6217 3.5834
theta_2 0.1039 0.1431 2.2842 2.6153 3.1751
@btime bayesian_result = dynamichmc_inference(prob1, Tsit5(), t, data, priors;
num_samples = 10_000)(posterior = @NamedTuple{parameters::Vector{Float64}, σ::Vector{Float64}}[(parameters = [0.9866285111019449, 2.4924235028877066], σ = [0.008735291466267031, 0.01084189998168222]), (parameters = [1.0044235813527567, 2.4757312642569107], σ = [0.007452931658981959, 0.012112872331162403]), (parameters = [1.0038921921773178, 2.4790444682692536], σ = [0.0068293466502139764, 0.012506345718633785]), (parameters = [1.0103985449459825, 2.485501970884515], σ = [0.006750115534980478, 0.013642604753319066]), (parameters = [0.9602720448411853, 2.5075496478513126], σ = [0.008860788234809396, 0.014520225285060154]), (parameters = [1.0029615424527945, 2.4771369967051253], σ = [0.012378329872284074, 0.01971636903464804]), (parameters = [0.9934010658426567, 2.4924125187265314], σ = [0.005439018876970695, 0.005255565969649078]), (parameters = [0.9912051421668614, 2.4851613807595077], σ = [0.005558696536570486, 0.005205211267749636]), (parameters = [0.995797501684715, 2.4910901035192086], σ = [0.005192520800588318, 0.006165280422727907]), (parameters = [0.990643322111553, 2.4843815196227057], σ = [0.005300039694939202, 0.006268443738684124]) … (parameters = [0.9651219520693896, 2.5140932950257064], σ = [0.007762739170161695, 0.010868077333232902]), (parameters = [0.9950145449545306, 2.482576180359437], σ = [0.009229087443561766, 0.009614490771068825]), (parameters = [0.992994596650999, 2.499047325883687], σ = [0.008511736582203763, 0.01024622741783435]), (parameters = [0.9824758393129263, 2.493080656771596], σ = [0.010540832367947998, 0.008335326569488823]), (parameters = [0.9909923800859891, 2.4894183959808935], σ = [0.0067399553255883425, 0.010402928279938]), (parameters = [0.9911582193810042, 2.485480780219747], σ = [0.011647915790794, 0.008343136476888571]), (parameters = [0.9898507612926479, 2.496254506226338], σ = [0.009691838051532576, 0.010467302181122825]), (parameters = [0.9786014945019131, 2.498356399658978], σ = [0.0057258498407105, 0.00868737406204171]), (parameters = [0.9756323641634377, 2.500116718300259], σ = [0.0056698028105690995, 0.009173555863490114]), (parameters = [0.982194857873756, 2.4923499656413868], σ = [0.011819813803654313, 0.007750112419326874])], posterior_matrix = [-0.013461692260351328 0.004413826075027321 … -0.024669439591343927 -0.017965560705545437; 0.9132555314652102 0.9065358126866756 … 0.9163374181044397 0.9132260267156806; -4.740383968320631 -4.899147812723178 … -5.172600939524641 -4.437978019780179; -4.524337023266581 -4.413486562822846 … -4.691430296504684 -4.860047930002492], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(50.815114146236475, 3, turning at positions -4:3, 0.9969390771881855, 7, DynamicHMC.Directions(0xc178770b)), DynamicHMC.TreeStatisticsNUTS(47.85741484272149, 3, turning at positions -7:0, 0.9220183262897909, 7, DynamicHMC.Directions(0xf053eb68)), DynamicHMC.TreeStatisticsNUTS(47.65067443281233, 3, turning at positions -6:1, 0.9700706226519743, 7, DynamicHMC.Directions(0xb451a309)), DynamicHMC.TreeStatisticsNUTS(46.60138829775463, 3, turning at positions -2:5, 0.9819498940167415, 7, DynamicHMC.Directions(0xad188c45)), DynamicHMC.TreeStatisticsNUTS(46.29218978733734, 3, turning at positions -8:-15, 0.9999433705813103, 15, DynamicHMC.Directions(0xd2a29970)), DynamicHMC.TreeStatisticsNUTS(45.82320398267571, 4, turning at positions -13:2, 0.985136485865738, 15, DynamicHMC.Directions(0x44322cd2)), DynamicHMC.TreeStatisticsNUTS(43.87156249584445, 4, turning at positions -2:-17, 0.9960002189961126, 31, DynamicHMC.Directions(0x7bf134ce)), DynamicHMC.TreeStatisticsNUTS(45.53385646016987, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x88ccbc12)), DynamicHMC.TreeStatisticsNUTS(42.52055048732671, 4, turning at positions -2:13, 0.7577636348104974, 15, DynamicHMC.Directions(0x092a6f4d)), DynamicHMC.TreeStatisticsNUTS(46.302224611716014, 2, turning at positions -3:0, 0.9604701440613873, 3, DynamicHMC.Directions(0x1f480c28)) … DynamicHMC.TreeStatisticsNUTS(46.15369211512939, 3, turning at positions -8:-11, 0.9861088657158601, 11, DynamicHMC.Directions(0x69538ef0)), DynamicHMC.TreeStatisticsNUTS(44.44554555514706, 3, turning at positions -6:-13, 0.9830053057544706, 15, DynamicHMC.Directions(0x599c4d92)), DynamicHMC.TreeStatisticsNUTS(50.76865071627392, 3, turning at positions 7:10, 0.9611201829341808, 15, DynamicHMC.Directions(0x05cb166a)), DynamicHMC.TreeStatisticsNUTS(51.15641491321499, 4, turning at positions 4:11, 0.9982105654249295, 23, DynamicHMC.Directions(0x9ee71dd3)), DynamicHMC.TreeStatisticsNUTS(50.87599887133985, 5, turning at positions -23:8, 0.9948753340993552, 31, DynamicHMC.Directions(0x5e9fb208)), DynamicHMC.TreeStatisticsNUTS(50.432533832219555, 5, turning at positions -1:30, 0.9882896031718997, 31, DynamicHMC.Directions(0xc0c70c7e)), DynamicHMC.TreeStatisticsNUTS(48.44553675792637, 3, turning at positions 7:14, 0.9880300381032101, 15, DynamicHMC.Directions(0x5a2503ce)), DynamicHMC.TreeStatisticsNUTS(47.903237789455694, 4, turning at positions -13:2, 0.962878621419239, 15, DynamicHMC.Directions(0xb3eda342)), DynamicHMC.TreeStatisticsNUTS(49.9615443889668, 3, turning at positions -3:4, 0.9886184420485995, 7, DynamicHMC.Directions(0x1bf71e3c)), DynamicHMC.TreeStatisticsNUTS(49.58456337337092, 4, turning at positions -1:14, 0.9858349783915741, 15, DynamicHMC.Directions(0x65e78b2e))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.022995310351715852, 0.021298850116413384, 0.26052758459039455, 0.26569785128090273], ϵ = 0.14545781093573146)