Implementation of the No U-Turn Sampler for MCMC.

Please read the documentation. For the impatient: you probably want to

1. define a log density problem (eg for Bayesian inference) using the LogDensityProblems

package, then

1. use it with mcmc_with_warmup.

Default maximum depth for trees.

Docstring for initial warmup arguments.

Shared docstring part for the MCMC API.

Maximum number of iterations next_direction supports.

Sentinel value for reaching maximum depth.

Significant digits to display for reporting.

Internal type implementing random directions. Draw a new value with rand, see next_direction.

Serves two purposes: a fixed value of Directions is useful for unit testing, and drawing a single bit flag collection economizes on the RNG cost.

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

Current state of adaptation for ϵ.

struct DynamicHMCError <: Exception

The error type used by this package. Debug information should be printed without truncation, with full precision.

• message

• debug_information

abstract type EuclideanKineticEnergy <: DynamicHMC.KineticEnergy

Euclidean kinetic energies (position independent).

struct EvaluatedLogDensity{T, S}

A log density evaluated at position q. The log densities and gradient are saved, so that they are not calculated twice for every leapfrog step (both as start- and endpoints).

Because of caching, a EvaluatedLogDensity should only be used with a specific Hamiltonian, preferably constructed with the evaluate_ℓ constructor.

In composite types and arguments, Q is usually used for this type.

Adaptation with fixed stepsize. Leaves ϵ unchanged.

GaussianKineticEnergy(N)
GaussianKineticEnergy(N, m⁻¹)

Gaussian kinetic energy with a diagonal inverse covariance matrix M⁻¹=m⁻¹*I.

struct GaussianKineticEnergy{T<:(AbstractMatrix), S<:(AbstractMatrix)} <: 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.

GaussianKineticEnergy(M⁻¹)

Gaussian kinetic energy with the given inverse covariance matrix M⁻¹.

GaussianKineticEnergy(M⁻¹)

Gaussian kinetic energy with the given inverse covariance matrix M⁻¹.

Statistics for the identification of turning points. See Betancourt (2017, appendix), and subsequent discussion of improvements at https://discourse.mc-stan.org/t/nuts-misses-u-turns-runs-in-circles-until-max-treedepth/9727/.

Momenta p₋ and p₊ are kept so that they can be added to ρ when combining turn statistics.

Turn detection is always done by combine_turn_statistics, which returns nothing in case of turning. A GeneralizedTurnStatistic should always correspond to a trajectory that is not turning (or a leaf node, where the concept does not apply).

struct InitialStepsizeSearch

Parameters for the search algorithm for the initial stepsize.

The algorithm finds an initial stepsize $ϵ$ so that the local log acceptance ratio $A(ϵ)$ is near params.log_threshold.

• initial_ϵ: The stepsize where the search is started.

• log_threshold: Log of the threshold that needs to be crossed.

• maxiter_crossing: Maximum number of iterations for crossing the threshold.

!!! NOTE

The algorithm is from Hoffman and Gelman (2014), default threshold modified to `0.8` following later practice in Stan.

Information about an invalid (sub)tree, using positions relative to the starting node.

1. When left < right, this tree was turning.

2. When left == right, this is a divergent node.

3. left == 1 && right == 0 is used as a sentinel value for reaching maximum depth without

