Bayesian Inference Examples
Stan
Like in the previous examples, we set up the Lotka-Volterra system and generate data.
f1 = @ode_def begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a b c d
p = [1.5,1.0,3.0,1.0]
u0 = [1.0,1.0]
tspan = (0.0,10.0)
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)Here we now give Stan an array of prior distributions for our parameters. Since the parameters of our differential equation must be positive, we utilize truncated Normal distributions to make sure that is satisfied in the result:
priors = [truncated(Normal(1.5,0.1),0,2),truncated(Normal(1.0,0.1),0,1.5),
truncated(Normal(3.0,0.1),0,4),truncated(Normal(1.0,0.1),0,2)]We then give these to the inference function.
bayesian_result = stan_inference(prob1,t,data,priors;
num_samples=100,num_warmup=500,
vars = (StanODEData(),InverseGamma(4,1)))InverseGamma(4,1) is our starting estimation for the variance hyperparameter of the default Normal distribution. The result is a Mamba.jl chain object. We can pull out the parameter values via:
theta1 = bayesian_result.chain_results[:,["theta.1"],:]
theta2 = bayesian_result.chain_results[:,["theta.2"],:]
theta3 = bayesian_result.chain_results[:,["theta.3"],:]
theta4 = bayesian_result.chain_results[:,["theta.4"],:]From these chains we can get our estimate for the parameters via:
mean(theta1.value[:,:,1])We can get more of a description via:
Mamba.describe(bayesian_result.chain_results)
# Result
Iterations = 1:100
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 100
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
lp__ -6.15472697 1.657551334 0.08287756670 0.18425029767 80.9314979
accept_stat__ 0.90165904 0.125913744 0.00629568721 0.02781181930 20.4968668
stepsize__ 0.68014975 0.112183047 0.00560915237 0.06468790087 3.0075188
treedepth__ 2.68750000 0.524911975 0.02624559875 0.10711170182 24.0159141
n_leapfrog__ 6.77000000 4.121841086 0.20609205428 0.18645821695 100.0000000
divergent__ 0.00000000 0.000000000 0.00000000000 0.00000000000 NaN
energy__ 9.12245750 2.518330231 0.12591651153 0.32894488320 58.6109941
sigma1.1 0.57164997 0.128579363 0.00642896816 0.00444242658 100.0000000
sigma1.2 0.58981422 0.131346442 0.00656732209 0.00397310122 100.0000000
theta1 1.50237077 0.008234095 0.00041170473 0.00025803930 100.0000000
theta2 0.99778276 0.009752574 0.00048762870 0.00009717115 100.0000000
theta3 3.00087782 0.009619775 0.00048098873 0.00020301023 100.0000000
theta4 0.99803569 0.008893244 0.00044466218 0.00040886528 100.0000000
theta.1 1.50237077 0.008234095 0.00041170473 0.00025803930 100.0000000
theta.2 0.99778276 0.009752574 0.00048762870 0.00009717115 100.0000000
theta.3 3.00087782 0.009619775 0.00048098873 0.00020301023 100.0000000
theta.4 0.99803569 0.008893244 0.00044466218 0.00040886528 100.0000000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
lp__ -10.11994750 -7.0569000 -5.8086150 -4.96936500 -3.81514375
accept_stat__ 0.54808912 0.8624483 0.9472840 0.98695850 1.00000000
stepsize__ 0.57975100 0.5813920 0.6440120 0.74276975 0.85282400
treedepth__ 2.00000000 2.0000000 3.0000000 3.00000000 3.00000000
n_leapfrog__ 3.00000000 7.0000000 7.0000000 7.00000000 15.00000000
divergent__ 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000
energy__ 5.54070300 7.2602200 8.7707000 10.74517500 14.91849500
sigma1.1 0.38135240 0.4740865 0.5533195 0.64092575 0.89713635
sigma1.2 0.39674703 0.4982615 0.5613655 0.66973025 0.88361407
theta1 1.48728600 1.4967650 1.5022750 1.50805500 1.51931475
theta2 0.97685115 0.9914630 0.9971435 1.00394250 1.01765575
theta3 2.98354100 2.9937575 3.0001450 3.00819000 3.02065950
theta4 0.97934128 0.9918495 0.9977415 1.00430750 1.01442975
theta.1 1.48728600 1.4967650 1.5022750 1.50805500 1.51931475
theta.2 0.97685115 0.9914630 0.9971435 1.00394250 1.01765575
theta.3 2.98354100 2.9937575 3.0001450 3.00819000 3.02065950
theta.4 0.97934128 0.9918495 0.9977415 1.00430750 1.01442975More extensive information about the distributions is given by the plots:
plot_chain(bayesian_result)Turing
This case we will build off of the Stan example. Note that turing_inference does not require the use of the @ode_def macro like Stan does, but it will still work with macro-defined functions. Thus, using the same setup as before, we simply give the setup to:
bayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500)The result is a MCMCChains.jl chains object. The chain for the first parameter is then given by:
bayesian_result["theta[1]"]Summary statistics can be also be accessed:
using StatsBase
describe(bayesian_result)The chain can be analysed by the trace plots and other plots obtained by:
using StatsPlots
plot(bayesian_result)DynamicHMC
We can use DynamicHMC.jl as the backend for sampling with the dynamic_inference function. It supports any DEProblem, priors can be passed as an array of Distributions.jl distributions, passing initial values is optional and in case where the user has a firm understanding of the domain the parameter values will lie in, transformations can be used to pass an array of constraints for the parameters as an array of Transformations.
bayesian_result_hmc = dynamichmc_inference(prob1, Tsit5(), t, data, [Normal(1.5, 1)], [bridge(ℝ, ℝ⁺, )])A tuple with summary statistics and the chain values is returned. The chain for the ith parameter is given by:
bayesian_result_hmc[1][i]For accessing the various summary statistics:
DynamicHMC.NUTS_statistics(bayesian_result_dynamic[2])Some details about the NUTS sampler can be obtained from:
bayesian_result_dynamic[3]In case of dynamic_inference the trace plots for the ith parameter can be obtained by:
plot(bayesian_result_hmc[1][i])For a better idea of the summary statistics and plotting you can take a look at the benchmarks.