PosteriorStats.jl

Interval estimation

PosteriorStats.jl provides the hdi and eti functions for computing highest density intervals and equal-tailed intervals, respectively. These are overloaded for FlexiChain objects in much the same way as Statistics.mean, Statistics.std, etc. (see the Summarising page for more information on those). Here is an example:

julia

using FlexiChains, DynamicPPL, Distributions, PosteriorStats

@model function g(z)
    x ~ Normal()
    y ~ Normal(x)
    z ~ Normal(y)
end
model = g(1.0)

chn = FlexiChains._make_prior_chain(model, 100, 2)
PosteriorStats.hdi(chn; prob=0.95, split_interval=true)

╭─FlexiSummary (2 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [hdi_lower, hdi_upper]                                             │
│                                                                              │
│ Parameters (2) ── VarName                                                    │
│  Float64  x, y                                                               │
│                                                                              │
│ Extras (3)                                                                   │
│  Float64  logprior, loglikelihood, logjoint                                  │
│                                                                              │
│ Summary                                                                      │
│   param  hdi_lower  hdi_upper                                                │
│       x    -2.2877     1.8821                                                │
│       y    -2.3141     2.6136                                                │
╰──────────────────────────────────────────────────────────────────────────────╯

LOO-CV

PosteriorStats.loo, which computes Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO), is also overloaded for FlexiChain objects.

This is most easily used in conjunction with Turing.jl, since you can use Turing's functionality to directly compute log-likelihood values.

As an example, let's use our favourite eight-schools model.

julia

using FlexiChains, DynamicPPL, Distributions, PosteriorStats, LinearAlgebra

y = [28, 8, -3, 7, -1, 1, 18, 12]
sigma = [15, 10, 16, 11, 9, 11, 10, 18]
@model function eight_schools(y, sigma)
    mu ~ Normal(0, 5)
    tau ~ truncated(Cauchy(0, 5); lower=0)
    theta ~ MvNormal(fill(mu, length(y)), tau^2 * I)
    y ~ MvNormal(theta, Diagonal(sigma .^ 2))
    return nothing
end
model = eight_schools(y, sigma)

chn = FlexiChains._make_posterior_chain(model, 500, 4)

╭─FlexiChain (500 iterations, 4 chains) ───────────────────────────────────────╮
│ ↓ iter  = 1:500                                                              │
│ → chain = 1:4                                                                │
│                                                                              │
│ Parameters (3) ── VarName                                                    │
│  Float64          mu, tau                                                    │
│  Vector{Float64}  theta (8,)                                                 │
│                                                                              │
│ Extras (3)                                                                   │
│  Float64  logprior, loglikelihood, logjoint                                  │
╰──────────────────────────────────────────────────────────────────────────────╯

(Instead of sampling with Turing, we use an internal function to draw samples from the posterior, in order to avoid incurring a docs dependency on Turing.)

Now, chn is a chain which contains posterior samples. You can directly get a LOO-CV result by calling:

julia

PosteriorStats.loo(model, chn; factorize=true)

NamesAndLOOResult
├─ param_names: VarName[y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8]]
└─ loo: PSISLOOResult with estimates
 elpd  se_elpd    p  se_p 
  -31      1.5  0.8  0.32

and PSISResult with 500 draws, 4 chains, and 8 parameters
Pareto shape (k) diagnostic values:
                    Count      Min. ESS
 (-Inf, 0.5]  good  6 (75.0%)  1301
  (0.5, 0.7]  okay  2 (25.0%)  991

However, it's worth taking the scenic route to understand what's going on. To go from posterior samples to a LOO-CV result, we need to compute the pointwise log-likelihood for each observation.

To obtain a chain of log-likelihood values, we can use DynamicPPL.pointwise_loglikelihoods. Note that, because we wrote y ~ MvNormal(...), if we simply called pointwise_likelihoods we would obtain a single log-likelihood for the entire vector y. Notice the Float64 type here:

julia

DynamicPPL.pointwise_loglikelihoods(model, chn)

╭─FlexiChain (500 iterations, 4 chains) ───────────────────────────────────────╮
│ ↓ iter  = 1:500                                                              │
│ → chain = 1:4                                                                │
│                                                                              │
│ Parameters (1) ── VarName                                                    │
│  Float64  y                                                                  │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
╰──────────────────────────────────────────────────────────────────────────────╯

To get the log-likelihood for each individual y[i], we can pass the factorize=true keyword argument, which uses PartitionedDistributions.jl to calculate the conditional log-likelihood for each observation. This gives us a Vector{Float64} instead:

julia

loglik_chn = DynamicPPL.pointwise_loglikelihoods(model, chn; factorize=true)

╭─FlexiChain (500 iterations, 4 chains) ───────────────────────────────────────╮
│ ↓ iter  = 1:500                                                              │
│ → chain = 1:4                                                                │
│                                                                              │
│ Parameters (1) ── VarName                                                    │
│  Vector{Float64}  y (8,)                                                     │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
╰──────────────────────────────────────────────────────────────────────────────╯

