using CSV, DataFrames, DataFramesMeta, Random
using StatsPlots, LaTeXStrings
using GLM # Linear Model
using ComponentArrays: ComponentArray, ComponentVector
Random.seed!(2)
default(fontfamily="Computer Modern", linewidth=1, label=nothing)Example of classical and Bayesian parameter estimation for ODE models
This notebook demonstrates parameter estimation for a system of ordinary differential equations (ODEs) using both classical optimization and Bayesian inference. The differential equation ecosystem uses SciML packages: DifferentialEquations.jl (OrdinaryDiffEq.jl), DiffEqParamEstim.jl, and Optimization.jl. The Bayesian inference is performed using Turing.jl.
The example focuses on a toy fermentation process model, estimating parameters from data.
I was inspired by this Tutorial: https://turinglang.org/docs/tutorials/bayesian-differential-equations/
Loading Data and Preliminary Analysis
Chargement et affichage des donnees
data = CSV.read("batch.txt", DataFrame; normalizenames=true, delim=',')
data[:, :N] = data[:, :N] / 1000
data[:, :dCO2dt] = data[:, :dCO2dt] / 100
Nbdata = nrow(data) # nombre de données
begin
pB = @df data scatter(:X, :N, xlabel="X", ylabel="N", c=:red)
pS = @df data scatter(:E, :S, xlabel="E", ylabel="S", c=:red)
plot(pB, pS)
end
Linear Regression: Estimating k₁ and k₂
Use linear regression to estimate parameters k₁ and k₂ from the data.
fit_1 = lm(@formula(-N ~ X), data)
c1_est, k1_est = coef(fit_1)
println("Estimated k₁ = ", k1_est)
println("R² for k₁ fit = ", r2(fit_1))
begin
@df data scatter(:X, :N, xlabel="X", ylabel="N", c=:red, label="data")
@df data plot!(:X, -k1_est * :X .- c1_est, c=:red, label="fit")
end
fit_2 = lm(@formula(-S ~ E), data)
c2_est, k2_est = coef(fit_2)
println("Estimated k₂ = ", k2_est)
println("R² for k₂ fit = ", r2(fit_2))
begin
@df data scatter(:E, :S, xlabel="E", ylabel="S", c=:red, label="data")
@df data plot!(:E, -k2_est * :E .- c2_est, c=:red, label="fit")
endEstimated k₁ = 0.050069879583955076
R² for k₁ fit = 0.8821567871700098
Estimated k₂ = 2.08478746168059
R² for k₂ fit = 0.9992787210363997
Parameter Estimation with Differential Equations
Estimate model parameters by fitting an ODE model to the data.
using Optimization, ForwardDiff, OptimizationOptimJL
using OrdinaryDiffEq
using DiffEqParamEstimODE Model Definition
Define the fermentation process as a system of differential equations.
function fermenteur!(du, u, p, t)
# Model parameters.
k₁, k₂, μ₁max, μ₂max, KN, KE, KS = p
# Current state.
X, N, E, S = u
# Calculation of μ₁(S)
μ₁ = μ₁max * N / (KN + N)
# Calculation of μ₂(E,S)
μ₂ = μ₂max * S / (KS + S) * KE / (KE + E)
# Evaluate differential equations.
du[1] = μ₁ * X # dX
du[2] = -k₁ * μ₁ * X # dN
du[3] = μ₂ * X # dE
du[4] = -k₂ * μ₂ * X # dS
return nothing
endfermenteur! (generic function with 1 method)Problem Definition
tspan = (0., 90.)
x0_data_val = Matrix(data[1:1, [:X, :N, :E, :S]]) |> vec
StartList_1 = ComponentVector(k₁=0.01, k₂=2, μ₁max=1.2, μ₂max=1.2, KN=1.6, KE=12, KS=0.03)
prob = ODEProblem(fermenteur!, x0_data_val, tspan, StartList_1)ODEProblem with uType Vector{Float64} and tType Float64. In-place: true Non-trivial mass matrix: false timespan: (0.0, 90.0) u0: 4-element Vector{Float64}: 0.0303083594 0.43247 0.4154047525 188.072160664
Parameter Estimation using classical Optimization
ts_obs = data[:, :time]
data_obs = Matrix(data[:, [:X, :N, :E, :S]]) |> permutedims4×11 Matrix{Float64}:
0.0303084 2.00686 4.07778 … 7.22633 7.52301 7.39643
0.43247 0.18387 0.10636 0.0 0.0 0.0
0.415405 5.82382 14.3733 90.3697 92.5345 92.8707
188.072 182.83 166.707 4.96662 0.135088 0.188063Cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(ts_obs, data_obs),
Optimization.AutoForwardDiff(),
maxiters=10000, verbose=false)
lower = 0.001 * ones(length(StartList_1))
upper = 20 * ones(length(StartList_1)
)7-element Vector{Float64}:
20.0
20.0
20.0
20.0
20.0
20.0
20.0Optimization problem
optprob = Optimization.OptimizationProblem(cost_function, StartList_1; lb=lower, ub=upper
)OptimizationProblem. In-place: true u0: ComponentVector{Float64}(k₁ = 0.01, k₂ = 2.0, μ₁max = 1.2, μ₂max = 1.2, KN = 1.6, KE = 12.0, KS = 0.03)
Solve with different algorithms
result_neldermead = solve(optprob, NelderMead())
result_bfgs = solve(optprob, BFGS())
println("BFGS objective: ", result_bfgs.objective)
println("NelderMead objective: ", result_neldermead.objective)BFGS objective: 49.659601041969196
NelderMead objective: 142.8683779757506Plotting the fitted model
so_bfgs = solve(prob, p=result_bfgs, reltol=1e-6)
so_neldermead = solve(prob, p=result_neldermead, reltol=1e-6)
begin
scatter(data[:, :time], Matrix(data[:, [:X, :N, :E, :S]]))
plot!(so_bfgs, c=permutedims(1:4), label=permutedims(String.([:X, :N, :E, :S])))
plot!(so_neldermead, c=permutedims(1:4), linestyle=:dash, label=:none)
end
Bayesian Inference with Turing.jl
Estimate parameters and uncertainty using Bayesian methods.
using Turing, LinearAlgebra, DistributionsBayesian Model Definition
all_param = [:σX, :σN, :σE, :σS, :k₁, :k₂, :μ₁max, :μ₂max, :KN, :KE, :KS]
all_param_Latex = [L"\sigma_X", L"\sigma_N", L"\sigma_E", L"\sigma_S", L"k_1", L"k_2", L"\mu_{1,max}", L"\mu_{2,max}", L"K_N", L"K_E", L"K_S"]
ode_param = all_param[5:end]
ode_param_Latex = all_param_Latex[5:end]7-element Vector{LaTeXString}:
L"$k_1$"
L"$k_2$"
L"$\mu_{1,max}$"
L"$\mu_{2,max}$"
L"$K_N$"
L"$K_E$"
L"$K_S$"Model definition
@model function fitVin(data, prob, ts)
# Priors for noise and parameters
σ ~ filldist(InverseGamma(3, 2), 4)
k₁ ~ truncated(Normal(prob.p[1], 0.02); lower=0.)
k₂ ~ truncated(Normal(prob.p[2], 0.2); lower=0)
μ₁max ~ truncated(Normal(prob.p[3], 0.2); lower=0)
μ₂max ~ truncated(Normal(prob.p[4], 0.2); lower=0)
KN ~ truncated(Normal(prob.p[5], 0.2); lower=0)
KE ~ truncated(Normal(prob.p[6], 1); lower=0)
KS ~ truncated(Normal(prob.p[7], 0.01); lower=0)
# Simulate ODE with sampled parameters
p = [k₁, k₂, μ₁max, μ₂max, KN, KE, KS]
prob_samp = remake(prob; p=p)
predicted = solve(prob_samp, Tsit5(); saveat=ts, abstol=1e-6)
# Likelihood: compare model to data
for i in eachindex(predicted)
data[:, i] ~ MvNormal(predicted[i], Diagonal(σ) .^ 2)
end
return nothing
end
probT = ODEProblem(fermenteur!, x0_data_val, tspan, [0.06, 2., 0.2, 0.94, 0.01, 19.7, 0.76])
obs = Matrix(data[:, [:X, :N, :E, :S]]) |> permutedims
model = fitVin(obs, probT, data[:, :time])DynamicPPL.Model{typeof(fitVin), (:data, :prob, :ts), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fermenteur!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}(fitVin, (data = [0.0303083594 2.006857401 … 7.5230052962 7.3964325141; 0.43247 0.18387 … 0.0 0.0; 0.4154047525 5.8238227786 … 92.5345073692 92.870674345; 188.072160664 182.8301080392 … 0.135087753 0.1880626952], prob = ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fermenteur!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}(ODEFunction{true, SciMLBase.AutoSpecialize, typeof(fermenteur!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}(fermenteur!, UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing, nothing, nothing), [0.0303083594, 0.43247, 0.4154047525, 188.072160664], (0.0, 90.0), [0.06, 2.0, 0.2, 0.94, 0.01, 19.7, 0.76], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}(), SciMLBase.StandardODEProblem()), ts = [0.0, 16.4, 20.9, 24.0, 26.1, 40.5, 45.7, 49.4, 65.4, 74.9, 88.0]), NamedTuple(), DynamicPPL.DefaultContext())# using Turing.DynamicPPL
# DynamicPPL.DebugUtils.model_warntype(model)Sampling from the Prior and Posterior
@time "Sampling from the Prior" sample_prior = sample(model, Prior(), 1000; progress=false)
@time "Sampling from the Posterior" sample_posterior = sample(model, NUTS(), 1000; progress=false)┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 Sampling from the Prior: 10.193871 seconds (7.08 M allocations: 289.146 MiB, 1.42% gc time, 41.99% compilation time) ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Info: Found initial step size └ ϵ = 0.00625 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 Sampling from the Posterior: 361.616857 seconds (114.88 M allocations: 21.808 GiB, 4.34% gc time, 3.89% compilation time)
Chains MCMC chain (1000×25×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 357.85 seconds
Compute duration = 357.85 seconds
parameters = σ[1], σ[2], σ[3], σ[4], k₁, k₂, μ₁max, μ₂max, KN, KE, KS
internals = 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, lp, logprior, loglikelihood
Use `describe(chains)` for summary statistics and quantiles.Results analysis
Posterior Analysis
plot(sample_posterior)
begin
plt = [stephist(sample_posterior[key], label=ode_param_Latex[i]) for (i, key) in enumerate(ode_param)]
[stephist!(plt[i], sample_prior[key], label="Prior") for (i, key) in enumerate(ode_param)]
[vline!(plt[i], [result_bfgs[i]]) for (i, k) in enumerate(ode_param)]
plt_σ = [stephist(sample_posterior[key], label=string(all_param[i])) for (i, key) in enumerate([Symbol("σ[1]"), Symbol("σ[2]"), Symbol("σ[3]"), Symbol("σ[4]")])]
[stephist!(plt_σ[i], sample_prior[key], label="Prior") for (i, key) in enumerate([Symbol("σ[1]"), Symbol("σ[2]"), Symbol("σ[3]"), Symbol("σ[4]")])]
plot(plt..., plt_σ...)
end
Parameter Correlations
heatmap(cor(Array(sample_posterior)), yflip=true, xticks=(1:11, all_param_Latex), yticks=(1:11, all_param_Latex), clims=(-1, 1), color=:balance)
Posterior Predictive Check
Simulate data from the posterior and compare to observed data.
Ys = map(eachrow(Array(sample_posterior[ode_param]))) do ps
sol_p = Array(solve(prob; p=ps, saveat=0.1, reltol=1e-6))
end |> stack
ts = 0:0.1:90
styles = [:auto :dot :auto :dash]
vars = String.([:X, :N, :E, :S])
begin
plt = plot(title="Data retrodiction")
for (i, label) in enumerate(vars)
errorline!(plt, ts, Ys[i, :, :], label=label, c=i, groupcolor=i, errortype=:percentile, percentiles=[0, 100], fillalpha=0.8)
end
plot!(so_bfgs, c=:black, s=:dash, label=:none, linewidth=2)
[scatter!(data[:, :time], (data[:, Symbol(label)]), c=:red, label=:none) for label in vars]
plt
end
Handling Missing Observations
You can fit the model even if some variables are not observed (e.g., missing N).
@model function fitVin_Missing(data, prob, ts)
# Prior distributions.
σ ~ filldist(InverseGamma(3, 2), 3)
k₁ ~ truncated(Normal(prob.p[1], 0.02); lower=0.)
k₂ ~ truncated(Normal(prob.p[2], 0.2); lower=0)
μ₁max ~ truncated(Normal(prob.p[3], 0.2); lower=0)
μ₂max ~ truncated(Normal(prob.p[4], 0.2); lower=0)
KN ~ truncated(Normal(prob.p[5], 0.2); lower=0)
KE ~ truncated(Normal(prob.p[6], 1); lower=0)
KS ~ truncated(Normal(prob.p[7], 0.01); lower=0)
# Simulate ODE model.
p = [k₁, k₂, μ₁max, μ₂max, KN, KE, KS]
prob_samp = remake(prob; p=p)
predicted = solve(prob_samp, Tsit5(); saveat=ts, abstol=1e-6, save_idxs=1:3)
# Observations (N is missing, so only X, E, S are used).
for i in eachindex(predicted)
data[:, i] ~ MvNormal(predicted[i], Diagonal(σ) .^ 2)
end
return nothing
endfitVin_Missing (generic function with 2 methods)Sampling with Missing Observations
obs_missing = Matrix(data[:, [:X, :E, :S]]) |> permutedims
model_missing = fitVin_Missing(obs_missing, probT, data[:, :time])
@time "Sampling Posterior with N missing" sample_posterior_missing = sample(model_missing, NUTS(), 1000; progress=false)
@time "Sampling Prior with N missing" sample_prior_missing = sample(model_missing, Prior(), 1000; progress=false)┌ Info: Found initial step size └ ϵ = 0.003125 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 Sampling Posterior with N missing: 195.582806 seconds (45.91 M allocations: 4.800 GiB, 0.82% gc time, 22.15% compilation time) ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 ┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\metivier\.julia\packages\SciMLBase\9sEkh\src\integrator_interface.jl:623 Sampling Prior with N missing: 1.453757 seconds (6.04 M allocations: 234.336 MiB, 1.10% gc time, 22.92% compilation time)
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 1.17 seconds
Compute duration = 1.17 seconds
parameters = σ[1], σ[2], σ[3], k₁, k₂, μ₁max, μ₂max, KN, KE, KS
internals = lp, logprior, loglikelihood
Use `describe(chains)` for summary statistics and quantiles.Posterior Analysis with Missing Observations
begin
plt_missing = [stephist(sample_posterior_missing[key], label=ode_param_Latex[i]) for (i, key) in enumerate(ode_param)]
[stephist!(plt_missing[i], sample_prior_missing[key], label="Prior") for (i, key) in enumerate(ode_param)]
[vline!(plt_missing[i], [result_bfgs[i]]) for (i, k) in enumerate(ode_param)]
plt_σ = [stephist(sample_posterior_missing[key], label=string(all_param[i])) for (i, key) in enumerate([Symbol("σ[1]"), Symbol("σ[2]"), Symbol("σ[3]")])]
[stephist!(plt_σ[i], sample_prior_missing[key], label="Prior") for (i, key) in enumerate([Symbol("σ[1]"), Symbol("σ[2]"), Symbol("σ[3]")])]
plot(plt_missing..., plt_σ...)
end
Correlation plot with missing observations
heatmap(cor(Array(sample_posterior_missing)), yflip=true, xticks=(1:10, all_param_Latex[[1:3; 5:11]]), yticks=(1:10, all_param_Latex[[1:3; 5:11]]), clims=(-1, 1), color=:balance)
Retrodiction plot for missing observations
Ys_missing = map(eachrow(Array(sample_posterior_missing[ode_param]))) do ps
Array(solve(probT; p=ps, saveat=0.1, reltol=1e-6))
end |> stack
begin
plt = plot(title="Data retrodiction (missing observations of N)", legend=false)
for (i, label) in enumerate(vars)
errorline!(plt, ts, Ys_missing[i, :, :], label=label, c=i, groupcolor=i, errortype=:percentile, percentiles=[0, 100], fillalpha=0.8)
end
[scatter!(data[:, :time], data[:, Symbol(label)], c=:red, label=:none) for label in vars]
plt
end
Julia, Computer and Packages settings
Code
using InteractiveUtils
InteractiveUtils.versioninfo()
import Pkg;
Pkg.status();Julia Version 1.11.4 Commit 8561cc3d68 (2025-03-10 11:36 UTC) Build Info: Official https://julialang.org/ release Platform Info: OS: Windows (x86_64-w64-mingw32) CPU: 28 × 13th Gen Intel(R) Core(TM) i7-13850HX WORD_SIZE: 64 LLVM: libLLVM-16.0.6 (ORCJIT, alderlake) Threads: 28 default, 0 interactive, 14 GC (on 28 virtual cores) Environment: JULIA_LOAD_PATH = @;@stdlib JULIA_PROJECT = @. Status `C:\Users\metivier\Dropbox\PC (2)\Documents\Simulations\julia_mwe\Bayesian_EqDiff\Project.toml` [336ed68f] CSV v0.10.15 [b0b7db55] ComponentArrays v0.15.29 [a93c6f00] DataFrames v1.7.1 [1313f7d8] DataFramesMeta v0.15.4 [1130ab10] DiffEqParamEstim v2.2.0 [31c24e10] Distributions v0.25.120 [f6369f11] ForwardDiff v1.1.0 [38e38edf] GLM v1.9.0 [b964fa9f] LaTeXStrings v1.4.0 [7f7a1694] Optimization v4.6.0 [36348300] OptimizationOptimJL v0.4.3 [1dea7af3] OrdinaryDiffEq v6.102.0 [f3b207a7] StatsPlots v0.15.7 [fce5fe82] Turing v0.40.2 [b77e0a4c] InteractiveUtils v1.11.0
This page was generated using Literate.jl.