This repository is a set of extension functionality for estimating the parameters of differential equations using Bayesian methods. It allows the choice of using [CmdStan.jl]((https://github.com/StanJulia/CmdStan.jl), Turing.jl, DynamicHMC.jl and ApproxBayes.jl to perform a Bayesian estimation of a differential equation problem specified via the DifferentialEquations.jl interface.

To begin you first need to add this repository using the following command.

Pkg.add("DiffEqBayes")
using DiffEqBayes

stan_inference(prob::DiffEqBase.DEProblem,t,data,priors = nothing;
               alg=:rk45, num_samples=1000, num_warmup=1000, 
               reltol=1e-3, abstol=1e-6, maxiter=Int(1e5),likelihood=Normal,
               vars=(StanODEData(),InverseGamma(3,3)),nchains=1, sample_u0 = false, 
               save_idxs = nothing, diffeq_string = nothing, printsummary = true)

stan_inference uses CmdStan.jl to perform the Bayesian inference. The Stan installation process is required to use this function. The first argument is a DEProblem, t is the array of time and data is the array where the first dimension (columns) corresponds to the array of system values. priors is an array of prior distributions for each parameter, specified via a Distributions.jl type. alg is a choice between :rk45 and :bdf, the two internal integrators of Stan. num_samples is the number of samples to take per chain, and num_warmup is the number of MCMC warmup steps. abstol and reltol are the keyword arguments for the internal integrator. likelihood is the likelihood distribution to use with the arguments from vars, and vars is a tuple of priors for the distributions of the likelihood hyperparameters. The special value StanODEData() in this tuple denotes the position that the ODE solution takes in the likelihood's parameter list. With the diffeq_string kwarg you can pass in a complex ODE specification if the need arises.

turing_inference(prob::DiffEqBase.DEProblem,alg,t,data,priors;
                    likelihood_dist_priors = [InverseGamma(2, 3)], 
                    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
                    num_samples=1000, sampler = Turing.NUTS(0.65),
                    syms = [Turing.@varname(theta[i]) for i in 1:length(priors)],
                    sample_u0 = false, save_idxs = nothing, progress = false, kwargs...)

turing_inference uses Turing.jl to perform its parameter inference. prob can be any DEProblem with a corresponding alg choice. t is the array of time points and data is the set of observations for the differential equation system at time point t[i] (or higher dimensional). priors is an array of prior distributions for each parameter, specified via a Distributions.jl type. num_samples is the number of samples per MCMC chain. The extra kwargs are given to the internal differential equation solver.

dynamichmc_inference(problem::DiffEqBase.DEProblem, algorithm, t, data,parameter_priors, 
                    parameter_transformations=as(Vector, asℝ₊, length(parameter_priors));
                    σ_priors = fill(Normal(0, 5), size(data, 1)),
                    rng = Random.GLOBAL_RNG, num_samples = 1000,
                    AD_gradient_kind = Val(:ForwardDiff),solve_kwargs = (), 
                    mcmc_kwargs = (initialization = (q = zeros(length(parameter_priors) + 2),),), sample_u0 = false)

dynamichmc_inference uses DynamicHMC.jl to perform the bayesian parameter estimation. prob can be any DEProblem, data is the set of observations for our model which is to be used in the Bayesian Inference process. priors represent the choice of prior distributions for the parameters to be determined, passed as an array of Distributions.jl distributions. t is the array of time points. parameter_transformations is an array of Tranformations imposed for constraining the parameter values to specific domains. rng is the random number generator used for MCMC. Defaults to the global one. num_samples is the number of MCMC draws (default: 1000) AD_gradient_kind is passed on to LogDensityProblems.ADgradient, make sure to importthe corresponding library. solve_kwargs is passed on to solve mcmc_kwargs are passed on as keyword arguments to DynamicHMC.mcmc_with_warmup

abc_inference(prob::DEProblem, alg, t, data, priors; ϵ=0.001,
     distancefunction = euclidean, ABCalgorithm = ABCSMC, progress = false,
     num_samples = 500, maxiterations = 10^5, kwargs...)

abc_inference uses ApproxBayes.jl which uses Approximate Bayesian Computation (ABC) to perform its parameter inference. prob can be any DEProblem with a corresponding alg choice. t is the array of time points and data[:,i] is the set of observations for the differential equation system at time point t[i] (or higher dimensional). priors is an array of prior distributions for each parameter, specified via a Distributions.jl type. num_samples is the number of posterior samples. ϵ is the target distance between the data and simulated data. distancefunction is a distance metric specified from the Distances.jl package, the default is euclidean. ABCalgorithm is the ABC algorithm to use, options are ABCSMC or ABCRejection from ApproxBayes.jl, the default is the former which is more efficient. maxiterations is the maximum number of iterations before the algorithm terminates. The extra kwargs are given to the internal differential equation solver.

using ParameterizedFunctions, OrdinaryDiffEq, RecursiveArrayTools, Distributions
f1 = @ode_def LotkaVolterra begin
 dx = a*x - x*y
 dy = -3*y + x*y
end a

p = [1.5]
u0 = [1.0,1.0]
tspan = (0.0,10.0)
prob1 = ODEProblem(f1,u0,tspan,p)

σ = 0.01                         # noise, fixed for now
t = collect(1.:10.)   # observation times
sol = solve(prob1,Tsit5())
priors = [Normal(1.5, 1)]
randomized = VectorOfArray([(sol(t[i]) + σ * randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

using CmdStan #required for using the Stan backend
bayesian_result_stan = stan_inference(prob1,t,data,priors)

bayesian_result_turing = turing_inference(prob1,Tsit5(),t,data,priors)

using DynamicHMC #required for DynamicHMC backend
bayesian_result_hmc = dynamichmc_inference(prob1, Tsit5(), t, data, priors)

bayesian_result_abc = abc_inference(prob1, Tsit5(), t, data, priors)

You don't always have data for all of the variables of the model. In case of certain latent variables you can utilise the save_idxs kwarg to declare the oberved variables and run the inference using any of the backends as shown below.

 sol = solve(prob1,Tsit5(),save_idxs=[1])
 randomized = VectorOfArray([(sol(t[i]) + σ * randn(1)) for i in 1:length(t)])
 data = convert(Array,randomized)

 using CmdStan #required for using the Stan backend
 bayesian_result_stan = stan_inference(prob1,t,data,priors,save_idxs=[1])

 bayesian_result_turing = turing_inference(prob1,Tsit5(),t,data,priors,save_idxs=[1])
 
 using DynamicHMC #required for DynamicHMC backend
 bayesian_result_hmc = dynamichmc_inference(prob1,Tsit5(),t,data,priors,save_idxs = [1])

 bayesian_result_abc = abc_inference(prob1,Tsit5(),t,data,priors,save_idxs=[1])