Merge pull request #40 from marcjwilliams1/abcinference

Added ABC inference using ApproxBayes.jl

README.md

@@ -4,13 +4,13 @@
 [![Coverage Status](https://coveralls.io/repos/JuliaDiffEq/DiffEqBayes.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/JuliaDiffEq/DiffEqBayes.jl?branch=master)
 [![codecov.io](http://codecov.io/github/JuliaDiffEq/DiffEqBayes.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaDiffEq/DiffEqBayes.jl?branch=master)
 
-This repository is a set of extension functionality for estimating the parameters of differential equations using Bayesian methods. It allows the choice of using [Stan.jl](https://github.com/goedman/Stan.jl), [Turing.jl](https://github.com/yebai/Turing.jl) and [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl) to perform a Bayesian estimation of a differential equation problem specified via the [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl) interface. 
+This repository is a set of extension functionality for estimating the parameters of differential equations using Bayesian methods. It allows the choice of using [Stan.jl](https://github.com/goedman/Stan.jl), [Turing.jl](https://github.com/yebai/Turing.jl), [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl) and [ApproxBayes.jl](https://github.com/marcjwilliams1/ApproxBayes.jl) to perform a Bayesian estimation of a differential equation problem specified via the [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl) interface.
 
 To begin you first need to add this repository using the following command.
 ```julia
 Pkg.add("DiffEqBayes")
 using DiffEqBayes
-``` 
+```
 
 ### stan_inference
 
@@ -63,15 +63,38 @@ dynamichmc_inference(prob::DEProblem,data,priors,t,transformations;
                       σ = 0.01,ϵ=0.001,initial=Float64[])
 ```
 
-`dynamichmc_inference` uses [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl) to 
- perform the bayesian parameter estimation. `prob` can be any `DEProblem`, `data` is the set 
- of observations for our model whihc 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]
- (https://juliastats.github.io/Distributions.jl/latest/) distributions. `t` is the array of time points. `transformations`
- is an array of [Tranformations](https://github.com/tpapp/ContinuousTransformations.jl) imposed for constraining the 
+`dynamichmc_inference` uses [DynamicHMC.jl](https://github.com/tpapp/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](https://juliastats.github.io/Distributions.jl/latest/) distributions. `t` is the array of time points. `transformations`
+ is an array of [Tranformations](https://github.com/tpapp/ContinuousTransformations.jl) imposed for constraining the
  parameter values to specific domains. `initial` values for the parameters can be passed, if not passed the means of the
  `priors` are used. `ϵ` can be used as a kwarg to pass the initial step size for the NUTS algorithm.
 
+### abc_inference
+
+```julia
+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](https://github.com/marcjwilliams1/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](https://juliastats.github.io/Distributions.jl/latest/)
+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](https://github.com/JuliaStats/Distances.jl)
+package, the default is `euclidean`. `ABCalgorithm` is the ABC algorithm to use, options are `ABCSMC` or `ABCRejection` from
+[ApproxBayes.jl](https://github.com/marcjwilliams1/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.
+
+
  ## Example
 
  ```julia
@@ -94,9 +117,12 @@ dynamichmc_inference(prob::DEProblem,data,priors,t,transformations;
 
  bayesian_result_stan = stan_inference(prob1,t,data,priors;num_samples=300,
                                  num_warmup=500,likelihood=Normal,
-                                 vars =(StanODEData(),InverseGamma(3,2))) 
+                                 vars =(StanODEData(),InverseGamma(3,2)))
 
  bayesian_result_turing = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500)
 
  bayesian_result_hmc = dynamichmc_inference(prob1, data, [Normal(1.5, 1)], t, [bridge(ℝ, ℝ⁺, )])
- ```
+
+ bayesian_result_abc = abc_inference(prob1, Tsit5(), t, data, [Normal(1.5, 1)];
+                                 num_samples=500)
+ ```

REQUIRE

@@ -13,3 +13,5 @@ Parameters
 DiffWrappers 0.0.1 0.1.0
 ContinuousTransformations 0.0.1 0.1.0
 DynamicHMC 0.0.1 0.1.0
+Distances
+ApproxBayes

src/DiffEqBayes.jl

@@ -4,13 +4,15 @@ using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
 using DynamicHMC, DiffWrappers, ContinuousTransformations
 using Parameters, Distributions, Optim
 using Compat
+using Distances, ApproxBayes
 
 include("stan_inference.jl")
 include("turing_inference.jl")
 include("stan_string.jl")
 include("utils.jl")
 include("dynamichmc_inference.jl")
+include("abc_inference.jl")
 
-export StanModel, stan_inference, turing_inference, stan_string, StanODEData, plot_chain, dynamichmc_inference, LotkaVolterraPosterior
+export StanModel, stan_inference, turing_inference, stan_string, StanODEData, plot_chain, dynamichmc_inference, LotkaVolterraPosterior, abc_inference
 
 end # module

src/abc_inference.jl

@@ -0,0 +1,23 @@
+function createabcfunction(prob, t, distancefunction, alg, kwargs...)
+    function simfunc(params, constants, targetdata)
+        sol = solve(problem_new_parameters(prob, params), alg, saveat = t, save_start = false, kwargs...)
+        data = convert(Array, sol)
+        distancefunction(targetdata, data), data
+    end
+end
+
+function abc_inference(prob::DEProblem, alg, t, data, priors; ϵ=0.001,
+     distancefunction = euclidean, ABCalgorithm = ABCSMC, progress = false,
+     num_samples = 500, maxiterations = 10^5, kwargs...)
+
+    abcsetup = ABCalgorithm(createabcfunction(prob, t, distancefunction, alg, kwargs...),
+     length(priors),
+     ϵ,
+     ApproxBayes.Prior(priors);
+     nparticles = num_samples,
+     maxiterations = maxiterations
+    )
+
+    abcresult = runabc(abcsetup, data, progress = progress)
+    return abcresult
+end

test/REQUIRE

@@ -1,3 +1,4 @@
 OrdinaryDiffEq
 RecursiveArrayTools
 ParameterizedFunctions
+StatsBase

test/abc.jl

@@ -0,0 +1,49 @@
+using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, Distances, StatsBase, RecursiveArrayTools
+using Base.Test
+
+println("One parameter case")
+f1 = @ode_def_nohes LotkaVolterraTest1 begin
+  dx = a*x - x*y
+  dy = -3y + x*y
+end a
+u0 = [1.0,1.0]
+tspan = (0.0,10.0)
+prob1 = ODEProblem(f1,u0,tspan,[1.5])
+sol = solve(prob1,Tsit5())
+t = collect(linspace(1,10,10))
+randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
+data = convert(Array,randomized)
+priors = [Normal(1.5,0.01)]
+
+bayesian_result = abc_inference(prob1,Tsit5(),t,data,priors;
+                                num_samples=500,ϵ = 0.001)
+
+@show mean(bayesian_result.parameters, weights(bayesian_result.weights))
+@test mean(bayesian_result.parameters, weights(bayesian_result.weights)) ≈ 1.5 atol=0.1
+
+println("Four parameter case")
+f1 = @ode_def_nohes LotkaVolterraTest4 begin
+  dx = a*x - b*x*y
+  dy = -c*y + d*x*y
+end a b c d
+u0 = [1.0,1.0]
+tspan = (0.0,10.0)
+p = [1.5,1.0,3.0,1.0]
+prob1 = ODEProblem(f1,u0,tspan,p)
+sol = solve(prob1,Tsit5())
+t = collect(linspace(1,10,10))
+randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
+data = convert(Array,randomized)
+priors = [Truncated(Normal(1.5,1),0,2),Truncated(Normal(1.0,1),0,1.5),
+          Truncated(Normal(3.0,1),0,4),Truncated(Normal(1.0,1),0,2)]
+
+bayesian_result = abc_inference(prob1,Tsit5(),t,data,priors; num_samples=500)
+
+meanvals = mean(bayesian_result.parameters, weights(bayesian_result.weights), 1)
+
+@show meanvals
+
+@test meanvals[1] ≈ 1.5 atol=3e-1
+@test meanvals[2] ≈ 1.0 atol=3e-1
+@test meanvals[3] ≈ 3.0 atol=3e-1
+@test meanvals[4] ≈ 1.0 atol=3e-1

test/runtests.jl

@@ -6,4 +6,5 @@ tic()
 @time @testset "Stan_String" begin include("stan_string.jl") end
 @time @testset "Stan" begin include("stan.jl") end
 @time @testset "Turing" begin include("turing.jl") end # Doesn't work on v0.6
+@time @testset "ABC" begin include("abc.jl") end
 toc()