`DynamicHMC.DynamicHMC`

— ModuleImplementation of the No U-Turn Sampler for MCMC.

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

- define a log density problem (eg for Bayesian inference) using the
`LogDensityProblems`

package, then

- use it with
`mcmc_with_warmup`

.

`DynamicHMC.DEFAULT_MAX_TREE_DEPTH`

— ConstantDefault maximum depth for trees.

`DynamicHMC.DOC_INITIAL_WARMUP_ARGS`

— ConstantDocstring for initial warmup arguments.

`DynamicHMC.DOC_MCMC_ARGS`

— ConstantShared docstring part for the MCMC API.

`DynamicHMC.MAX_DIRECTIONS_DEPTH`

— ConstantMaximum number of iterations `next_direction`

supports.

`DynamicHMC.REACHED_MAX_DEPTH`

— ConstantSentinel value for reaching maximum depth.

`DynamicHMC.REPORT_SIGDIGITS`

— ConstantSignificant digits to display for reporting.

`DynamicHMC.Directions`

— TypeInternal 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.

`DynamicHMC.DualAveraging`

— Type`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

`DynamicHMC.DualAveragingState`

— TypeCurrent state of adaptation for `ϵ`

.

`DynamicHMC.DynamicHMCError`

— Type`struct DynamicHMCError <: Exception`

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

`message`

`debug_information`

`DynamicHMC.EuclideanKineticEnergy`

— Type`abstract type EuclideanKineticEnergy <: DynamicHMC.KineticEnergy`

Euclidean kinetic energies (position independent).

`DynamicHMC.EvaluatedLogDensity`

— Type`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.

`DynamicHMC.FixedStepsize`

— TypeAdaptation with fixed stepsize. Leaves `ϵ`

unchanged.

`DynamicHMC.GaussianKineticEnergy`

— Type```
GaussianKineticEnergy(N)
GaussianKineticEnergy(N, m⁻¹)
```

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

.

`DynamicHMC.GaussianKineticEnergy`

— Type`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.

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

`DynamicHMC.GaussianKineticEnergy`

— Method```
GaussianKineticEnergy(M⁻¹)
```

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

.

`DynamicHMC.GaussianKineticEnergy`

— Method```
GaussianKineticEnergy(M⁻¹)
```

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

.

`DynamicHMC.GeneralizedTurnStatistic`

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

`DynamicHMC.InitialStepsizeSearch`

— Type`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.`

`DynamicHMC.InvalidTree`

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

When

`left < right`

, this tree was*turning*.When

`left == right`

, this is a*divergent*node.`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.

`DynamicHMC.KineticEnergy`

— Type`abstract type KineticEnergy`

Kinetic energy specifications. Implements the methods

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.

`DynamicHMC.LogMCMCReport`

— Type`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`

).

`DynamicHMC.LogProgressReport`

— Type`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.

`DynamicHMC.MCMCSteps`

— Type`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`

.

`DynamicHMC.NUTS`

— Type`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.

`DynamicHMC.NoProgressReport`

— Type`struct NoProgressReport`

A placeholder type for not reporting any information.

`DynamicHMC.PhasePoint`

— Type`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.

`DynamicHMC.ProgressMeterReport`

— Type`struct ProgressMeterReport`

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

.

Example usage:

`julia> ProgressMeterReport()`

`DynamicHMC.SamplingLogDensity`

— Type`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.

`DynamicHMC.TrajectoryNUTS`

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

`DynamicHMC.TreeStatisticsNUTS`

— Type`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).

`DynamicHMC.TuningNUTS`

— Type`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

`Diagonal`

for diagonal metric (the default),`Symmetric`

for a dense metric,`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.

`DynamicHMC.WarmupState`

— Type`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

`DynamicHMC._doubling_warmup_stages`

— Method```
_doubling_warmup_stages(
M,
stepsize_adaptation,
middle_steps,
doubling_stages
)
```

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

