Summarising

In general a FlexiChain contains data in matrices of size (niters, nchains). Often it is useful to summarise this data along one or both dimensions.

FlexiChains therefore allows you to calculate one or more statistics for each variable stored in a FlexiChain. The result is a FlexiSummary object, which can be indexed into in a very similar way to FlexiChains: see the Indexing page for full details, or the examples on this page.

Unsupported data types

Before we launch into the available summary statistics, it is worth mentioning one point about data types. Since FlexiChains allows for storage of completely arbitrary data types, it can contain data for which the mean (or other statistic) is not defined. For example, the mean of String values is not defined. Thus, when calculating mean(chain), any String-valued parameters will be dropped from the result.

If multiple summary statistics are requested (e.g. with summarystats), the key is only dropped if all of them fail. If at least one statistic is successfully computed for a key, that key will be included in the result, with missing values for the statistics that failed.

To give a flavour of how this works, here is an example of a model that generates parameters of different types:

using FlexiChains, Turing

@model function f()
    f ~ Normal()               # float
    v ~ MvNormal(zeros(2), I)  # vector
    s := "a string"            # string
end

chain = sample(f(), MH(), MCMCThreads(), 20, 3; chain_type=VNChain)

╭─FlexiChain (20 iterations, 3 chains) ────────────────────────────────────────╮
│ ↓ iter  = 1:20                                                               │
│ → chain = 1:3                                                                │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64          f                                                          │
│  Vector{Float64}  v (2,)                                                     │
│  String           s                                                          │
│                                                                              │
│ Extras (4)                                                                   │
│  Bool     accepted                                                           │
│  Float64  logprior, loglikelihood, logjoint                                  │
╰──────────────────────────────────────────────────────────────────────────────╯

Over the course of this page we will see what happens to each of these parameters when we try to compute summary statistics.

Overall summary statistics

If you want a quick overview of what's in your chain, summarystats provides a handy selection of commonly used statistics:

StatsBase.summarystats Method

StatsBase.summarystats(
    chain::FlexiChain{TKey};
    split_varnames::Bool=true,
    warn::Bool=true,
) where {TKey}

Compute a standard set of summary statistics for each key in the chain. The statistics include:

• mean (using Statistics.mean)

• standard deviation (Statistics.std)

• Monte Carlo standard error (MCMCDiagnosticTools.mcse)

• bulk effective sample size (MCMCDiagnosticTools.ess)

• tail effective sample size

• R-hat diagnostic (MCMCDiagnosticTools.rhat)

• 5th, 50th (median), and 95th percentiles (Statistics.quantile)

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. The splitting is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types, this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

If any of the statistics cannot be computed for a key, a missing value is returned. If none of the statistics can be computed for a key, that key will be dropped from the resulting FlexiSummary, and a warning issued. The warning can be suppressed by setting warn=false.

source

st = summarystats(chain)

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64  f, v[1], v[2]                                                      │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param     mean     std    mcse  ess_bulk  ess_tail    rhat       q5  …     │
│       f  -0.0267  0.9731  0.1061   92.5113   64.4285  0.9879  -1.6506  …     │
│    v[1]   0.0014  0.8604  0.1010   71.1998   37.5486  0.9941  -1.4380  …     │
│    v[2]   0.2013  1.0990  0.1360   68.9064   23.2807  1.0017  -1.6600  …     │
╰──────────────────────────────────────────────────────────────────────────────╯

