Accessing samples

FlexiChains stores data in a dict-of-array format (a column-oriented form, similar to dataframes). This means that it is very easy to access all the data for a given variable, but accessing all the data for a given sample (i.e., iteration / chain number) is a bit more involved.

Consider the following example:

using FlexiChains, Turing

@model function f()
    x ~ Normal()
    y ~ Uniform(1, 2)
end

chn = sample(f(), Prior(), 10; chain_type=VNChain, progress=false)

╭─FlexiChain (10 iterations, 1 chain) ─────────────────────────────────────────╮
│ ↓ iter  = 1:10                                                               │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (2) ── VarName                                                    │
│  Float64  x, y                                                               │
│                                                                              │
│ Extras (3)                                                                   │
│  Float64  logprior, loglikelihood, logjoint                                  │
╰──────────────────────────────────────────────────────────────────────────────╯

Here, we have two variables x and y. FlexiChains makes it very easy to access all the samples for x, or all the samples for y, and so on:

chn[@varname(x)]

┌ 10×1 DimArray{Float64, 2} Parameter(x) ┐
├────────────────────────────────────────┴───────────── dims ┐
  ↓ iter Sampled{Int64} 1:10 ForwardOrdered Regular Points,
  → chain Sampled{Int64} 1:1 ForwardOrdered Regular Points
└────────────────────────────────────────────────────────────┘
  ↓ →   1
  1     0.0438116
  2     0.0206356
  3    -0.0249099
  4     0.903717
  ⋮    
  7     0.828689
  8     0.693876
  9    -1.32693
 10    -0.752804

However, suppose you want to access the samples for a given iteration and chain number. For example, you might want to know what happened at the fifth iteration.

One way is to subset a chain. The issue is that that gives you just another chain back:

subsetted = chn[iter=5, chain=1]

╭─FlexiChain (1 iteration, 1 chain) ───────────────────────────────────────────╮
│ ↓ iter  = 5:5                                                                │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (2) ── VarName                                                    │
│  Float64  x, y                                                               │
│                                                                              │
│ Extras (3)                                                                   │
│  Float64  logprior, loglikelihood, logjoint                                  │
╰──────────────────────────────────────────────────────────────────────────────╯

and to access the data you will still need to index into it, which is rather awkward:

(x=only(subsetted[@varname(x)]), y=only(subsetted[@varname(y)]))

(x = 0.15469950740626526, y = 1.9748221242807371)

More generally speaking, the operation above can be viewed as a transformation from column-oriented data to row-oriented data, or from a dict-of-arrays to an array-of-dicts. FlexiChains provides a few convenient ways for you to do this transformation.

values_at

The function values_at is at the core of this transformation. Given iteration and chain indices (as keyword arguments), it returns some container that holds all the values for a given iteration.

vs = FlexiChains.values_at(chn; iter=5, chain=1)

ParamsWithStats
 ├─ params
 │  VarNamedTuple
 │  ├─ x => 0.15469950740626526
 │  └─ y => 1.9748221242807371
 └─ stats
    ├─ logprior = -0.9309045020005433
    ├─ loglikelihood = 0.0
    └─ logjoint = -0.9309045020005433

In the case of a chain sampled from Turing, the returned container is DynamicPPL.ParamsWithStats, which separately stores the parameters and the stats as a VarNamedTuple and NamedTuple respectively. This is a high-fidelity representation of the data, and is exactly what you get when sampling with Turing.jl (for example, if you call sample(...; chain_type=Any), you will get an array of ParamsWithStats objects).

Info

This is accomplished by storing a skeletal VarNamedTuple for each sample in the chain; if you are interested, see the DynamicPPL docs for more info.

The main benefit of this is that you can feed this right back into Turing's API. For example, to initialise MCMC sampling from the fifth sample, you can write:

pws = FlexiChains.values_at(chn; iter=5, chain=1)
sample(f(), NUTS(), 10; initial_params=InitFromParams(pws), progress=false);

╭─FlexiChain (10 iterations, 1 chain) ─────────────────────────────────────────╮
│ ↓ iter  = 6:15                                                               │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (2) ── VarName                                                    │
│  Float64  x, y                                                               │
│                                                                              │
│ Extras (14)                                                                  │
│  Int64    n_steps, tree_depth                                                │
│  Bool     is_accept, numerical_error                                         │
│  Float64  acceptance_rate, log_density, hamiltonian_energy,                  │
│           hamiltonian_energy_error, max_hamiltonian_energy_error, step_size, │
│           nom_step_size, logprior, loglikelihood, logjoint                   │
╰──────────────────────────────────────────────────────────────────────────────╯

Sometimes, though, one might want different output formats. There is some support for converting to NamedTuple or AbstractDict, by passing an optional type argument.

Warning

Note that conversion to NamedTuple is lossy if you have VarNames that contain indexing or field access syntax (e.g., x[1] or x.a).

FlexiChains.values_at(chn, NamedTuple; iter=5, chain=1)

(x = 0.15469950740626526, y = 1.9748221242807371, logprior = -0.9309045020005433, loglikelihood = 0.0, logjoint = -0.9309045020005433)