encountering any invalid trees (see REACHED_MAX_DEPTH. All other left > right values are disallowed.

abstract type KineticEnergy

Kinetic energy specifications. Implements the methods

• Base.size

• kinetic_energy

• calculate_p♯

• ∇kinetic_energy

• rand_p

For all subtypes, it is implicitly assumed that kinetic energy is symmetric in the momentum p,

kinetic_energy(κ, p, q) == kinetic_energy(κ, .-p, q)

When the above is violated, the consequences are undefined.

mutable struct LogMCMCReport{T}

A composite type for tracking the state for which the last log message was emitted, for MCMC reporting with a given total number of steps (see make_mcmc_reporter.

Fields

• log_progress_report: The progress report sink.

• total_steps: Total steps for this stage.

• last_reported_step: Index of the last reported step.

• last_reported_time_ns: The last time a report was logged (determined using time_ns).

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.

• step_interval: Always report progress past step_interval of the last report.

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

struct MCMCSteps{TR, TA, TH, TE}

A composite type for performing MCMC stepwise after warmup.

The type is not part of the API, see mcmc_steps and mcmc_next_step.

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.

struct NoProgressReport

A placeholder type for not reporting any information.

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.

struct ProgressMeterReport

Report progress via a progress bar, using ProgressMeter.jl.

Example usage:

julia> ProgressMeterReport()

struct SamplingLogDensity{R, L, O, S}

A log density bundled with an RNG and options for sampling. Contains the parts of the problem which are not changed during warmup.

Fields

• rng: Random number generator.

• ℓ: Log density.

• algorithm: Algorithm used for sampling, also contains the relevant parameters that are not affected by adaptation. See eg NUTS.

• reporter: Reporting warmup information and chain progress.

Representation of a trajectory, ie a Hamiltonian with a discrete integrator that also checks for divergence.

struct TreeStatisticsNUTS

Diagnostic information for a single tree built with the No-U-turn sampler.

Fields

Accessing fields directly is part of the API.

• π: Log density (negative energy).

• depth: Depth of the tree.

• termination: Reason for termination. See InvalidTree and REACHED_MAX_DEPTH.

• acceptance_rate: Acceptance rate statistic.

• steps: Number of leapfrog steps evaluated.

• directions: Directions for tree doubling (useful for debugging).

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:

• posterior_matrix, a matrix of position vectors, indexes by [parameter_index, draw_index]

• 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.

struct WarmupState{TQ<:DynamicHMC.EvaluatedLogDensity, Tκ<:DynamicHMC.KineticEnergy, Tϵ<:Union{Nothing, Real}}

Representation of a warmup state. Not part of the API.

Fields

• Q: phasepoint

• κ: kinetic energy

• ϵ: stepsize

_doubling_warmup_stages(
    M,
    stepsize_adaptation,
    middle_steps,
    doubling_stages
)

Helper function for constructing the “middle” doubling warmup stages in default_warmup_stages.

_empty_posterior_matrix(Q, N)

Create an empty posterior matrix, based on Q (a logdensity evaluated at a position).

_error(message; kwargs...)

Throw a DynamicHMCError with given message, keyword arguments used for debug information.

_is_turning(p♯₋, p♯₊, ρ)

Internal test for turning. See Betancourt (2017, appendix).

_log_meta(chain_id, meta)

Assemble log message metadata.

Currently, it adds chain_id iff it is not nothing.

_warmup(sampling_logdensity, stages, initial_warmup_state)

Helper function for implementing warmup.

acceptance_rate(A)

Return the acceptance rate (a Real betwen 0 and 1).

adapt_stepsize(parameters, A, a)

Update the adaptation A of log stepsize logϵ with average Metropolis acceptance rate a over the whole visited trajectory, using the dual averaging algorithm of Gelman and Hoffman (2014, Algorithm 6). Return the new adaptation state.

result, v = adjacent_tree(rng, trajectory, z, i, depth, is_forward)

Traverse the tree of given depth adjacent to point z in trajectory.

is_forward specifies the direction, rng is used for random numbers in combine_proposals. i is an integer position relative to the initial node (0).

The first value is either

1. an InvalidTree, indicating the first divergent node or turning subtree that was

encounteted and invalidated this tree.

1. a tuple of `(ζ, ω, τ, z′, i′), with

   • ζ: the proposal from the tree.

   • ω: the log weight of the subtree that corresponds to the proposal

   • τ: turn statistics

   • z′: the last node of the tree

   • i′: the position of the last node relative to the initial node.

The second value is always the visited node statistic.

biased_progressive_logprob2(bias, ω₁, ω₂)
biased_progressive_logprob2(bias, ω₁, ω₂, ω)

Given (relative) log probabilities ω₁ and ω₂, return the log probabiliy of drawing a sample from the second (logprob2).

When bias, biases towards the second argument, introducing anti-correlations.

calculate_logprob2(trajectory, is_doubling::Bool, ω₁, ω₂, ω)

Calculate the log probability if selecting the subtree corresponding to ω₂. Being the log of a probability, it is always ≤ 0, but implementations are allowed to return and accept values > 0 and treat them as 0.

When is_doubling, the tree corresponding to ω₂ was obtained from a doubling step (this can be relevant eg for biased progressive sampling).

The value ω = logaddexp(ω₁, ω₂) is provided for avoiding redundant calculations.

See biased_progressive_logprob2 for an implementation.

calculate_p♯(κ, p)
calculate_p♯(κ, p, q)

Return $p♯ = M⁻¹⋅p$, used for turn diagnostics.

combine_proposals(rng, trajectory, ζ₁, ζ₂, logprob2::Real, is_forward::Bool)

Combine two proposals ζ₁, ζ₂ on trajectory, with log probability logprob2 for selecting ζ₂.

ζ₁ is before ζ₂ iff is_forward.

combine_turn_statistics(trajectory, τ₁, τ₂)

Combine turn statistics on trajectory. Implementation can assume that the trees that correspond to the turn statistics have the same ordering.

When

τ = combine_turn_statistics(trajectory, τ₁, τ₂)
is_turning(trajectory, τ)

the combined turn statistic τ is guaranteed not to escape the caller, so it can eg change type.

combine_turn_statistics_in_direction(
    trajectory,
    τ₁,
    τ₂,
    is_forward
)

Combine turn statistics with the given direction. When is_forward, τ₁ is before τ₂, otherwise after.

Internal helper function.

combine_visited_statistics(trajectory, v₁, v₂)

Combine visited node statistics for adjacent trees trajectory. Implementation should be invariant to the ordering of v₁ and v₂ (ie the operation is commutative).

current_ϵ(A)
current_ϵ(A, tuning)

Return the stepsize ϵ for the next HMC step while adapting.

default_reporter(; kwargs...)

Return a default reporter, taking the environment into account. Keyword arguments are passed to constructors when applicable.

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.

evaluate_ℓ(ℓ, q; strict)

Evaluate log density and gradient and save with the position. Preferred interface for creating EvaluatedLogDensity instances.

Non-finite elements in q always throw an error.

Non-finite and not -Inf elements in the log density throw an error if strict, otherwise replace the log density with -Inf.

Non-finite elements in the gradient throw an error if strict, otherwise replace the log density with -Inf.

final_ϵ(A)
final_ϵ(A, tuning)

Return the final stepsize ϵ after adaptation.

find_initial_stepsize(parameters, A)

Find an initial stepsize that matches the conditions of parameters (see InitialStepsizeSearch).

A is the local log acceptance ratio (uncapped). Cf local_log_acceptance_ratio.

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.

initial_adaptation_state(_, ϵ)

Return an initial adaptation state for the adaptation method and a stepsize ϵ.

initialize_warmup_state(rng, ℓ; q, κ, ϵ)

Create an initial warmup state from a random position.

Keyword arguments

• 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.

is_turning(trajectory, τ)

Test if the turn statistics τ indicate that the corresponding tree is turning.

Will only be called on nontrivial trees (at least two nodes).

kinetic_energy(κ, p)
kinetic_energy(κ, p, q)

Return kinetic energy κ, at momentum p.

ζωτ_or_nothing, v = leaf(trajectory, z, is_initial)

Information for a tree made of a single node. When is_initial == true, this is the first node.

The first value is either

1. nothing for a divergent node,

2. a tuple containing the proposal ζ, the log weight (probability) of the node ω, the

turn statistics τ (never tested with is_turning for leaf nodes).

The second value is the visited node information.

leaf_acceptance_statistic(Δ, is_initial)

Acceptance statistic for a leaf. The initial leaf is considered not to be visited.

leapfrog(H, z, ϵ)

Take a leapfrog step of length ϵ from z along the Hamiltonian H.

Return the new phase point.

The leapfrog algorithm uses the gradient of the next position to evolve the momentum. If this is not finite, the momentum won't be either, logdensity above will catch this and return an -Inf, making the point divergent.

local_log_acceptance_ratio(H, z)

Return a function of the stepsize ($ϵ$) that calculates the local log acceptance ratio for a single leapfrog step around z along the Hamiltonian H. Formally, let

A(ϵ) = logdensity(H, leapfrog(H, z, ϵ)) - logdensity(H, z)

Note that the ratio is not capped by 0, so it is not a valid (log) probability per se.

logdensity(H, z)

Log density for Hamiltonian H at point z.

If ℓ(q) == -Inf (rejected), skips the kinetic energy calculation.

Non-finite values (incl NaN, Inf) are automatically converted to -Inf. This can happen if

1. the log density is not a finite value,

2. the kinetic energy is not a finite value (which usually happens when NaN or Inf got

mixed in the leapfrog step, leading to an invalid position).

make_mcmc_reporter(
    reporter,
    total_steps;
    currently_warmup,
    meta...
)

Return a reporter which can be used for progress reports with a known number of total_steps. May return the same reporter, or a related object. Will display meta as key-value pairs.

Arguments:

• reporter::NoProgressReport: the original reporter
• total_steps: total number of steps

Keyword arguments:

• currently_warmup::Bool: true if we are currently doing warmup; false if we are currently doing MCMC
• meta: key-value pairs that will be displayed by the reporter

mcmc(sampling_logdensity, N, warmup_state)

Markov Chain Monte Carlo for sampling_logdensity, with the adapted warmup_state.

Return a NamedTuple of

• posterior_matrix, a matrix of position vectors, indexes by [parameter_index, draw_index]

• tree_statistics, a vector of tree statistics for each sample

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 posterior_matrix and tree_statistics, see mcmc_with_warmup.

• sampling_logdensity, which contains information that is invariant to warmup

Arguments

• rng: the random number generator, eg Random.default_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.

mcmc_next_step(mcmc_steps, Q)

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

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

   • rng, the random number generator,
   • algorithm, see mcmc_with_warmup,
   • κ, the (adapted) metric,
   • ℓ, the log density callable (see mcmc_with_warmup,
   • ϵ, the stepsize.

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

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

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

• posterior_matrix, a matrix of position vectors, indexes by [parameter_index, draw_index]

• tree_statistics, a vector of tree statistics for each sample

• κ and ϵ, the adapted metric and stepsize.

Arguments

• rng: the random number generator, eg Random.default_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))

move(trajectory, z, is_forward)

Move along the trajectory in the specified direction. Return the new position.

next_direction(directions)

Return the next direction flag and the new state of directions. Results are undefined for more than MAX_DIRECTIONS_DEPTH updates.

pool_posterior_matrices(results)

Given a vector of results, each containing a property posterior_matrix (eg obtained from mcmc_with_warmup with the same sample length), return a lazy view as an array indexed by [parameter_index, pooled_draw_index].

This is useful for posterior analysis after diagnostics (see eg Base.eachcol).

rand_bool_logprob(rng, logprob)

Random boolean which is true with the given probability exp(logprob), which can be ≥ 1 in which case no random value is drawn.

rand_p(rng, κ)
rand_p(rng, κ, q)

Generate a random momentum from a kinetic energy at position q.

random_position(rng, N)

Helper function to create random starting positions in the [-2,2]ⁿ box.

regularize_M⁻¹(Σ, λ)

Adjust the inverse metric estimated from the sample, using an ad-hoc shrinkage method.

report(reporter, step; meta...)

Report to the given reporter.

The second argument can be

1. a string, which is displayed as is (this is supported by all reporters).

2. or a step in an MCMC chain with a known number of steps for progress reporters (see

make_mcmc_reporter).

meta arguments are key-value pairs.

In this context, a step is a NUTS transition, not a leapfrog step.

sample_M⁻¹(_, posterior_matrix)

Estimate the inverse metric from the chain.

In most cases, this should be regularized, see regularize_M⁻¹.

sample_trajectory(rng, trajectory, z, max_depth, directions)

Sample a trajectory starting at z, up to max_depth. directions determines the tree expansion directions.

Return the following values

• ζ: proposal from the tree

• v: visited node statistics

• termination: an InvalidTree (this includes the last doubling step turning, which is technically a valid tree) or REACHED_MAX_DEPTH when all subtrees were valid and no turning happens.

• depth: the depth of the tree that was sampled from. Doubling steps that lead to an invalid adjacent tree do not contribute to depth.

sample_tree(rng, algorithm, H, Q, ϵ; p, directions)

No-U-turn Hamiltonian Monte Carlo transition, using Hamiltonian H, starting at evaluated log density position Q, using stepsize ϵ. Parameters of algorithm govern the details of tree construction.

Return two values, the new evaluated log density position, and tree statistics.

stack_posterior_matrices(results)

Given a vector of results, each containing a property posterior_matrix (eg obtained from mcmc_with_warmup with the same sample length), return a lazy view as an array indexed by [draw_index, chain_index, parameter_index].

This is useful as an input for eg MCMCDiagnosticTools.ess_rhat.

warmup(sampling_logdensity, warmup_stage, warmup_state)

Return the results and the next warmup state after warming up/adapting according to warmup_stage, starting from warmup_state.

Use nothing for a no-op.

warmup(sampling_logdensity, tuning, warmup_state)

Perform a warmup on a given sampling_logdensity, using the specified tuning, starting from warmup_state.

Return two values. The first is either nothing, or a NamedTuple of

• posterior_matrix, a matrix of position vectors, indexes by [parameter_index, draw_index]

• tree_statistics, a vector of tree statistics for each sample

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

The second is the warmup state.

∇kinetic_energy(κ, p)
∇kinetic_energy(κ, p, q)

Calculate the gradient of the logarithm of kinetic energy in momentum p.

Acceptance quantiles for TreeStatisticsSummary diagnostic summary.

struct LeapfrogTrajectory{TH, TZ, TF, Tϵ}

Implements an iterator on a leapfrog trajectory until the first non-finite log density.

Fields

• H: Hamiltonian

• z₀: Initial position

• π₀: Negative energy at initial position.

• ϵ: Stepsize (negative: move backward).

struct TreeStatisticsSummary{T<:Real, C<:NamedTuple}

Storing the output of NUTS_statistics in a structured way, for pretty printing. Currently for internal use.

Fields

• N: Sample length.

• a_mean: average_acceptance

• a_quantiles: acceptance quantiles

• termination_counts: termination counts

• depth_counts: depth counts (first element is for 0)

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

_position_information(lft, z, i)

Position information returned by leapfrog_trajectory, see documentation there. Internal function.

count_depths(tree_statistics)

Count depths in tree statistics.

count_terminations(tree_statistics)

Count termination reasons in tree_statistics, return as a NamedTuple.

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.

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.

summarize_tree_statistics(tree_statistics)

Summarize tree statistics. Mostly useful for NUTS diagnostics.