User interface

Sampling

The main entry point for sampling is

DynamicHMC.mcmc_with_warmupFunction
mcmc_with_warmup(rng, ℓ, N; initialization, warmup_stages, algorithm, reporter)

Perform MCMC with NUTS, including warmup which is not returned. Return a NamedTuple of

  • chain, a vector of positions from the posterior

  • tree_statistics, a vector of tree statistics

  • κ and ϵ, the adapted metric and stepsize.

Arguments

  • rng: the random number generator, eg Random.GLOBAL_RNG.

  • : the log density, supporting the API of the LogDensityProblems package

  • N: the number of samples for inference, after the warmup.

Keyword arguments

  • initialization: see below.

  • warmup_stages: a sequence of warmup stages. See default_warmup_stages and fixed_stepsize_warmup_stages; the latter requires an ϵ in initialization.

  • algorithm: see NUTS. It is very unlikely you need to modify this, except perhaps for the maximum depth.

  • reporter: how progress is reported. By default, verbosely for interactive sessions using the log message mechanism (see LogProgressReport, and no reporting for non-interactive sessions (see NoProgressReport).

Initialization

The initialization keyword argument should be a NamedTuple which can contain the following fields (all of them optional and provided with reasonable defaults):

  • q: initial position. Default: random (uniform [-2,2] for each coordinate).

  • κ: kinetic energy specification. Default: Gaussian with identity matrix.

  • ϵ: a scalar for initial stepsize, or nothing for heuristic finders.

Usage examples

Using a fixed stepsize:

mcmc_with_warmup(rng, ℓ, N;
                 initialization = (ϵ = 0.1, ),
                 warmup_stages = fixed_stepsize_warmup_stages())

Starting from a given position q₀ and kinetic energy scaled down (will still be adapted):

mcmc_with_warmup(rng, ℓ, N;
                 initialization = (q = q₀, κ = GaussianKineticEnergy(5, 0.1)))

Using a dense metric:

mcmc_with_warmup(rng, ℓ, N;
                 warmup_stages = default_warmup_stages(; M = Symmetric))

Disabling the initial stepsize search (provided explicitly, still adapted):

mcmc_with_warmup(rng, ℓ, N;
                 initialization = (ϵ = 1.0, ),
                 warmup_stages = default_warmup_stages(; stepsize_search = nothing))
source

Warmup

Warmup can be customized.

Default warmup sequences

A warmup sequence is just a tuple of warmup building blocks. Two commonly used sequences are predefined.

DynamicHMC.default_warmup_stagesFunction
default_warmup_stages(; stepsize_search, M, stepsize_adaptation, init_steps, middle_steps, doubling_stages, terminating_steps)

A sequence of warmup stages:

  1. select an initial stepsize using stepsize_search (default: based on a heuristic),

  2. tuning stepsize with init_steps steps

  3. tuning stepsize and covariance: first with middle_steps steps, then repeat with twice the steps doubling_stages times

  4. tuning stepsize with terminating_steps steps.

M (Diagonal, the default or Symmetric) determines the type of the metric adapted from the sample.

This is the suggested tuner of most applications.

Use nothing for stepsize_adaptation to skip the corresponding step.

source
DynamicHMC.fixed_stepsize_warmup_stagesFunction
fixed_stepsize_warmup_stages(; M, middle_steps, doubling_stages)

A sequence of warmup stages for fixed stepsize, only tuning covariance: first with middle_steps steps, then repeat with twice the steps doubling_stages times

Very similar to default_warmup_stages, but omits the warmup stages with just stepsize tuning.

source

Warmup building blocks

DynamicHMC.InitialStepsizeSearchType
struct InitialStepsizeSearch

Parameters for the search algorithm for the initial stepsize.

The algorithm finds an initial stepsize $ϵ$ so that the local acceptance ratio $A(ϵ)$ satisfies

\[a_\text{min} ≤ A(ϵ) ≤ a_\text{max}\]

This is achieved by an initial bracketing, then bisection.

  • a_min

    Lowest local acceptance rate.

  • a_max

    Highest local acceptance rate.

  • ϵ₀

    Initial stepsize.

  • C

    Scale factor for initial bracketing, > 1. Default: 2.0.

  • maxiter_crossing

    Maximum number of iterations for initial bracketing.

  • maxiter_bisect

    Maximum number of iterations for bisection.

Note

Cf. Hoffman and Gelman (2014), which does not ensure bounds for the acceptance ratio, just that it has crossed a threshold. This version seems to work better for some tricky posteriors with high curvature.

source
DynamicHMC.DualAveragingType
struct DualAveraging{T}

Parameters for the dual averaging algorithm of Gelman and Hoffman (2014, Algorithm 6).

To get reasonable defaults, initialize with DualAveraging().

Fields

  • δ

    target acceptance rate

  • γ

    regularization scale

  • κ

    relaxation exponent

  • t₀

    offset

source
DynamicHMC.TuningNUTSType
struct TuningNUTS{M, D}

Tune the step size ϵ during sampling, and the metric of the kinetic energy at the end of the block. The method for the latter is determined by the type parameter M, which can be

  1. Diagonal for diagonal metric (the default),

  2. Symmetric for a dense metric,

  3. Nothing for an unchanged metric.

Results

A NamedTuple with the following fields:

  • chain, a vector of position vectors

  • tree_statistics, a vector of tree statistics for each sample

  • ϵs, a vector of step sizes for each sample

Fields

  • N

    Number of samples.

  • stepsize_adaptation

    Dual averaging parameters.

  • λ

    Regularization factor for normalizing variance. An estimated covariance matrix Σ is rescaled by λ towards $σ²I$, where $σ²$ is the median of the diagonal. The constructor has a reasonable default.

source
DynamicHMC.GaussianKineticEnergyType
struct GaussianKineticEnergy{T<:(AbstractArray{T,2} where T), S<:(AbstractArray{T,2} where T)} <: DynamicHMC.EuclideanKineticEnergy

Gaussian kinetic energy, with $K(q,p) = p ∣ q ∼ 1/2 pᵀ⋅M⁻¹⋅p + log|M|$ (without constant), which is independent of $q$.

The inverse covariance $M⁻¹$ is stored.

Note

Making $M⁻¹$ approximate the posterior variance is a reasonable starting point.

source

Progress report

Progress reports can be explicit or silent.

DynamicHMC.LogProgressReportType
struct LogProgressReport{T}

Report progress into the Logging framework, using @info.

For the information reported, a step is a NUTS transition, not a leapfrog step.

Fields

  • chain_id

    ID of chain. Can be an arbitrary object, eg nothing. Default: nothing

  • step_interval

    Always report progress past step_interval of the last report. Default: 100

  • time_interval_s

    Always report progress past this much time (in seconds) after the last report. Default: 1000.0

source

Algorithm and parameters

You probably won't need to change these options with normal usage, except possibly increasing the maximum tree depth.

DynamicHMC.NUTSType
struct NUTS{S}

Implementation of the “No-U-Turn Sampler” of Hoffman and Gelman (2014), including subsequent developments, as described in Betancourt (2017).

Fields

  • max_depth

    Maximum tree depth.

  • min_Δ

    Threshold for negative energy relative to starting point that indicates divergence.

  • turn_statistic_configuration

    Turn statistic configuration. Currently only Val(:generalized) (the default) is supported.

source

Inspecting warmup

Note

The warmup interface below is not considered part of the exposed API, and may change with just minor version bump. It is intended for interactive use; the docstrings and the field names of results should be informative.

DynamicHMC.mcmc_keep_warmupFunction
mcmc_keep_warmup(rng, ℓ, N; initialization, warmup_stages, algorithm, reporter)

Perform MCMC with NUTS, keeping the warmup results. Returns a NamedTuple of

  • initial_warmup_state, which contains the initial warmup state

  • warmup, an iterable of NamedTuples each containing fields

    • stage: the relevant warmup stage

    • results: results returned by that warmup stage (may be nothing if not applicable, or a chain, with tree statistics, etc; see the documentation of stages)

    • warmup_state: the warmup state after the corresponding stage.

  • final_warmup_state, which contains the final adaptation after all the warmup

  • inference, which has chain and tree_statistics, see mcmc_with_warmup.

  • sampling_logdensity, which contains information that is invariant to warmup

Warning

This function is not (yet) exported because the the warmup interface may change with minor versions without being considered breaking. Recommended for interactive use.

Arguments

  • rng: the random number generator, eg Random.GLOBAL_RNG.

  • : the log density, supporting the API of the LogDensityProblems package

  • N: the number of samples for inference, after the warmup.

Keyword arguments

  • initialization: see below.

  • warmup_stages: a sequence of warmup stages. See default_warmup_stages and fixed_stepsize_warmup_stages; the latter requires an ϵ in initialization.

  • algorithm: see NUTS. It is very unlikely you need to modify this, except perhaps for the maximum depth.

  • reporter: how progress is reported. By default, verbosely for interactive sessions using the log message mechanism (see LogProgressReport, and no reporting for non-interactive sessions (see NoProgressReport).

Initialization

The initialization keyword argument should be a NamedTuple which can contain the following fields (all of them optional and provided with reasonable defaults):

  • q: initial position. Default: random (uniform [-2,2] for each coordinate).

  • κ: kinetic energy specification. Default: Gaussian with identity matrix.

  • ϵ: a scalar for initial stepsize, or nothing for heuristic finders.

source

Stepwise sampling

Note

The stepwise sampling interface below is not considered part of the exposed API, and may change with just minor version bump.

An experimental interface is available to users who wish to do MCMC one step at a time, eg until some desired criterion about effective sample size or mixing is satisfied. See the docstrings below for an example.

DynamicHMC.mcmc_stepsFunction
mcmc_steps(rng, algorithm, κ, ℓ, ϵ)

Return a value which can be used to perform MCMC stepwise, eg until some criterion is satisfied about the sample. See mcmc_next_step.

Two constructors are available:

  1. Explicitly providing

  2. Using the fields sampling_logdensity and warmup_state, eg from mcmc_keep_warmup (make sure you use eg final_warmup_state).

Example

# initialization
results = DynamicHMC.mcmc_keep_warmup(RNG, ℓ, 0; reporter = NoProgressReport())
steps = mcmc_steps(results.sampling_logdensity, results.final_warmup_state)
Q = results.final_warmup_state.Q

# a single update step
Q, tree_stats = mcmc_next_step(steps, Q)

# extract the position
Q.q
source
DynamicHMC.mcmc_next_stepFunction
mcmc_next_step(mcmc_steps, Q)

Given Q (an evaluated log density at a position), return the next Q and tree statistics.

source

Diagnostics

Note

Strictly speaking the Diagnostics submodule API is not considered part of the exposed interface, and may change with just minor version bump. It is intended for interactive use.

DynamicHMC.Diagnostics.explore_log_acceptance_ratiosFunction
explore_log_acceptance_ratios(ℓ, q, log2ϵs; rng, κ, N, ps)

From an initial position, calculate the uncapped log acceptance ratio for the given log2 stepsizes and momentums ps, N of which are generated randomly by default.

source
DynamicHMC.Diagnostics.leapfrog_trajectoryFunction
leapfrog_trajectory(ℓ, q, ϵ, positions; rng, κ, p)

Calculate a leapfrog trajectory visiting positions (specified as a UnitRange, eg -5:5) relative to the starting point q, with stepsize ϵ. positions has to contain 0, and the trajectories are only tracked up to the first non-finite log density in each direction.

Returns a vector of NamedTuples, each containin

  • z, a PhasePoint object,

  • position, the corresponding position,

  • Δ, the log density + the kinetic energy relative to position 0.

source
DynamicHMC.Diagnostics.EBFMIFunction
EBFMI(tree_statistics)

Energy Bayesian fraction of missing information. Useful for diagnosing poorly chosen kinetic energies.

Low values (≤ 0.3) are considered problematic. See Betancourt (2016).

source
DynamicHMC.PhasePointType
struct PhasePoint{T<:DynamicHMC.EvaluatedLogDensity, S}

A point in phase space, consists of a position (in the form of an evaluated log density at q) and a momentum.

source