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:500
Thinning interval = 1
Chains            = 1, 2, 3
Samples per chain = 500
parameters        = a, b

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64

           a    0.4962    0.2906     0.0075    0.0065   1235.9560    0.9994
           b    0.4932    0.2883     0.0074    0.0062   1409.4492    0.9989

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

           a    0.0302    0.2356    0.4861    0.7590    0.9729
           b    0.0221    0.2452    0.4831    0.7434    0.9728

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:100
Thinning interval = 1
Chains            = 1, 2
Samples per chain = 100
parameters        = param_1, param_2

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

     param_1    0.5009    0.2875     0.0203    0.0063   222.4027    0.9923
     param_2    0.5089    0.2832     0.0200    0.0190   224.7700    0.9967

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

     param_1    0.0151    0.2510    0.5056    0.7440    0.9477
     param_2    0.0352    0.2769    0.5048    0.7558    0.9515

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.0167485  0.601748
 0.945817   0.84651
 0.84343    0.371065
 0.163738   0.195013
 0.416416   0.866167

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

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

Iterations        = 1:100
Thinning interval = 1
Chains            = 1, 2
Samples per chain = 100
parameters        = param_1, param_2

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

     param_1    0.5009    0.2875     0.0203    0.0063   222.4027    0.9923
     param_2    4.0000    0.0000     0.0000    0.0000        NaN       NaN

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

     param_1    0.0151    0.2510    0.5056    0.7440    0.9477
     param_2    4.0000    4.0000    4.0000    4.0000    4.0000

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:100
Thinning interval = 1
Chains            = 1, 2
Samples per chain = 100
parameters        = a, b, c, d

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

           a    0.5072    0.2901     0.0205    0.0084   208.1973    0.9918
           b    0.5405    0.2757     0.0195    0.0077   191.8219    0.9928
           c    0.4861    0.2925     0.0207    0.0015   202.8022    0.9978
           d    0.5011    0.2662     0.0188    0.0240   241.3182    0.9989

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

           a    0.0296    0.2669    0.5043    0.7608    0.9777
           b    0.0310    0.3093    0.5425    0.7917    0.9884
           c    0.0153    0.2115    0.5185    0.7421    0.9679
           d    0.0397    0.2821    0.4925    0.7177    0.9553

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:100
Thinning interval = 1
Chains            = 1, 2
Samples per chain = 100
parameters        = a, b
internals         = c, d

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

           a    0.5072    0.2901     0.0205    0.0084   208.1973    0.9918
           b    0.5405    0.2757     0.0195    0.0077   191.8219    0.9928

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

           a    0.0296    0.2669    0.5043    0.7608    0.9777
           b    0.0310    0.3093    0.5425    0.7917    0.9884

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 or MCMCChains.gelmandiag 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{String,Symbol})

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

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.726892202151076; 0.7897987447542649; … ; 0.7007065543867941; 0.28230852428601994],)

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.726892202151076
 0.7897987447542649
 0.8764192497161967
 0.4298823940662093
 0.7007065543867941
 0.28230852428601994

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

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

(a = [0.726892202151076; 0.7897987447542649; … ; 0.7007065543867941; 0.28230852428601994],
 b = [0.5217728481429613; 0.516762893315768; … ; 0.06530586170606023; 0.7861388159619587],)

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

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.0376467  0.285629
 0.494245   0.669787
 0.622376   0.442597
 0.962889   0.574631
 0.0843298  0.327158
 0.709288   0.362666

Array(chn, [:parameters])

6×2 Matrix{Float64}:
 0.0376467  0.285629
 0.494245   0.669787
 0.622376   0.442597
 0.962889   0.574631
 0.0843298  0.327158
 0.709288   0.362666

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.03764674908160326 0.28562924993877314; 0.49424522608710353 0.6697872641396161; 0.622375581088799 0.44259717400617427]
 [0.9628890746643262 0.5746307472028429; 0.08432978813107983 0.3271577375512067; 0.7092878338443243 0.36266551726565144]

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

A[2]

3×2 Matrix{Float64}:
 0.962889   0.574631
 0.0843298  0.327158
 0.709288   0.362666

Similarly, for DataFrames:

using DataFrames

DataFrame(chn)

6 rows × 4 columns

See also ?DataFrame and ?Array for more help.

Sampling Chains

MCMCChains overloads several sample methods as defined in StatsBase:

sample(chn::Chains, [wv::AbstractWeights,] n; replace=true, kwargs...)
sample(rng::Random.AbstractRNG, chn::Chains, [wv::AbstractWeights,] n; kwargs...)

subset(chn::Chains, samples)

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.