# ABC / MSM with Turing.jl

**URL:** <https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189>\
**Category:** Statistics\
**Tags:** turing\
**Created:** [February 10, 2022, 6:29pm UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189 "2022-02-10T18:29:08Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [February 10, 2022, 6:29pm UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189/1 "2022-02-10T18:29:08Z")

</div>

I believe that ABC (approximate Bayesian computing) and/or Bayesian MSM (method of simulated moments) are not yet well supported by Turing.jl. I have come up with a simple example that does the basic steps, and it was surprisingly fun and straightforward. In this example, inference is based on a vector of statistics, as is done for certain types of ABC, and as is always done with MSM. The distance measure (likelihood) is an estimate of the asymptotic Gaussian distribution of the statistics, computed by simulation. Here’s the example, which might be of interest for more fully developing ABC/MSM with Turing.

```julia
# second example, with a vector of parameters
using Turing, Statistics, Distributions, Random, LinearAlgebra, StatsPlots
using AdvancedMH

#@views function main()
θ = [2.; 3.]
n = 100 # sample size
S = 100 # number of simulation draws
# the data generating process
function dgp(θ, n)
    [rand.(Exponential(θ[1]),n) rand.(Poisson(θ[2]),n)] 
end
    # summary statistics for estimation
function moments(y)
    sqrt(n) .* [mean(y, dims=1)[:]; std(y, dims=1)[:]]
end
y = dgp(θ, n)
z = moments(y)
zs = zeros(S,size(z,1))

@model function abc(z)
    # create the prior: the product of the following array of marginal priors
    θ ~ arraydist([LogNormal(1.,1.); LogNormal(1.,1.)])
    # sample from the model, at the trial parameter value, and compute statistics
    @inbounds for i = 1:S
        y .= dgp(θ, n)
        zs[i,:] .= moments(y) # simulated summary statistics
    end
    # the asymptotic Gaussian distribution of the statistics
    m = mean(zs, dims=1)[:]
    Σ = Symmetric((1. + 1/S)*cov(zs))
    z ~ MvNormal(m, Σ)
end;

# sample chains, 4 chains of length 1000, with 100 burnins dropped
length = 1000
burnin = 100
chain = sample(abc(z), 
    MH(:θ => AdvancedMH.RandomWalkProposal(MvNormal(zeros(2), 0.25*I))),
    MCMCThreads(), length+burnin, 4)
chain = chain[burnin+1:end,:,:]
display(chain)
p = plot(chain)
display(p)
#end
#main()

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [February 10, 2022, 7:35pm UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189/2 "2022-02-10T19:35:10Z")

</div>

Thanks for the example. I would be interested in more support for ABC methods.

---

<div class="post-metadata">

**Author:** ![jbrea](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbrea/32/3879_2.png) [@jbrea](https://discourse.julialang.org/u/jbrea)\
**Post date:** [February 10, 2022, 8:12pm UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189/3 "2022-02-10T20:12:42Z")

</div>

Nice, thanks for sharing! We discussed at some point similar approaches [here](https://github.com/JuliaApproxInference/ApproxInferenceBase.jl/issues/1), but so far I haven’t found the time to continue this work. Some of it has been done, however, in [KissABC](https://github.com/francescoalemanno/KissABC.jl).

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [February 11, 2022, 7:03am UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189/4 "2022-02-11T07:03:05Z")

</div>

That proposed structure to organize methods looks very useful.

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [February 11, 2022, 3:48pm UTC](https://discourse.julialang.org/t/abc-msm-with-turing-jl/76189/5 "2022-02-11T15:48:49Z")

</div>

I put a few scripts at [https://github.com/mcreel/ABC-MSM-Turing](https://github.com/mcreel/ABC-MSM-Turing)