Notice that the string-valued parameter s has been dropped from the result, since no statistics could be computed for it. (If you run this in the terminal, you will see a warning about this, so it is not completely silent; it's just not shown in the docs.)

You can index with a variable name (or names!) and the stat dimension:

st[@varname(v[1]), stat=:mean]  # Mean of first element of vector v

0.0013661554201467344

Named statistic shorthand

Passing a Symbol to stat is shorthand for a named lookup: stat=:mean is equivalent to stat=At(:mean). Other selectors follow DimensionalData.jl's usual behaviour; for example, stat=1 selects the first statistic by position.

For more details on indexing, please see the Indexing page.

In the result above, the vector-valued v has been broken up into its individual elements v[1] and v[2]. This happens automatically for chains with VarName keys; you can disable this behaviour by passing split_varnames=false:

st2 = summarystats(chain; split_varnames=false)

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (2) ── VarName                                                    │
│  Float64                          f                                          │
│  Union{Missing, Vector{Float64}}  v                                          │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param       mean        std     mcse  ess_bulk  ess_tail     rhat  …       │
│       f    -0.0267     0.9731   0.1061   92.5113   64.4285   0.9879  …       │
│       v  [0.0,0.2]  [0.9,1.1]  missing   missing   missing  missing  …       │
╰──────────────────────────────────────────────────────────────────────────────╯

Now, the summary statistics are calculated with each vector v being a single entity. The key v is still present in the result (since mean and std could be computed for it). However, notice that other statistics such as ess are no longer defined, and so return missing values. Just like before, the string s is dropped since no statistics could be computed for it.

Individual statistics

Sometimes you may only want to calculate a single statistic.

The following functions are all overloaded to accept FlexiChain objects. In all cases, they can be called with dims=:both, dims=:iter, or dims=:chain to specify the dimension over which to compute the statistic; the default is dims=:both.

All of these functions return a FlexiSummary where the :stat dimension has already been collapsed. That means that if you want to access the mean of a variable @varname(a) you don't need to further use the stat dimension:

mn = mean(chain)
# Not needed: mean_f = m[@varname(f), stat=:mean]
# Just do:
mean_f = mn[@varname(f)]

-0.02668494343427702

In fact, behind the scenes, the actual name of the statistic is retained. This means that, for example, if you convert the result to a DataFrame, the column will be named mean rather than stat:

using DataFrames;
DataFrame(mn)

3×2 DataFrame

Row	param	stat
	VarName	Float64
1	f	-0.0266849
2	v[1]	0.00136616
3	v[2]	0.201306

Statistics.mean Method

Statistics.mean(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the mean across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Statistics.mean; please see its documentation for details of supported keyword arguments.

source

Statistics.median Method

Statistics.median(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the median across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Statistics.median; please see its documentation for details of supported keyword arguments.

source

Statistics.std Method

Statistics.std(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the standard deviation across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Statistics.std; please see its documentation for details of supported keyword arguments.

source

Statistics.var Method

Statistics.var(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the variance across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Statistics.var; please see its documentation for details of supported keyword arguments.

source

Statistics.quantile Method

Statistics.quantile(
    chain::FlexiChain{TKey},
    p;
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the quantile across all iterations and chains for each key in the chain. If it cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse.

• :iter: collapse the iteration dimension only

• :chain: collapse the chain dimension only

• :both: collapse both the iteration and chain dimensions (default)

The argument p specifies the quantile to compute, and is forwarded to Statistics.quantile, along with any other keyword arguments.

source

Base.minimum Method

Base.minimum(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the minimum across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Base.minimum; please see its documentation for details of supported keyword arguments.

source

Base.maximum Method

Base.maximum(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the maximum across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Base.maximum; please see its documentation for details of supported keyword arguments.

source

Base.sum Method

Base.sum(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the sum across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Base.sum; please see its documentation for details of supported keyword arguments.

source

Base.prod Method

Base.prod(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the product across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to Base.prod; please see its documentation for details of supported keyword arguments.

source

MCMCDiagnosticTools.ess Method

MCMCDiagnosticTools.ess(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the effective sample size across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to MCMCDiagnosticTools.ess; please see its documentation for details of supported keyword arguments.

source

MCMCDiagnosticTools.rhat Method

MCMCDiagnosticTools.rhat(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the R-hat diagnostic across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to MCMCDiagnosticTools.rhat; please see its documentation for details of supported keyword arguments.

source

MCMCDiagnosticTools.mcse Method

MCMCDiagnosticTools.mcse(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the Monte Carlo standard error across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to MCMCDiagnosticTools.mcse; please see its documentation for details of supported keyword arguments.

source

StatsBase.mad Method

StatsBase.mad(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the median absolute deviation across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to StatsBase.mad; please see its documentation for details of supported keyword arguments.

source

StatsBase.geomean Method

StatsBase.geomean(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the geometric mean across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to StatsBase.geomean; please see its documentation for details of supported keyword arguments.

source

StatsBase.harmmean Method

StatsBase.harmmean(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the harmonic mean across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to StatsBase.harmmean; please see its documentation for details of supported keyword arguments.

source

StatsBase.iqr Method

StatsBase.iqr(
    chain::FlexiChain{TKey};
    dims::Symbol=:both,
    warn::Bool=true,
    split_varnames::Bool=true,
    kwargs...
) where {TKey}

Calculate the interquartile range across all iterations and chains for each key in the chain. If the statistic cannot be computed for a key, that key is skipped and a warning is issued (which can be suppressed by setting warn=false).

The dims keyword argument specifies which dimensions to collapse. The default value of :both collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The split_varnames keyword argument, if true, will first split up variables in the chain such that each key corresponds to a single scalar value. This is only supported for chains with TKey<:VarName or TKey==Symbol; for other key types this will be a no-op as long as the data already contain scalar values for every key, but will error if the data contain non-scalar values.

Other keyword arguments are forwarded to StatsBase.iqr; please see its documentation for details of supported keyword arguments.

source

There are also similar overloads for PosteriorStats.hdi and PosteriorStats.eti, which calculate highest density intervals and equal-tailed intervals, respectively. Please see the PosteriorStats integration section for details.

Custom statistics

There are two scenarios where the above are not enough: 2. you want to calculate a specific set of statistics that is not the same as what summarystats does; or

3. you want to calculate a completely custom statistic, which is not implemented above.

In both cases, you can directly use [FlexiChains.collapse] to achieve this. (But in the latter case, please do also consider opening an issue so that we can implement it!)

FlexiChains.collapse Function

FlexiChains.collapse(
    chain::FlexiChain,
    funcs::AbstractVector;
    dims::Symbol=:both,
    warn::Bool=true,
    drop_stat_dim::Bool=false,
)

Low-level function to collapse one or both dimensions of a FlexiChain by applying a list of summary functions.

The funcs argument must be a vector which contains either:

• tuples of the form (statistic_name::Symbol, func::Function); or

• just functions, in which case the statistic name is obtained from the function name.

The dims keyword argument specifies which dimensions to collapse. By default, dims is :both, which collapses both the iteration and chain dimensions. Other valid values are :iter or :chain, which respectively collapse only the iteration or chain dimension.

The functions in funcs must map a vector to a single value. For example, both Statistics.mean and Statistics.std satisfy this:

using FlexiChains: collapse
using Statistics: mean, std

collapse(chn, [mean, std]; dims=:both)

If dims=:iter or dims=:chain are selected, then the functions are automatically applied to each column or row as appropriate. No adjustment to the functions is necessary:

collapse(chn, [mean, std]; dims=:iter)
collapse(chn, [mean, std]; dims=:chain)

For dims=:both, the function is applied to all the samples stacked together as a single vector.

Sometimes, for more complicated functions like quantile, you have to pass an anonymous function (such as x -> quantile(x, 0.05) or a closure (such as Base.Fix2(quantile, 0.05)). In this case, to get a sensible statistic name, instead of just passing the function you can pass a tuple of the form (statistic_name::Symbol, func::Function).

collapse(chn, [
    mean,
    std,
    (:q5, x -> quantile(x, 0.05)),
    (:q95, x -> quantile(x, 0.95)),
])

If a statistic function errors when applied to a key, that key is skipped and a warning is issued. The warning can be suppressed by setting warn=false.

If the drop_stat_dim keyword argument is true and only one function is provided in funcs, then the resulting FlexiSummary will have the stat dimension dropped. This allows for easier indexing into the result when only one statistic is computed. It is an error to set drop_stat_dim=true when more than one function is provided.

source

As an example, suppose you have a statistic that calculates the sum of the mean and standard deviation. (This is of course quite contrived: if you have a real example, again, please do open an issue!)

We start by defining our own function:

using Statistics

function mean_std_sum(x::AbstractVector{<:Real})
    return mean(x) + std(x)
end

mean_std_sum (generic function with 1 method)

As noted in the docstring of FlexiChains.collapse, the function you provide must accept a vector argument and return a single value. Of course, it can also have other methods, but this is the one which collapse uses.

Now we can use collapse to apply this function to all variables in the chain. The second argument is a vector, which in this case will only contain our one function:

custom_stat = FlexiChains.collapse(chain, [mean_std_sum]; dims=:both)

╭─FlexiSummary (1 statistic) ──────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean_std_sum]                                                     │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64  f, v[1], v[2]                                                      │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param  mean_std_sum                                                        │
│       f        0.9464                                                        │
│    v[1]        0.8618                                                        │
│    v[2]        1.3003                                                        │
╰──────────────────────────────────────────────────────────────────────────────╯

There are two things worth mentioning, which we will note in passing here without demonstrating (since they are also covered in the docstring): 2. If there is only one function provided, you can additionally pass drop_stat_dim=true to remove the :stat dimension from the result, much like what mean(chain) et al. do.

3. The name of the statistic is inferred from the function. Sometimes this doesn't work out nicely, for example if you pass an anonymous function. In this case you can provide a tuple of (:name, function) instead of just the function.

Merging summaries

As an alternative to using collapse, you can also calculate separate summaries and then merge them.

mn = mean(chain)
sd = std(chain)

merge(mn, sd)

╭─FlexiSummary (2 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std]                                                        │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64  f, v[1], v[2]                                                      │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param     mean     std                                                     │
│       f  -0.0267  0.9731                                                     │
│    v[1]   0.0014  0.8604                                                     │
│    v[2]   0.2013  1.0990                                                     │
╰──────────────────────────────────────────────────────────────────────────────╯

Base.merge Method

Base.merge(
    s1::FlexiSummary{TKey1},
    s2::FlexiSummary{TKey2},
) where {TKey1,TKey2}

Merge two FlexiSummarys along both the key and statistic dimensions.

The two summaries must have the same iteration and chain indices.

The merged result contains the union of keys and the union of stat names from both summaries. For each (key, stat) pair, the value from s2 takes priority over s1. Pairs that exist in neither summary are filled with missing.

If the key types are different, the resulting FlexiSummary will have a promoted key type, and a warning will be issued.

source

Density intervals

It is worth just noting briefly here that FlexiChains has a PosteriorStats extension which allows you to calculate highest density intervals and equal-tailed intervals. Please see the PosteriorStats integration section for details.

MCMC diagnostics

There are a number of functions in MCMCDiagnosticTools.jl which do not nicely fit the notion of 'collapse over one or more dimensions'. For example, quantities like the Gelman diagnostic are calculated over the entire dataset.

For these functions, you can simply pass a FlexiChain to them. For example:

gelmandiag(chain)

╭─FlexiSummary (2 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [psrf, psrfci]                                                     │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64  f, v[1], v[2]                                                      │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
│                                                                              │
│ Summary                                                                      │
│   param    psrf  psrfci                                                      │
│       f  0.9879  0.9990                                                      │
│    v[1]  0.9986  1.0403                                                      │
│    v[2]  1.0394  1.1790                                                      │
╰──────────────────────────────────────────────────────────────────────────────╯

The available functions are as follows (and FlexiChains reexports all of these too, for convenience).

MCMCDiagnosticTools.gelmandiag Function

MCMCDiagnosticTools.gelmandiag(
    chain::FlexiChain{TKey};
    warn::Bool=true,
    kwargs...
) where {TKey}

Compute the Gelman–Rubin–Brooks diagnostic (Potential Scale Reduction Factor, PSRF) for each parameter in the chain.

Returns a FlexiSummary with two statistics per parameter: :psrf (the point estimate) and :psrfci (the upper confidence limit).

The FlexiChain must contain at least 2 chains. Array-valued parameters are automatically split into their constituent scalars. Non-Real-valued keys are skipped with a warning (which can be suppressed via warn=false).

Other keyword arguments are forwarded to MCMCDiagnosticTools.gelmandiag.

source

MCMCDiagnosticTools.gelmandiag_multivariate Function

MCMCDiagnosticTools.gelmandiag_multivariate(
    chain::FlexiChain{TKey};
    warn::Bool=true,
    kwargs...
) where {TKey}

Compute the multivariate Gelman–Rubin–Brooks diagnostic for the chain. This requires at least 2 parameters and at least 2 chains.

Returns a NamedTuple with two fields:

• summary: a FlexiSummary with :psrf and :psrfci statistics (same as what gelmandiag returns)

• psrf_multivariate: the multivariate potential scale reduction factor (Float64)

Array-valued parameters are automatically split into their constituent scalars. Non-Real-valued keys are skipped with a warning (which can be suppressed via warn=false).

Other keyword arguments are forwarded to MCMCDiagnosticTools.gelmandiag_multivariate.

source

MCMCDiagnosticTools.discretediag Function

MCMCDiagnosticTools.discretediag(
    chain::FlexiChain{TKey};
    warn::Bool=true,
    kwargs...
) where {TKey}

Compute the discrete diagnostic for each parameter in the chain. This diagnostic is designed for discrete (categorical) MCMC samples.

Returns a NamedTuple with two fields:

• between: a FlexiSummary with between-chain diagnostics (:stat, :df, :pvalue) per parameter

• within: a FlexiSummary with within-chain diagnostics (:stat, :df, :pvalue) per parameter and per chain

The FlexiChain must contain at least 2 chains. Array-valued parameters are automatically split into their constituent scalars. Non-Integer-valued keys are skipped with a warning (which can be suppressed via warn=false).

Other keyword arguments are forwarded to MCMCDiagnosticTools.discretediag.

source

MCMCDiagnosticTools.bfmi Function

MCMCDiagnosticTools.bfmi(
    chain::FlexiChain,
    energy_key
)

Calculate the Bayesian fraction of missing information (BFMI) from the given chain, using the specified energy_key to identify the key in the chain that corresponds to the Hamiltonian energy. Returns a DimVector of BFMI values, one per chain (note that even if there is only one chain, the result will still be a vector of length 1).

For chains sampled with Turing.jl's HMC/NUTS, the energy key is :hamiltonian_energy.

source

Flattening summaries

Finally, FlexiSummary objects can be converted to DimArrays or DataFrames, just like FlexiChains. This is useful for various analyses and visualisations.

st = summarystats(chain)
DimArray(st)

┌ 9×3 DimArray{Float64, 2} ┐
├──────────────────────────┴──────────────────────── dims ┐
  ↓ stat Categorical{Symbol} [:mean, …, :q95] Unordered,
  → param Categorical{VarName} [f, …, v[2]] Unordered
└─────────────────────────────────────────────────────────┘
 ↓ →           f           v[1]         v[2]
  :mean      -0.0266849   0.00136616   0.201306
  :std        0.973082    0.860395     1.09896
  :mcse       0.106082    0.101039     0.135994
  :ess_bulk  92.5113     71.1998      68.9064
  :ess_tail  64.4285     37.5486      23.2807
  :rhat       0.98788     0.994136     1.00172
  :q5        -1.65059    -1.43801     -1.66004
  :q50        0.133813   -0.0701984    0.256264
  :q95        1.40124     1.24931      2.13617

using DataFrames

DataFrame(st)

3×10 DataFrame

Row	param	mean	std	mcse	ess_bulk	ess_tail	rhat	q5	q50	q95
	VarName	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	f	-0.0266849	0.973082	0.106082	92.5113	64.4285	0.98788	-1.65059	0.133813	1.40124
2	v[1]	0.00136616	0.860395	0.101039	71.1998	37.5486	0.994136	-1.43801	-0.0701984	1.24931
3	v[2]	0.201306	1.09896	0.135994	68.9064	23.2807	1.00172	-1.66004	0.256264	2.13617