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.0607877 -0.357317 0.10934 … -0.0164937 -0.0180744 0.0185937
-1.19721 0.342442 0.295411 -0.0491622 0.0172401 -0.0020523Let'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; num_samples = 10_000,
syms = [:omega, :L])Chains MCMC chain (10000×15×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 15.31 seconds
Compute duration = 15.31 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.0600 0.1857 0.0027 5461.3167 4912.8795 1.0000 ⋯
L 2.5235 0.2085 0.0030 4992.9217 4816.8630 1.0000 ⋯
σ[1] 0.1574 0.0363 0.0005 5273.7809 3783.5333 1.0000 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
omega 0.7596 0.9327 1.0369 1.1588 1.4867
L 2.1129 2.3914 2.5215 2.6500 2.9450
σ[1] 0.1004 0.1312 0.1527 0.1783 0.2414
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], 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 = 8.79 seconds
Compute duration = 8.79 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.0608 0.1889 0.0028 5081.0485 4705.9705 1.0001 ⋯
L 2.5256 0.2162 0.0030 5283.0416 4738.3048 0.9999 ⋯
σ[1] 0.1584 0.0369 0.0005 4671.9189 4922.6057 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
omega 0.7651 0.9345 1.0375 1.1572 1.4880
L 2.1167 2.3928 2.5178 2.6513 2.9811
σ[1] 0.1018 0.1327 0.1531 0.1778 0.2477
@btime bayesian_result = stan_inference(prob1, t, data, priors; 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 ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
sigma1.1 0.2606 0.0769 0.0009 9382.8135 6176.2259 1.0001 ⋯
sigma1.2 0.2796 0.0858 0.0011 6950.8465 5706.0265 0.9999 ⋯
theta_1 1.0923 0.2879 0.0035 7477.7311 6252.1427 1.0001 ⋯
theta_2 2.5818 0.3488 0.0042 7143.7481 5965.8497 1.0002 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
sigma1.1 0.1507 0.2063 0.2470 0.3007 0.4463
sigma1.2 0.1578 0.2196 0.2653 0.3210 0.4912
theta_1 0.6498 0.8935 1.0456 1.2415 1.7927
theta_2 1.9433 2.3602 2.5639 2.7737 3.3694
@btime bayesian_result = dynamichmc_inference(prob1, Tsit5(), t, data, priors;
num_samples = 10_000)(posterior = @NamedTuple{parameters::Vector{Float64}, σ::Vector{Float64}}[(parameters = [0.997864876500599, 2.5323783235581256], σ = [0.008584614123674816, 0.013990441359426396]), (parameters = [1.0250519193098198, 2.49399492962621], σ = [0.014623594470234667, 0.013709929978250586]), (parameters = [1.013778418975111, 2.506771715809517], σ = [0.0120612022499371, 0.01087911738150916]), (parameters = [1.0079375332423817, 2.500529444219809], σ = [0.009164644064916705, 0.0212487055441134]), (parameters = [1.0079837375127454, 2.538727142082563], σ = [0.011423156562762092, 0.02049419981657292]), (parameters = [1.0028255210239152, 2.4995089037786453], σ = [0.011815724732708428, 0.016132429839124874]), (parameters = [1.019910206137813, 2.4986137332687965], σ = [0.014446942315026617, 0.012455712818824509]), (parameters = [1.0364121789946699, 2.5051697722682933], σ = [0.01296650771588427, 0.014505158943589452]), (parameters = [1.0097893185598203, 2.4985751755435373], σ = [0.01821710689282607, 0.013311808737292331]), (parameters = [1.0168511082677874, 2.5017960578236322], σ = [0.007086446211780923, 0.013495227388960824]) … (parameters = [1.0316557435162574, 2.501715784100377], σ = [0.014292190141865212, 0.011530128958378854]), (parameters = [1.0308984172653077, 2.4749341418247868], σ = [0.010942962196692661, 0.013445150386602697]), (parameters = [1.0193337959071853, 2.4878607796034835], σ = [0.010202475938889205, 0.013555506245331911]), (parameters = [1.0252780488890993, 2.5176589840635653], σ = [0.009474778174676643, 0.013186552356776739]), (parameters = [1.028259754050354, 2.514175348473979], σ = [0.012319400525613342, 0.009953766768374414]), (parameters = [1.0289796588755848, 2.4788621101143002], σ = [0.011798279827530806, 0.01033128374521914]), (parameters = [1.027183777999011, 2.5107074879326405], σ = [0.013239549133824732, 0.010574943089214512]), (parameters = [1.0244323438359597, 2.516210265011663], σ = [0.013846764685621654, 0.014430257871545282]), (parameters = [1.0137088008107757, 2.4919820720834984], σ = [0.01035962954132547, 0.015126935139985778]), (parameters = [1.0247126545948921, 2.497788582203807], σ = [0.011344817193327096, 0.019183127091970103])], posterior_matrix = [-0.002137406125284091 0.024743264292742276 … 0.01361568523769786 0.024412236308021806; 0.929158910007434 0.9138858142269931 … 0.9130784067113261 0.9154057732953186; -4.757783733380232 -4.225118995081899 … -4.56983880134947 -4.478994274293368; -4.269380942594633 -4.2896348927695644 … -4.191278340109702 -3.9537241836103427], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(42.97985406248991, 2, turning at positions -2:1, 0.9027782153063262, 3, DynamicHMC.Directions(0x1a128fad)), DynamicHMC.TreeStatisticsNUTS(37.62055273745401, 4, turning at positions -10:5, 0.7755235326465495, 15, DynamicHMC.Directions(0x2e979eb5)), DynamicHMC.TreeStatisticsNUTS(44.75316710954127, 4, turning at positions -9:6, 0.9824512677269851, 15, DynamicHMC.Directions(0x8ce96c26)), DynamicHMC.TreeStatisticsNUTS(39.03678828438563, 4, turning at positions -5:10, 0.7372974347895868, 15, DynamicHMC.Directions(0xd2ed159a)), DynamicHMC.TreeStatisticsNUTS(41.28579754375001, 4, turning at positions -4:11, 0.8671920130797819, 15, DynamicHMC.Directions(0x54d7f1cb)), DynamicHMC.TreeStatisticsNUTS(42.28738043331191, 3, turning at positions -11:-14, 0.9429233619218368, 15, DynamicHMC.Directions(0x948a60a1)), DynamicHMC.TreeStatisticsNUTS(43.81164955620392, 4, turning at positions -10:5, 0.998725070321277, 15, DynamicHMC.Directions(0xd05311c5)), DynamicHMC.TreeStatisticsNUTS(44.12740509201883, 3, turning at positions -2:-9, 0.887074417619566, 15, DynamicHMC.Directions(0xe1a5cd56)), DynamicHMC.TreeStatisticsNUTS(43.69617592814777, 3, turning at positions -3:-10, 0.9934851133685837, 15, DynamicHMC.Directions(0x9d9a2935)), DynamicHMC.TreeStatisticsNUTS(43.009457675603045, 4, turning at positions -2:13, 0.9893535015381176, 15, DynamicHMC.Directions(0x6e78e36d)) … DynamicHMC.TreeStatisticsNUTS(44.01947541257551, 3, turning at positions -7:-10, 0.9980049752726005, 11, DynamicHMC.Directions(0x227a8d51)), DynamicHMC.TreeStatisticsNUTS(43.51184695384702, 3, turning at positions 10:13, 0.8713488550890663, 15, DynamicHMC.Directions(0xfd2b2a8d)), DynamicHMC.TreeStatisticsNUTS(43.65989721758397, 2, turning at positions 2:5, 0.9162908276974993, 7, DynamicHMC.Directions(0xce8b560d)), DynamicHMC.TreeStatisticsNUTS(45.177144656974946, 2, turning at positions -3:0, 0.9647948580515159, 3, DynamicHMC.Directions(0xb915a6a4)), DynamicHMC.TreeStatisticsNUTS(42.573044358522175, 4, turning at positions 14:21, 0.8662267488311584, 23, DynamicHMC.Directions(0xfcc27fdd)), DynamicHMC.TreeStatisticsNUTS(44.00312823129502, 4, turning at positions -12:3, 1.0, 15, DynamicHMC.Directions(0xfae1d7e3)), DynamicHMC.TreeStatisticsNUTS(44.19922826347932, 4, turning at positions -6:9, 0.9944893414356132, 15, DynamicHMC.Directions(0x19dffa99)), DynamicHMC.TreeStatisticsNUTS(44.72403202020232, 4, turning at positions -3:12, 1.0, 15, DynamicHMC.Directions(0x873c62ac)), DynamicHMC.TreeStatisticsNUTS(44.69693406924829, 3, turning at positions -3:-10, 0.9987215728610956, 15, DynamicHMC.Directions(0xbe725235)), DynamicHMC.TreeStatisticsNUTS(43.40868838664127, 3, turning at positions 7:14, 0.9584237094046205, 15, DynamicHMC.Directions(0xe509636e))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.025506967767455785, 0.021126478632016926, 0.24968631202776592, 0.2564037710196256], ϵ = 0.21386010318173798)