# Setting up multiple chains in Turing.jl

**URL:** <https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806>\
**Category:** General Usage\
**Created:** [February 15, 2019, 5:03am UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806 "2019-02-15T05:03:34Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![mkarikom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkarikom/32/7096_2.png) [@mkarikom](https://discourse.julialang.org/u/mkarikom)\
**Post date:** [February 15, 2019, 5:03am UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/1 "2019-02-15T05:03:35Z")

</div>

I’ve run a couple Metropolis chains and produced a Chains object for convergence diagnosis, but I’m pretty sure there is a better way to do this:

fit1 = sample(SBM(Y = Y, Pi = π₀, B = B₀),MH(1000,))  
fit2 = sample(SBM(Y = Y, Pi = π₀, B = B₀),MH(1000,))  
fit = Chains(cat( hcat(fit1[:z]…)‘, hcat(fit1[:z]…)’, dims = 3))

---

<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 15, 2019, 8:46am UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/2 "2019-02-15T08:46:51Z")

</div>

I opened a issue to request internal handling multiple chains in parallel. In the meantime, I have been using `pmap` to run multiple chains.

```julia
using Distributed
Nchains = 4
addprocs(Nchains)
@everywhere begin
    #put all of your supporting code here
    push!(LOAD_PATH,"path/myfile.jl")
    using Turing,myModule
end

@everywhere @model mymodel() begin = 
    #stuff here
end

chain = pmap(x->sample(SBM(Y = Y, Pi = π₀, B = B₀),MH(1000,)),1:Nchains)
chains = removeBurnin(chain)
r̂ = gelmandiag(chains)

```

There is probably a more elegant solution for handling the chains and removing the burnin samples, but this is what I use:

```julia
function removeBurnin(ch::Array{Chain{AbstractRange{Int64}},1},Nadapt=0)
    nch = Chain()
    Nchains = length(ch)
    v = ch[1].value
    dims = (size(v,1)-Nadapt,size(v,2),Nchains)
    nch.value = fill(0.0,dims)
    value2 = []
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(ch)
        nch.value[:,:,i] = c.value[rng,:,:]
        push!(value2,c.value2[rng])
    end
    nch.value2 =reshape(vcat(value2...),dims[1],dims[3])
    nch.names = ch[1].names
    nch.range = 1:dims[1]
    nch.chains = 1:Nchains
    return nch
end

```

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [February 15, 2019, 4:24pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/3 "2019-02-15T16:24:13Z")

</div>

Won’t simply creating a new Chains object work, i.e.:

```julia
draws = adapt_cycles+1:samples
chn2 = MCMCChain.Chains(chn.value[draws,:,:], names=chn.names)

```

where `chn` is the uncorrected chn?

It would also be interesting to know if anyone has experimented with the number of required `draws`. I’m wondering if the numbers have gotten inflated because the results - as in describe(chn) - of the uncorrected chn were usually off. Slowly I’m testing that by reducing the number of draws in TuringModels.jl.

---

<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 15, 2019, 6:04pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/4 "2019-02-15T18:04:34Z")

</div>

I believe this correct for the simple case with one chain. The problem is that Turing does not handle multiple chains internally. My function restructures an array of chains from `pmap` into a single chain. At some point, all the seemingly superfluous information was necessary to generate the plots (but not the diagnostics). That may have changed since I created the function.

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [February 15, 2019, 7:29pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/5 "2019-02-15T19:29:12Z")

</div>

Thanks, for now I will keep a copy of your removeBurnin around if you don’t mind.

---

<div class="post-metadata">

**Author:** ![mkarikom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkarikom/32/7096_2.png) [@mkarikom](https://discourse.julialang.org/u/mkarikom)\
**Post date:** [February 15, 2019, 7:44pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/6 "2019-02-15T19:44:51Z")

</div>

This is very helpful, thanks

---

<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 16, 2019, 7:44am UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/7 "2019-02-16T07:44:52Z")

</div>

Please do, Rob. I owe you one… or ten.

You are welcome mkarikom. I’m glad it’s helpful.

---

<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:** [March 6, 2019, 12:33pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/8 "2019-03-06T12:33:52Z")

</div>

