Getting started

Chains type

MCMCChains.ChainsType
Chains

Parameters:

  • value: An AxisArray object with axes iter × var × chains
  • logevidence : A field containing the logevidence.
  • name_map : A NamedTuple mapping each variable to a section.
  • info : A NamedTuple containing miscellaneous information relevant to the chain.

The info field can be set using setinfo(c::Chains, n::NamedTuple).

source

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:

MCMCChains.replacenamesFunction
replacenames(chains::Chains, dict::AbstractDict)

Replace parameter names by creating a new Chains object that shares the same underlying data.

Examples

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
source

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

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

Return the parameters with the same name sym, but have a different index. Bracket indexing format in the form of :sym[index] is assumed by default. Use index_type=:dot for parameters with dot indexing, i.e. :sym.index.

If the chain contains a parameter of name :sym it will be returned as well.

Example

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")
source

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)
6×4 DataFrame
Rowiterationchainab
Int64Int64Float64Float64
1110.5914240.0238958
2210.637220.117082
3310.3847310.152856
4120.5343290.0840422
5220.0200470.0846894
6320.1960540.576018

See also ?DataFrame and ?Array for more help.

Sampling Chains

MCMCChains overloads several sample methods as defined in StatsBase:

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

Sample n samples from the pooled (!) chain chn.

The keyword arguments replace and ordered determine whether sampling is performed with replacement and whether the sample is ordered, respectively. If specified, sampling probabilities are proportional to weights wv.

Note

If chn contains multiple chains, they are pooled (i.e., appended) before sampling. This ensures that even in this case exactly n samples are returned:

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

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

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.