GitHub

Getting started

Chains type

MCMCChains.Chains — Type
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

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

           a    0.4959    0.2925    0.0084   1200.9445   1349.5633    1.0004   ⋯
           b    0.5049    0.2956    0.0077   1459.2812   1388.6217    1.0014   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           a    0.0245    0.2336    0.5015    0.7479    0.9695
           b    0.0245    0.2399    0.5077    0.7758    0.9768

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

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

     param_1    0.4926    0.2994    0.0219   190.3982   127.9210    1.0154     ⋯
     param_2    0.5139    0.2784    0.0189   216.5128   180.6551    0.9962     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

     param_1    0.0191    0.2334    0.5001    0.7764    0.9529
     param_2    0.0352    0.2718    0.5279    0.7469    0.9725

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

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

     param_1    0.4926    0.2994    0.0219   190.3982   127.9210    1.0154     ⋯
     param_2    4.0000    0.0000       NaN        NaN        NaN       NaN     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

     param_1    0.0191    0.2334    0.5001    0.7764    0.9529
     param_2    4.0000    4.0000    4.0000    4.0000    4.0000

Rename Parameters

Parameter names can be changed with the function replacenames:

MCMCChains.replacenames — Function
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

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           a    0.4758    0.2866    0.0215   165.9821   130.1124    1.0002     ⋯
           b    0.4675    0.2940    0.0169   277.5899   220.9838    0.9932     ⋯
           c    0.4679    0.2858    0.0193   193.5564   185.0956    1.0016     ⋯
           d    0.5162    0.3041    0.0271   130.3717   179.7268    1.0379     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           a    0.0217    0.2177    0.4505    0.6937    0.9651
           b    0.0175    0.2070    0.4482    0.7156    0.9819
           c    0.0243    0.2259    0.4518    0.6724    0.9740
           d    0.0138    0.2743    0.4781    0.7673    0.9823

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

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           a    0.4758    0.2866    0.0215   165.9821   130.1124    1.0002     ⋯
           b    0.4675    0.2940    0.0169   277.5899   220.9838    0.9932     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           a    0.0217    0.2177    0.4505    0.6937    0.9651
           b    0.0175    0.2070    0.4482    0.7156    0.9819

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.namesingroup — Function
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.sample — Method
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.