Rob, I was wondering if you or someone else might now how to adapt removeBurnin to the new Chains object? Here is what I have so far:

```julia
function removeBurnin(chainArray::Array{<:Chains,1},Nadapt)
    Nchains = length(chainArray)
    ch = chainArray[1]
    v = ch.value
    dims = (size(v,1)-Nadapt,size(v,2),Nchains)
    value = Array{Union{Missing,Real}}(undef,dims...)
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(chainArray)
        value[:,:,i] = c.value[rng,:,:]
    end
    nch = Chains(value,ch.name_map.parameters)
    return nch
end

```

Error Message:

```julia
ArgumentError: the length of each axis must match the corresponding size of data
AxisArray(::Array{Union{Missing, Float64},3}, ::Tuple{Axis{:iter,StepRange{Int64,Int64}},Axis{:var,Array{String,1}},Axis{:chain,Array{Symbol,1}}}) at core.jl:222
AxisArray(::Array{Union{Missing, Float64},3}, ::Axis{:iter,StepRange{Int64,Int64}}, ::Vararg{Union{Axis, AbstractArray{T,1} where T},N} where N) at core.jl:215
#Chains#12(::Int64, ::Int64, ::Float64, ::NamedTuple{(),Tuple{}}, ::Type, ::Array{Union{Missing, Float64},3}, ::Array{String,1}, ::Dict{Symbol,Array{Any,1}}) at chains.jl:93
Type at chains.jl:28 [inlined]
Chains(::Array{Union{Missing, Float64},3}, ::Array{String,1}) at chains.jl:28
top-level scope at none:0

```

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 2:50pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/9 "2019-03-06T14:50:58Z")

</div>

Hi Chris, I’ll take a look this morning.

---

<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:** [March 6, 2019, 3:04pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/10 "2019-03-06T15:04:46Z")

</div>

Thanks, Rob. From what I can tell, the issue is that value contains samples from the parameters and info from the MCMC sampler, such as lp.

As a workaround I removed the MCMC info:

```julia
function removeBurnin(chainArray::Array{<:Chains,1},Nadapt)
    Nchains = length(chainArray)
    ch = chainArray[1]
    v = ch.value
    parms = ch.name_map.parameters
    dims = (size(v,1)-Nadapt,length(parms),Nchains)
    value = Array{Union{Missing,Real}}(undef,dims...)
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(chainArray)
        value[:,:,i] = c.value[rng,parms,:]
    end
    nch = Chains(value,parms)
    return nch
end

```

The downside is that the MCMC info is not included. I usually don’t use it, but there could be situations in which it might be needed for diagnostics. So it might be nice to have the info.

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 4:55pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/11 "2019-03-06T16:55:47Z")

</div>

Hi Chris,

Given that you don’t use the :internals section that often, would the following work (provided the `chn` below is defined as the output of a single run)?

```julia
function removeBurnin(chainArray::Array{<:Chains,1},Nadapt, ::Val{:parameters})
    Nchains = length(chainArray)
    ch = chainArray[1]
    v = ch.value
    parms = ch.name_map.parameters
    dims = (size(v,1)-Nadapt,length(parms),Nchains)
    value = Array{Union{Missing,Real}}(undef,dims...)
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(chainArray)
        value[:,:,i] = c.value[rng,parms,:]
    end
    nch = Chains(value,parms)
    return nch
end

function removeBurnin(chainArray::Array{<:Chains,1},Nadapt, ::Val{:internals})
    Nchains = length(chainArray)
    ch = chainArray[1]
    v = ch.value
    parms = ch.name_map.internals
    dims = (size(v,1)-Nadapt,length(parms),Nchains)
    value = Array{Union{Missing,Real}}(undef,dims...)
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(chainArray)
        value[:,:,i] = c.value[rng,parms,:]
    end
    nch = Chains(value,parms)
    return nch
end

ca = Vector{MCMCChains.Chains}(undef, 4)
for i in 1:4
  ca[i] = chn
end

describe(ca[2])
describe(ca[2], section=:internals)

parameters_chain = removeBurnin(ca, 1000, Val(:parameters))

describe(parameters_chain)

internals_chain = removeBurnin(ca, 1000, Val(:internals))

describe(internals_chain)

```