.

`DynamicHMC._empty_posterior_matrix`

— Method```
_empty_posterior_matrix(Q, N)
```

Create an empty posterior matrix, based on `Q`

(a logdensity evaluated at a position).

`DynamicHMC._error`

— Method```
_error(message; kwargs...)
```

Throw a `DynamicHMCError`

with given message, keyword arguments used for debug information.

`DynamicHMC._is_turning`

— Method```
_is_turning(p♯₋, p♯₊, ρ)
```

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

`DynamicHMC._log_meta`

— Method```
_log_meta(chain_id, meta)
```

Assemble log message metadata.

Currently, it adds `chain_id`

*iff* it is not `nothing`

.

`DynamicHMC._warmup`

— Method```
_warmup(sampling_logdensity, stages, initial_warmup_state)
```

Helper function for implementing warmup.

Changes may imply documentation updates in `mcmc_keep_warmup`

.

`DynamicHMC.acceptance_rate`

— Method```
acceptance_rate(A)
```

Return the acceptance rate (a `Real`

betwen `0`

and `1`

).

`DynamicHMC.adapt_stepsize`

— Method```
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.

`DynamicHMC.adjacent_tree`

— Method`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

- an
`InvalidTree`

, indicating the first divergent node or turning subtree that was

encounteted and invalidated this tree.

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.

`DynamicHMC.biased_progressive_logprob2`

— Function```
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.

`DynamicHMC.calculate_logprob2`

— Function`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.

`DynamicHMC.calculate_p♯`

— Function```
calculate_p♯(κ, p)
calculate_p♯(κ, p, q)
```

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

`DynamicHMC.combine_proposals`

— Function`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`

.

`DynamicHMC.combine_turn_statistics`

— Function`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.

`DynamicHMC.combine_turn_statistics_in_direction`

— Method```
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.

`DynamicHMC.combine_visited_statistics`

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

`DynamicHMC.current_ϵ`

— Function```
current_ϵ(A)
current_ϵ(A, tuning)
```

Return the stepsize `ϵ`

for the next HMC step while adapting.

`DynamicHMC.default_reporter`

— Method```
default_reporter(; kwargs...)
```

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

`DynamicHMC.default_warmup_stages`

— Method```
default_warmup_stages(
;
stepsize_search,
M,
stepsize_adaptation,
init_steps,
middle_steps,
doubling_stages,
terminating_steps
)
```

A sequence of warmup stages:

select an initial stepsize using

`stepsize_search`

(default: based on a heuristic),tuning stepsize with

`init_steps`

stepstuning stepsize and covariance: first with

`middle_steps`

steps, then repeat with twice the steps`doubling_stages`

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

`DynamicHMC.evaluate_ℓ`

— Method```
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`

.

`DynamicHMC.final_ϵ`

— Function```
final_ϵ(A)
final_ϵ(A, tuning)
```

Return the final stepsize `ϵ`

after adaptation.

`DynamicHMC.find_initial_stepsize`

— Method```
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`

.

`DynamicHMC.fixed_stepsize_warmup_stages`

— Method```
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.

`DynamicHMC.initial_adaptation_state`

— Method```
initial_adaptation_state(_, ϵ)
```

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

.

`DynamicHMC.initialize_warmup_state`

— Method```
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.

`DynamicHMC.is_turning`

— Function`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).

`DynamicHMC.kinetic_energy`

— Function```
kinetic_energy(κ, p)
kinetic_energy(κ, p, q)
```

Return kinetic energy `κ`

, at momentum `p`

.

`DynamicHMC.leaf`

— Function`ζωτ_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

`nothing`

for a divergent node,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.

`DynamicHMC.leaf_acceptance_statistic`

— Method```
leaf_acceptance_statistic(Δ, is_initial)
```

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

`DynamicHMC.leapfrog`

