README.md

DiffEqBayesStan.jl

This repository is a set of extension functionality for estimating the parameters of differential equations using Stan-based Bayesian methods as available in StanSample.jl to perform a Bayesian estimation of a differential equation problem specified via the DifferentialEquations.jl interface.

This repository is a simplification of DiffEqBayes.jl. While DiffEqBayes provides and shows how to run the same problem on multiple mcmc implementations available in Julia, this packages only supports Stan.

Version v3.0.0 is a breaking change with v2.x.x in that is assumes cmdstan has been compiled with STAN_THREADS=true. By default 8 CPP threads are used and 4 CPP chains.

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

Pkg.add("DiffEqBayesStan")
using DiffEqBayesStan

Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

Example

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 StanSample                      #required for using the Stan backend
bayesian_result_stan = stan_inference(prob1,t,data,priors)

Using save_idxs to declare observables

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 StanSample #required for using the Stan backend
 bayesian_result_stan = stan_inference(prob1,t,data,priors,save_idxs=[1])