It would not be too hard to recreate the normal structure (i.e. a Chains object with both sections present,  
just let us know what your preference is, I stuck to how you did this as you are familiar with it). I’ve copied the master ( @cpfiffer )!

Other ways of removing the adaptation draws you can see in examples such as [m10.4t.jl](https://statisticalrethinkingjulia.github.io/TuringModels.jl/latest/10/m10.04t/):

```julia
draws = 1001:2000
posterior2 = Chains(posterior[draws,:,:], :parameters)

```

you would have to apply to the elements of the chain array.

I love these sections! Especially in multilevel models, e.g. [m10.4s.jl](https://github.com/StanJulia/StanMCMCChains.jl/blob/master/examples/m10.4s/m10.4s.jl).

---

<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:** [March 6, 2019, 5:28pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/12 "2019-03-06T17:28:39Z")

</div>

Thanks, Rob! I think both of these look good.

---

<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:** [March 6, 2019, 5:30pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/13 "2019-03-06T17:30:44Z")

</div>

Of course, in practice, it might be convenient to have both sections present.

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 6:17pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/14 "2019-03-06T18:17:34Z")

</div>

Agreed, I’ll make a new version.

---

<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:** [March 6, 2019, 9:09pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/15 "2019-03-06T21:09:13Z")

</div>

By the way, I noticed that in some cases removeBurnin does not assign the values to the correct parameters. Here is a simple example:

```julia
using Distributed
addprocs(4)
@everywhere using Turing,Distributions

@everywhere @model model() = begin
    a ~ Normal(-10,.1)
    b ~ Normal(-8,.1)
    c ~ Normal(-6,.1)
    d ~ Normal(-4,.1)
    return nothing
end

Nsamples = 2000
Nadapt = 1000
δ = .85
specs = NUTS(Nsamples,Nadapt,δ)
chainArray = pmap(x->sample(model(),specs),1:4)
chain = removeBurnin(chainArray,Nadapt)

```

Output:

```julia
parameters
    Mean SD Naive SE MCSE ESS
**a -6.0009 0.0989 0.0016 0.0015 1000**
b -8.0006 0.0996 0.0016 0.0013 1000
**c -9.9994 0.1007 0.0016 0.0016 1000**
d -4.0009 0.1007 0.0016 0.0015 1000

```

I think it has something to do with slicing the array c.value.

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 9:18pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/16 "2019-03-06T21:18:58Z")

</div>

Yip, I think in the new MCMCChains this mapping is via AxisArrays and can be tricky. I just finished another, somewhat similar exercise to split an existing section in 2 sections and ran into this as well. Now looking at removeBurnin().

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 9:45pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/17 "2019-03-06T21:45:27Z")

</div>

Hi Chris, can you have a look at this setup?

```julia
using Distributed
addprocs(4)
@everywhere using Turing

@everywhere @model model() = begin
    a ~ Normal(-10,.1)
    b ~ Normal(-8,.1)
    c ~ Normal(-6,.1)
    d ~ Normal(-4,.1)
    return nothing
end

function flatten_name_map(chn::MCMCChains.AbstractChains)
  pn = String[]
  parms = values(chn.name_map)
  for i in 1:length(parms)
    for j in 1:length(parms[i])
      append!(pn, [String(parms[i][j])])
    end
  end
  pn
end

function construct_a3d(chn::MCMCChains.AbstractChains)
  cnames = flatten_name_map(chn)
  d, p, c = size(chn.value)
  a3d = fill(0.0, d, p, c);

  for (i, par) in enumerate(cnames)
    a3d[:, i, 1] = reshape(chn[par], d)
  end
  a3d
end

function removeBurnin(chainArray::Array{<:Chains,1},Nadapt)
    Nchains = length(chainArray)
    ch = chainArray[1]
    v = ch.value
    parms = flatten_name_map(ch)
    dims = (size(v,1)-Nadapt,length(parms),Nchains)
    value = Array{Union{Missing,Real}}(undef,dims...)
    rng = (Nadapt+1):size(v,1)
    for (i,c) in enumerate(chainArray)
        a3d = construct_a3d(c)
        value[:,:,i] = a3d[rng, :, :]
    end
    MCMCChains.Chains(value,
      Symbol.(parms),
      Dict(
        :parameters => Symbol.(values(ch.name_map.parameters)),
        :internals => Symbol.(values(ch.name_map.internals))
      )
    )
end

Nsamples = 2000
Nadapt = 1000
δ = .85
specs = NUTS(Nsamples,Nadapt,δ)
chainArray = pmap(x->sample(model(),specs),1:4)
chain = removeBurnin(chainArray,Nadapt)

```