FlexiChains.values_at(chn, Dict; iter=5, chain=1)

Dict{Union{FlexiChains.Parameter{var"#s77"}, FlexiChains.Extra} where var"#s77"<:VarName, Any} with 5 entries:
  Extra(:logjoint)      => -0.930905
  Extra(:loglikelihood) => 0.0
  Parameter(y)          => 1.97482
  Parameter(x)          => 0.1547
  Extra(:logprior)      => -0.930905

Parameters only

Often the stats are not of very much interest, and you just want the parameters. In that case, you can use parameters_at:

For Turing-sampled chains, this returns a VarNamedTuple, which is just the same as the params field of the ParamsWithStats object.

FlexiChains.parameters_at(chn; iter=5, chain=1)

VarNamedTuple
├─ x => 0.15469950740626526
└─ y => 1.9748221242807371

Arrays of samples

To get more than one sample, you can pass arrays of indices.

FlexiChains.parameters_at(chn; iter=[5, 6], chain=1)

┌ 2-element DimArray{DynamicPPL.VarNamedTuples.VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
  ↓ iter Sampled{Int64} [5, 6] ForwardOrdered Irregular Points
└──────────────────────────────────────────────────────────────────────────────┘
 5  VarNamedTuple(x = 0.1547, y = 1.97482)
 6  VarNamedTuple(x = -0.153802, y = 1.8557)

All of DimensionalData's selectors also work:

FlexiChains.parameters_at(chn; iter=Not(1..8), chain=1)

┌ 2-element DimArray{DynamicPPL.VarNamedTuples.VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
  ↓ iter Sampled{Int64} [9, 10] ForwardOrdered Irregular Points
└──────────────────────────────────────────────────────────────────────────────┘
  9  VarNamedTuple(x = -1.32693, y = 1.53784)
 10  VarNamedTuple(x = -0.752804, y = 1.05915)

This also means that you can convert a FlexiChain back into an Array of ParamsWithStats objects by passing : for both the iteration and chain indices, which is essentially the inverse of what sample does (it bundles the array into a FlexiChain). In fact, : is the default for both of these keyword arguments, so you can just write:

FlexiChains.values_at(chn)

┌ 10×1 DimArray{DynamicPPL.ParamsWithStats{DynamicPPL.VarNamedTuples.VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}, 2} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
  ↓ iter Sampled{Int64} 1:10 ForwardOrdered Regular Points,
  → chain Sampled{Int64} 1:1 ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────┘
  ↓ →  …  1
  1        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.0438116, y = 1.24075), (logprior = -0.919898, loglikelihood = 0.0, logjoint = -0.919898))
  2        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.0206356, y = 1.41781), (logprior = -0.919151, loglikelihood = 0.0, logjoint = -0.919151))
  3        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = -0.0249099, y = 1.9249), (logprior = -0.919249, loglikelihood = 0.0, logjoint = -0.919249))
  4        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.903717, y = 1.19199), (logprior = -1.32729, loglikelihood = 0.0, logjoint = -1.32729))
  ⋮    ⋱  
  7        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.828689, y = 1.7057), (logprior = -1.2623, loglikelihood = 0.0, logjoint = -1.2623))
  8        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.693876, y = 1.61436), (logprior = -1.15967, loglikelihood = 0.0, logjoint = -1.15967))
  9        ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = -1.32693, y = 1.53784), (logprior = -1.79931, loglikelihood = 0.0, logjoint = -1.79931))
 10    …   ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = -0.752804, y = 1.05915), (logprior = -1.2023, loglikelihood = 0.0, logjoint = -1.2023))

Drawing random samples

To get a random sample from the chain, you can use rand:

rand(chn)

ParamsWithStats
 ├─ params
 │  VarNamedTuple
 │  ├─ x => 0.04381160928715135
 │  └─ y => 1.240750618186261
 └─ stats
    ├─ logprior = -0.9198982617588378
    ├─ loglikelihood = 0.0
    └─ logjoint = -0.9198982617588378

This method follows Julia's conventions closely, so you can have an optional rng first argument, and you can also specify the number of samples to draw:

using Random: Xoshiro
rand(Xoshiro(468), chn, 2, 2)

2×2 Matrix{DynamicPPL.ParamsWithStats{DynamicPPL.VarNamedTuples.VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}}:
 ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.0438116, y = 1.24075), (logprior = -0.919898, loglikelihood = 0.0, logjoint = -0.919898))  …  ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = -0.153802, y = 1.8557), (logprior = -0.930766, loglikelihood = 0.0, logjoint = -0.930766))
 ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = 0.0438116, y = 1.24075), (logprior = -0.919898, loglikelihood = 0.0, logjoint = -0.919898))     ParamsWithStats{VarNamedTuple{(:x, :y), Tuple{Float64, Float64}}, @NamedTuple{logprior::Float64, loglikelihood::Float64, logjoint::Float64}}(VarNamedTuple(x = -1.32693, y = 1.53784), (logprior = -1.79931, loglikelihood = 0.0, logjoint = -1.79931))

If you only want parameters, pass the parameters_only=true keyword argument:

rand(chn, parameters_only=true)