This chain has the same structure as chn, but instead of containing posterior samples, it contains log-likelihood values for each observation in y.

julia

loglik_chn[@varname(y), iter=1:5, chain=1]

┌ 5-element DimArray{Vector{Float64}, 1} Parameter(y) ┐
├─────────────────────────────────────────────────────┴ dims ┐
  ↓ iter Sampled{Int64} 1:5 ForwardOrdered Regular Points
└────────────────────────────────────────────────────────────┘
 1  …  [-4.46891, -3.2224, -3.94225, -3.32547, -3.65321, -3.53664, -3.69061, -3.82967]
 2     [-5.10187, -3.22955, -3.97979, -3.34008, -3.34934, -3.45859, -3.84206, -3.87705]
 3     [-3.99387, -3.44055, -3.8421, -3.70315, -3.11691, -3.34709, -3.81406, -3.88125]
 4     [-4.51281, -3.2299, -3.86443, -3.32038, -3.39126, -3.34639, -3.75048, -3.83942]
 5     [-4.44078, -3.44819, -3.69721, -3.32011, -3.61941, -3.35365, -3.63907, -3.93675]

Then you can pass this to PosteriorStats.loo directly, without the model:

julia

PosteriorStats.loo(loglik_chn)

NamesAndLOOResult
├─ param_names: VarName[y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8]]
└─ loo: PSISLOOResult with estimates
 elpd  se_elpd    p  se_p 
  -31      1.5  0.8  0.32

and PSISResult with 500 draws, 4 chains, and 8 parameters
Pareto shape (k) diagnostic values:
                    Count      Min. ESS
 (-Inf, 0.5]  good  6 (75.0%)  1301
  (0.5, 0.7]  okay  2 (25.0%)  991

This returns a struct n::NamesAndLOOResult, where n.param_names contains the VarNames corresponding to each observation, and n.loo contains the actual PosteriorStats.PSISLOOResult.

Docstrings

PosteriorStats.hdi Method

julia

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

Calculate the highest density interval 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 PosteriorStats.hdi; please see its documentation for details of supported keyword arguments.

Splitting intervals

The split_interval keyword argument, if set to true, causes the output FlexiSummary to have one statistic column per interval bound, i.e., hdi_lower and hdi_upper. Note that this will only be valid if you request a single interval (i.e., method=:unimodal, which is the default in PosteriorStats.jl).

If you specify method=:multimodal, the returned FlexiSummary will have a single statistic column named hdi that contains a vector of intervals, and the split_interval argument will be ignored.

source

PosteriorStats.eti Method

julia

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

Calculate the equal-tailed interval 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).

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

Splitting intervals

The split_interval keyword argument, if set to true, causes the output FlexiSummary to have one statistic column per interval bound, i.e., eti_lower and eti_upper.

source

PosteriorStats.loo Method

julia

PosteriorStats.loo(chn::FlexiChains; kwargs...)

Calculates the leave-one-out cross-validation (LOO) statistic for a FlexiChain object. The chain must contain only log-likelihood values, and they must all be stored as parameters (not extras). Parameters that map to arrays of log-likelihood values are supported: they will be flattened before being passed to PosteriorStats.loo.

Returns a struct with the following fields:

param_names::Vector: A vector of parameter names whose log-likelihood values were used.
loo::PosteriorStats.PSISLOOResult: The return value of PosteriorStats.loo applied to the log-likelihood values extracted from the FlexiChain. This contains the statistics of interest.

Additional keyword arguments are forwarded to PosteriorStats.loo.

source

PosteriorStats.loo Method

julia

PosteriorStats.loo(
    model::DynamicPPL.Model,
    posterior_chn::FlexiChains;
    factorize::Bool=false,
    kwargs...
)

Calculates the leave-one-out cross-validation (LOO) statistic, given a model plus a posterior chain. This first uses the model and posterior chain to calculate pointwise log-likelihoods, and then uses those to calculate the LOO statistic.

Returns a struct with the following fields:

param_names::Vector: A vector of parameter names whose log-likelihood values were used.
loo::PosteriorStats.PSISLOOResult: The return value of PosteriorStats.loo applied to the log-likelihood values extracted from the FlexiChain. This contains the statistics of interest.

The factorize keyword argument is passed to DynamicPPL.pointwise_loglikelihoods. If factorize=true, factorised log-densities will be calculated for distributions that can be partitioned into blocks (e.g. MvNormal). Please see the docstring of DynamicPPL.pointwise_loglikelihoods for more details.

Additional keyword arguments are forwarded to PosteriorStats.loo.

source

PosteriorStats.jl ​

Interval estimation ​

LOO-CV ​

Docstrings ​

PosteriorStats.jl

Interval estimation

LOO-CV

Docstrings