---

<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:** [March 6, 2019, 10:07pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/18 "2019-03-06T22:07:21Z")

</div>

Rob, this looks good as far as I can tell. Changing to AxisArrays seems to have complicated the manipulation of the chain object quite a bit. Thanks for your help!

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [March 6, 2019, 10:23pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/19 "2019-03-06T22:23:31Z")

</div>

I think you guys have a lot of moving machinery here. MCMCChains should support all this behavior. I added a `chainscat` function that adds chains together a little while ago — you should be able to do something like

```julia
chns = reduce(chainscat, pmap(x->sample(model(),specs),1:4))
subset = chns[Nadapt:Nsamples, :, :]

```

Here’s the output from `describe(subset)`. Guess I should fix that `quantile` warning too, while I’m at it.

```julia
Log evidence = 0.0
Iterations = 1000:2000
Thinning interval = 1
Chains = Chain1, Chain2, Chain3, Chain4
Samples per chain = 1001
parameters = c, b, a, d

┌ Warning: `quantile(v::AbstractArray{<:Real})` is deprecated, use `quantile(v, [0.0, 0.25, 0.5, 0.75, 1.0])` instead.
│ caller = (::getfield(MCMCChains, Symbol("##102#104")){Chains{Union{Missing, Float64},Float64,NamedTuple{(:parameters,),Tuple{Array{String,1}}},NamedTuple{(:hashedsummary,),Tuple{Base.RefValue{Tuple{UInt64,MCMCChains.ChainSummaries}}}}}})(::String) at none:0
└ @ MCMCChains ./none:0
Empirical Posterior Estimates
───────────────────────────────────────────
parameters
    Mean SD Naive SE MCSE ESS
a -9.9998 0.0980 1.5×10⁻³ 0.0013 1001
b -7.9995 0.1019 1.6×10⁻³ 0.0016 1001
c -5.9986 0.1016 1.6×10⁻³ 0.0016 1001
d -3.9992 0.0992 1.6×10⁻³ 0.0017 1001

Quantiles
───────────────────────────────────────────
parameters
    2.5% 25.0% 50.0% 75.0% 97.5% 
a -10.3729 -10.0661 -9.9998 -9.9317 -9.6298
b -8.3329 -8.0687 -8.0031 -7.9298 -7.6591
c -6.3259 -6.0663 -5.9990 -5.9311 -5.6473
d -4.3360 -4.0632 -3.9993 -3.9342 -3.6764

```

And here’s the full code I used to make that:

```julia
using Distributed
addprocs(4)
@everywhere using Turing,Distributions

@everywhere @model model() = begin
    a ~ Normal(-10,.1)
    b ~ Normal(-8,.1)
    c ~ Normal(-6,.1)
    d ~ Normal(-4,.1)
    return nothing
end

Nsamples = 2000
Nadapt = 1000
δ = .85
specs = NUTS(Nsamples,Nadapt,δ)
chns = reduce(chainscat, pmap(x->sample(model(),specs),1:4))
subset = chns[Nadapt:Nsamples, :, :]

```

---

<div class="post-metadata">

**Author:** ![goedman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goedman/32/217_2.png) [@goedman](https://discourse.julialang.org/u/goedman)\
**Post date:** [March 6, 2019, 11:00pm UTC](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806/20 "2019-03-06T23:00:32Z")

</div>

Definitely looks simpler! I had looked at those cat functions but hadn’t figured out when to use them.

Will study them again! Thanks for your help!

[Next page](https://discourse.julialang.org/t/setting-up-multiple-chains-in-turing-jl/20806.md?page=2)