— Method`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.

`DynamicHMC.local_log_acceptance_ratio`

— Method```
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*.

`DynamicHMC.logdensity`

— Method```
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

the log density is not a finite value,

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

`DynamicHMC.make_mcmc_reporter`

— Method```
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

`DynamicHMC.mcmc`

— Method```
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

`DynamicHMC.mcmc_keep_warmup`

— Method```
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`NamedTuple`

s 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

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

`DynamicHMC.mcmc_next_step`

— Method```
mcmc_next_step(mcmc_steps, Q)
```

Given `Q`

(an evaluated log density at a position), return the next `Q`

and tree statistics.

`DynamicHMC.mcmc_steps`

— Method```
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:

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.

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
```

`DynamicHMC.mcmc_with_warmup`

— Method```
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))
```

`DynamicHMC.move`

— Function`move(trajectory, z, is_forward)`

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

`DynamicHMC.next_direction`

— Method```
next_direction(directions)
```

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

updates.

`DynamicHMC.pool_posterior_matrices`

— Method```
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`

).

`DynamicHMC.rand_bool_logprob`

— Method```
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.

`DynamicHMC.rand_p`

— Function```
rand_p(rng, κ)
rand_p(rng, κ, q)
```

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

.

`DynamicHMC.random_position`

— Method```
random_position(rng, N)
```

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

box.

`DynamicHMC.regularize_M⁻¹`

— Method```
regularize_M⁻¹(Σ, λ)
```

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

`DynamicHMC.report`

— Method```
report(reporter, step; meta...)
```

Report to the given `reporter`

.

The second argument can be

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

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

`meta`

arguments are key-value pairs.

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

`DynamicHMC.sample_M⁻¹`

— Method```
sample_M⁻¹(_, posterior_matrix)
```

Estimate the inverse metric from the chain.

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

.

`DynamicHMC.sample_trajectory`

— Method```
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`

.

`DynamicHMC.sample_tree`

— Method```
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.

`DynamicHMC.stack_posterior_matrices`

— Method```
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`

.

The ordering is not compatible with MCMCDiagnostictools version < 0.2.

`DynamicHMC.warmup`

— Method```
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.

`DynamicHMC.warmup`

— Method```
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.

`DynamicHMC.∇kinetic_energy`

— Function```
∇kinetic_energy(κ, p)
∇kinetic_energy(κ, p, q)
```

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

.

`DynamicHMC.Diagnostics.ACCEPTANCE_QUANTILES`

— ConstantAcceptance quantiles for `TreeStatisticsSummary`

diagnostic summary.

`DynamicHMC.Diagnostics.LeapfrogTrajectory`

— Type`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).

`DynamicHMC.Diagnostics.TreeStatisticsSummary`

— Type`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`

)

`DynamicHMC.Diagnostics.EBFMI`

— Method```
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).

`DynamicHMC.Diagnostics._position_information`

— Method```
_position_information(lft, z, i)
```

Position information returned by `leapfrog_trajectory`

, see documentation there. Internal function.

`DynamicHMC.Diagnostics.count_depths`

— Method```
count_depths(tree_statistics)
```

Count depths in tree statistics.

`DynamicHMC.Diagnostics.count_terminations`

— Method```
count_terminations(tree_statistics)
```

Count termination reasons in `tree_statistics`

, return as a `NamedTuple`

.

`DynamicHMC.Diagnostics.explore_log_acceptance_ratios`

— Method```
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.

`DynamicHMC.Diagnostics.leapfrog_trajectory`

— Method```
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 `NamedTuple`

s, each containin

`z`

, a`PhasePoint`

object,`position`

, the corresponding position,`Δ`

, the log density + the kinetic energy relative to position`0`

.

`DynamicHMC.Diagnostics.summarize_tree_statistics`

— Method```
summarize_tree_statistics(tree_statistics)
```

Summarize tree statistics. Mostly useful for NUTS diagnostics.