VarNamedTuple
├─ x => 0.15469950740626526
└─ y => 1.9748221242807371

Using VarNamedTuples

All of the above functions, when run with Turing-sampled chains, return VarNamedTuples or containers thereof. To access the values, you need to index with a VarName:

vnt = rand(chn, parameters_only=true)
vnt[@varname(x)]

0.15469950740626526

Please see the Turing docs for more details on working with VarNamedTuples.

Docstrings

FlexiChains.values_at Function

FlexiChains.values_at(chn::FlexiChain; iter=:, chain=:)
FlexiChains.values_at(chn::FlexiChain, ::Type{Tout}; iter=:, chain=:)

Extract all values from the chain corresponding to a particular set of MCMC iterations(s).

The iter and chain keyword arguments can be anything used to index into the respective dimensions of a FlexiChain, such as an integer, a vector of integers, a range, or a DimensionalData.jl selector. Both default to : (i.e. all iterations / all chains).

In particular, you can convert the entire chain into a DimMatrix of the desired output type by calling this function with no keyword arguments.

To get only the parameter keys, use FlexiChains.parameters_at.

The output type can be specified with an (optional) positional argument. Possible options are:

• unspecified: uses the structure stored in the chain to determine the output type. Specifically, for chains sampled with Turing.jl, the stored structure will be VarNamedTuple, and consequently the output of parameters_at will be a DynamicPPL.ParamsWithStats per iteration (which stores the parameters as a VarNamedTuple, and the stats/extras as a NamedTuple). For chains that do not have a stored structure, the default output will be an OrderedDict.

• Tout <: AbstractDict: returns a dictionary mapping ParameterOrExtra{TKey} to their values. This is the most faithful representation of the data in the chain.

• Tout = NamedTuple, or Tout = ComponentArray: attempts to convert every key name to a Symbol, which is used as the field name in the output NamedTuple or ComponentArray.

Using NamedTuple or ComponentArray

This will throw an error if any key cannot be converted to a Symbol, or if there are duplicate key names after conversion. If you have parameter names that convert to the same Symbol, you can either use OrderedDict, subset the chain before calling this function, or rename your parameters. Furthermore, please be aware that this is a lossy conversion as it does not retain information about whether a key is a parameter or an extra.

For order-sensitive output types, such as OrderedDict, the keys are returned in the same order as they are stored in the FlexiChain. This also corresponds to the order returned by keys(chn).

source

FlexiChains.parameters_at Function

FlexiChains.parameters_at(chn::FlexiChain; iter=:, chain=:)
FlexiChains.parameters_at(chn::FlexiChain, ::Type{Tout}; iter=:, chain=:)

Extract all parameter values from the chain corresponding to a particular set of MCMC iteration(s), discarding non-parameter (i.e. Extra) keys.

The iter and chain keyword arguments can be anything used to index into the respective dimensions of a FlexiChain, such as an integer, a vector of integers, a range, or a DimensionalData.jl selector. Both default to : (i.e. all iterations / all chains).

In particular, you can convert the entire chain into a DimMatrix of the desired output type by calling this function with no keyword arguments.

To get all keys (not just parameters), use FlexiChains.values_at.

The output type for each iteration can be specified with the first (optional) positional argument. Possible options are:

• unspecified: uses the structure stored in the chain to determine the output type. Specifically, for chains sampled with Turing.jl, the stored structure will be VarNamedTuple, and consequently the output of parameters_at will be a VarNamedTuple per iteration. For chains that do not have a stored structure, the default output will be an OrderedDict.

• Tout <: AbstractDict: returns a dictionary mapping TKey to their values.

• Tout = NamedTuple or Tout <: ComponentArray: attempts to convert every parameter name to a Symbol, which is used as the field name in the output NamedTuple or ComponentArray.

Using NamedTuple or ComponentArray

This will throw an error if any key cannot be converted to a Symbol, or if there are duplicate key names after conversion. If you have parameter names that convert to the same Symbol, you can either use OrderedDict, subset the chain before calling this function, or rename your parameters. Furthermore, please be aware that this is a lossy conversion as it does not retain information about whether a key is a parameter or an extra.

For order-sensitive output types, such as OrderedDict, the parameters are returned in the same order as they are stored in the FlexiChain. This also corresponds to the order returned by FlexiChains.parameters(chn).

source

Base.rand Function

Base.rand(rng::Random.AbstractRNG, chn::FlexiChain, dims::Int...; parameters_only=false)

Sample uniformly from a FlexiChain with replacement.

rand([rng,] chn) returns a single sample, whereas rand([rng,] chn, dims...) returns an Array of samples with dimensions dims.... The parameters_only keyword specifies whether the returned samples include only parameters, or both parameters and extras.

In general, the return type of rand is the same as the return type of FlexiChains.parameters_at or FlexiChains.values_at (but wrapped in an Array if dims is not empty). This means that the return type can depend on the 'structure' stored in the FlexiChain.

For example, if chn::FlexiChain{VarName} was constructed using Turing.jl, rand(chn) will return a DynamicPPL.ParamsWithStats, and rand(chn, parameters_only=true) will return a DynamicPPL.VarNamedTuple.

source