Getting started

Chains type

Chains

Indexing and parameter Names

Chains can be constructed with parameter names. For example, to create a chains object with

• 500 samples,
• 2 parameters (named a and b)
• 3 chains

use

val = rand(500, 2, 3)
chn = Chains(val, [:a, :b])

Chains MCMC chain (500×2×3 Array{Float64, 3}):

Iterations        = 1:1:500
Number of chains  = 3
Samples per chain = 500
parameters        = a, b

Use `describe(chains)` for summary statistics and quantiles.

By default, parameters will be given the name param_i, where i is the parameter number:

chn = Chains(rand(100, 2, 2))

Chains MCMC chain (100×2×2 Array{Float64, 3}):

Iterations        = 1:1:100
Number of chains  = 2
Samples per chain = 100
parameters        = param_1, param_2

Use `describe(chains)` for summary statistics and quantiles.

We can set and get indexes for parameter 2:

chn_param2 = chn[1:5,2,:];

2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 1:1:5
    :chain, 1:2
And data, a 5×2 Matrix{Float64}:
 0.162927  0.388083
 0.716976  0.334046
 0.417752  0.286943
 0.920146  0.552454
 0.647637  0.735593

chn[:,2,:] = fill(4, 100, 1, 2)
chn

Chains MCMC chain (100×2×2 Array{Float64, 3}):

Iterations        = 1:1:100
Number of chains  = 2
Samples per chain = 100
parameters        = param_1, param_2

Use `describe(chains)` for summary statistics and quantiles.

Rename Parameters

Parameter names can be changed with the function replacenames:

replacenames(chains::Chains, dict::AbstractDict)

julia> chn = Chains(rand(100, 2, 2), ["one", "two"]);

julia> chn2 = replacenames(chn, "one" => "A");

julia> names(chn2)
2-element Vector{Symbol}:
 :A
 :two

julia> chn3 = replacenames(chn2, Dict("A" => "one", "two" => "B"));

julia> names(chn3) 
2-element Vector{Symbol}:
 :one
 :B

Sections

Chains parameters are sorted into sections that represent groups of parameters, see MCMCChains.group. By default, every chain contains a parameters section, to which all unassigned parameters are assigned to. Chains can be assigned a named map during construction:

chn = Chains(rand(100, 4, 2), [:a, :b, :c, :d])

Chains MCMC chain (100×4×2 Array{Float64, 3}):

Iterations        = 1:1:100
Number of chains  = 2
Samples per chain = 100
parameters        = a, b, c, d

Use `describe(chains)` for summary statistics and quantiles.

The MCMCChains.set_section function returns a new Chains object:

chn2 = set_section(chn, Dict(:internals => [:c, :d]))

Chains MCMC chain (100×4×2 Array{Float64, 3}):

Iterations        = 1:1:100
Number of chains  = 2
Samples per chain = 100
parameters        = a, b
internals         = c, d

Use `describe(chains)` for summary statistics and quantiles.

Note that only a and b are being shown. You can explicity retrieve an array of the summary statistics and the quantiles of the :internals section by calling describe(chn; sections = :internals), or of all variables with describe(chn; sections = nothing). Many functions such as MCMCChains.summarize support the sections keyword argument.

Groups of parameters

You can access the names of all parameters in a chain that belong to the group name by using

namesingroup(chains::Chains, sym::Union{AbstractString,Symbol}; index_type::Symbol=:bracket)

julia> chn = Chains(rand(100, 2, 2), ["A[1]", "A[2]"]);

julia> namesingroup(chn, :A)
2-element Vector{Symbol}:
 Symbol("A[1]")
 Symbol("A[2]")

julia> # Also works for specific elements.
       namesingroup(chn, Symbol("A[1]"))
1-element Vector{Symbol}:
 Symbol("A[1]")

julia> chn = Chains(rand(100, 3, 2), ["A.1", "A.2", "B"]);

julia> namesingroup(chn, :A; index_type=:dot)
2-element Vector{Symbol}:
 Symbol("A.1")
 Symbol("A.2")

The get Function

MCMCChains also provides a get function designed to make it easier to access parameters:

val = rand(6, 3, 1)
chn = Chains(val, [:a, :b, :c]);

x = get(chn, :a)

(a = [0.4406973184362195; 0.06413541593814787; … ; 0.6564549998197146; 0.4260695808953028;;],)

You can also access the variables via getproperty:

x.a

2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 1:1:6
    :chain, 1:1
And data, a 6×1 Matrix{Float64}:
 0.4406973184362195
 0.06413541593814787
 0.45084693164962364
 0.7032412293454546
 0.6564549998197146
 0.4260695808953028

get also accepts vectors of things to retrieve, so you can call

x = get(chn, [:a, :b])

(a = [0.4406973184362195; 0.06413541593814787; … ; 0.6564549998197146; 0.4260695808953028;;], b = [0.94708496660485; 0.20744026514507408; … ; 0.8847254154251538; 0.33010442589688427;;])

Saving and Loading Chains

Like any Julia object, a Chains object can be saved using Serialization.serialize and loaded back by Serialization.deserialize as identical as possible. Note, however, that in general this process will not work if the reading and writing are done by different versions of Julia, or an instance of Julia with a different system image. You might want to consider JLSO for saving metadata such as the Julia version and the versions of all packages installed as well.

using Serialization

serialize("chain-file.jls", chn)
chn2 = deserialize("chain-file.jls")

The MCMCChainsStorage.jl package also provides the ability to serialize/deserialize a chain to an HDF5 file across different versions of Julia and/or different system images.

Exporting Chains

A few utility export functions have been provided to convert Chains objects to either an Array or a DataFrame:

chn = Chains(rand(3, 2, 2), [:a, :b])

Array(chn)

6×2 Matrix{Float64}:
 0.591424  0.0238958
 0.63722   0.117082
 0.384731  0.152856
 0.534329  0.0840422
 0.020047  0.0846894
 0.196054  0.576018

Array(chn, [:parameters])

6×2 Matrix{Float64}:
 0.591424  0.0238958
 0.63722   0.117082
 0.384731  0.152856
 0.534329  0.0840422
 0.020047  0.0846894
 0.196054  0.576018

By default chains are appended. This can be disabled by using the append_chains keyword argument:

A = Array(chn, append_chains=false)

2-element Vector{Matrix{Float64}}:
 [0.5914241067613842 0.023895751073906513; 0.6372198453275196 0.11708206509997343; 0.3847313444636996 0.15285574309518535]
 [0.5343285749143705 0.08404220391148931; 0.020047016741082002 0.08468937966306267; 0.19605406287287575 0.5760176829441074]

which will return a matrix for each chain. For example, for the second chain:

A[2]

3×2 Matrix{Float64}:
 0.534329  0.0840422
 0.020047  0.0846894
 0.196054  0.576018

Similarly, for DataFrames:

using DataFrames

DataFrame(chn)

See also ?DataFrame and ?Array for more help.

Sampling Chains

MCMCChains overloads several sample methods as defined in StatsBase:

sample([rng,] chn::Chains, [wv::AbstractWeights,] n; replace=true, ordered=false)

julia> chn = Chains(randn(11, 4, 3));

julia> size(sample(chn, 7)) == (7, 4, 1)
true

See ?sample for additional help on sampling. Alternatively, you can construct and sample from a kernel density estimator using KernelDensity.jl, see test/sampling_tests.jl.