API: Turing.Inference

Conditional SMC, an alias for PG.

ESS

Elliptical slice sampling algorithm.

Examples

julia> @model function gdemo(x)
           m ~ Normal()
           x ~ Normal(m, 0.5)
       end
gdemo (generic function with 2 methods)

julia> sample(gdemo(1.0), ESS(), 1_000) |> mean
Mean

│ Row │ parameters │ mean     │
│     │ Symbol     │ Float64  │
├─────┼────────────┼──────────┤
│ 1   │ m          │ 0.824853 │

ESSThresholdResampler(threshold, scheme = StratifiedResampler())

Resample with scheme, but only when the effective sample size drops below threshold * nparticles. This is the default for SMC and PG.

Emcee(n_walkers::Int, stretch_length=2.0)

Affine-invariant ensemble sampling algorithm.

Reference

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. (2013). emcee: The MCMC Hammer. Publications of the Astronomical Society of the Pacific, 125 (925), 306. https://doi.org/10.1086/670067

ExternalSampler{Unconstrained,S<:AbstractSampler,AD<:ADTypes.AbstractADType}

Represents a sampler that does not have a custom implementation of AbstractMCMC.step(rng, ::DynamicPPL.Model, spl).

The Unconstrained type-parameter is to indicate whether the sampler requires unconstrained space.

Fields

• sampler::AbstractMCMC.AbstractSampler: the sampler to wrap

• adtype::ADTypes.AbstractADType: the automatic differentiation (AD) backend to use

Turing.jl's interface for external samplers

If you implement a new MySampler <: AbstractSampler and want it to work with Turing.jl models, there are two options:

1. Directly implement the AbstractMCMC.step methods for DynamicPPL.Model. That is to say, implement AbstractMCMC.step(rng::Random.AbstractRNG, model::DynamicPPL.Model, sampler::MySampler; kwargs...) and related methods. This is the most powerful option and is what Turing.jl's in-house samplers do. Implementing this means that you can directly call sample(model, MySampler(), N).

2. Implement a generic AbstractMCMC.step method for AbstractMCMC.LogDensityModel (the same signature as above except that model::AbstractMCMC.LogDensityModel). This struct wraps an object that obeys the LogDensityProblems.jl interface, so your step implementation does not need to know anything about Turing.jl or DynamicPPL.jl. To use this with Turing.jl, you will need to wrap your sampler: sample(model, externalsampler(MySampler()), N).

This section describes the latter.

MySampler must implement the following methods:

• AbstractMCMC.step (the main function for taking a step in MCMC sampling; this is documented in AbstractMCMC.jl). This function must return a tuple of two elements, a 'transition' and a 'state'.

• AbstractMCMC.step_warmup (optional; documented in AbstractMCMC.jl). If your sampler adapts, implement this and Turing.jl will route the num_warmup iterations through it. Samplers that do not implement it fall back to AbstractMCMC.step, as before.

• AbstractMCMC.getparams(external_state): How to extract the parameters from the state returned by your sampler (i.e., the second return value of step). For your sampler to work with Turing.jl, this function should return a Vector of parameter values. Note that this function does not need to perform any linking or unlinking; Turing.jl will take care of this for you. You should return the parameters exactly as your sampler sees them.

• AbstractMCMC.getstats(external_state): Extract sampler statistics corresponding to this iteration from the state returned by your sampler (i.e., the second return value of step). For your sampler to work with Turing.jl, this function should return a NamedTuple. If there are no statistics to return, return NamedTuple().

  Note that getstats should not include log-probabilities as these will be recalculated by Turing automatically for you.

Notice that both of these functions take the state as input, not the transition. In other words, the transition is completely useless for the external sampler interface. This is in line with long-term plans for removing transitions from AbstractMCMC.jl and only using states.

There are a few more optional functions which you can implement to improve the integration with Turing.jl:

• AbstractMCMC.requires_unconstrained_space(::MySampler): If your sampler requires unconstrained space, you should return true. This tells Turing to perform linking on the VarInfo before evaluation, and ensures that the parameter values passed to your sampler will always be in unconstrained (Euclidean) space.

• Turing.Inference.supports_gibbs(::MySampler): If you want to disallow your sampler from a component in Turing's Gibbs sampler, you should make this evaluate to false. Note that the default is true, so you should only need to implement this in special cases.

• Turing.Inference.allow_discrete_variables(::MySampler): Return false if your sampler needs every variable to be continuous, as the gradient-based ones do. sample checks the model against this before it starts.

• Turing.Inference.gibbs_get_stats(::MyState): The statistics of the step that produced this state, as a NamedTuple, for a Gibbs chain to carry. Gibbs drops component transitions, so they cannot come from there. The wrapper implements this for you from AbstractMCMC.getstats.

• Turing.Inference.post_sample_hook(chain, ::MySampler): Anything to report once sampling has finished, such as a warning about numerical errors. Returns nothing by default.

• Turing.Inference.init_strategy(::MySampler): The DynamicPPL.AbstractInitStrategy your sampler starts from when the caller gives no initial_params. The default is InitFromPrior(); Gibbs asks each component for its own.

supports_gibbs, allow_discrete_variables, init_strategy and post_sample_hook are the whole of Turing's own interface for running a sampler on a model, together with allow_varying_dimension, which is on the list for a sampler you write yourself but is not forwarded through this wrapper, for the reason below. One that implements AbstractMCMC.step for DynamicPPL.Model (option 1 above) and overrides whichever of them do not match its defaults needs no wrapping at all; externalsampler exists only to supply the step method option 2 leaves out. Serving as a Gibbs component takes two further methods, Turing.Inference.gibbs_get_parameter_values and Turing.Inference.gibbs_update_state!!, which the wrapper implements for you – and which are also why allow_varying_dimension is not on the list above. allow_varying_dimension is not forwarded, and the wrapper implements no gibbs_update_state!! for a ReshapedBlock: the wrapped state is opaque, so there is no way to tell what within it is shaped like the block. A wrapped sampler therefore cannot own a Gibbs block whose set of variables changes, whatever it declares; one that needs to has to take option 1.

One method the wrapper cannot implement for you:

• AbstractMCMC.setparams!!(model::AbstractMCMC.LogDensityModel, state, params): required to serve as a Gibbs component. Gibbs re-conditions the model between sweeps, so this must recompute any log-density the state caches, not only write params into it. Define the three-argument form: AbstractMCMC's fallback drops model and calls a two-argument setparams!!(state, params), which cannot recompute anything, and a state that caches a log-density would then start each step from the previous conditioning's value.

Gibbs

A type representing a Gibbs sampler.

Constructors

Gibbs takes pairs of variable names and samplers. A variable name is a Symbol or a VarName, and an iterable of them assigns several variables to one component:

Gibbs(:x => NUTS(), :y => MH())
Gibbs((@varname(x), :y) => NUTS(), :z => MH())

Every variable the model reaches must belong to a component, and several components may share one. What they may not do is split a value the model stores as a unit: x[1] ~ Normal(); x[2] ~ Normal() gives each element its own key and splits cleanly, while x ~ MvNormal(...) is one key that Gibbs cannot free part of.

A component's step may reach into another block, because a variable in one block can decide something about a variable in another without sampling it. Only one of those is refused:

A component may never change the dimension of a variable belonging to another component, or whether that variable exists at all, unless the two share it.

Sharing is the remedy: put the deciding variable and what it decides in one block, or write a component that samples both.

The third row is the ordinary hierarchical case and carries no risk: m ~ Normal(0, 1) in one block and x ~ Normal(m, 1) in another changes x's distribution every sweep but not the set it is drawn from.

Examples

x decides whether z exists, so they share a block, given to a component that can sample a set of variables that changes between sweeps:

@model function f()
    x ~ Normal()
    y ~ Normal()
    if x > 0
        z ~ Normal()
    end
end

sample(f(), Gibbs(@varname(y) => MH(), (@varname(x), @varname(z)) => PG(20)), 1000)

The first row of the table is what Gibbs refuses, and the shape to recognise is a variable deciding how many of another exist:

@model function bad_dimension()
    b ~ Bernoulli(0.5)
    θ = Vector{Float64}(undef, b ? 2 : 1)
    for i in eachindex(θ)
        θ[i] ~ Normal(0, 1)
    end
    return 1.0 ~ Normal(sum(θ), 1.0)
end

Gibbs(@varname(b) => MH(), @varname(θ) => MH()) is refused, with

The variable θ[2] stopped existing during a step of the component sampling b, which does not
sample it. The component sampling θ does.

While b's component steps, θ is conditioned and cannot move, so if b flips the sweep holds a θ of the wrong length and whichever component conditions on it next is conditioning on a state that does not exist. The chain comes back biased rather than merely slow. Gibbs((@varname(b), @varname(θ)) => PG(20)) samples it correctly, because PG redraws whatever the model reaches each sweep. Declaring that is allow_varying_dimension.

@model function bad()
    b ~ Bernoulli(0.5)
    if b
        x ~ Uniform(0.0, 1.0)
    else
        x ~ Uniform(2.0, 3.0)
    end
    return 0.5 ~ Normal(x, 5.0)
end

Each component initialises the variables it samples with its own default strategy, so an HMC component starts from its InitFromUniform. A user-supplied initial_params overrides that and applies to the model as a whole; it cannot yet be set per component.

Fields

• varnames::NTuple{N, AbstractVector{<:AbstractPPL.VarName}} where N: varnames representing variables for each sampler

• samplers::NTuple{N, Any} where N: samplers for each entry in varnames

GibbsConditional(get_cond_dists)

A Gibbs component sampler that samples variables according to user-provided analytical conditional posterior distributions.

When using Gibbs sampling, sometimes one may know the analytical form of the posterior for a given variable, given the conditioned values of the other variables. In such cases one can use GibbsConditional as a component sampler to to sample from these known conditionals directly, avoiding any MCMC methods. One does so with

sampler = Gibbs(
    (@varname(var1), @varname(var2)) => GibbsConditional(get_cond_dists),
    other samplers go here...
)

Here get_cond_dists(vnt::VarNamedTuple) should be a function that takes a VarNamedTuple that contains the values of all other variables (apart from var1 and var2), and returns the conditional posterior distributions for var1 and var2.

VarNamedTuples behave very similarly to Dict{VarName,Any}s, but are more efficient and more general: you can obtain values simply by using, e.g. vnt[@varname(var3)]. See https://turinglang.org/docs/usage/varnamedtuple/ for more details on VarNamedTuples.

You may, of course, have any number of variables being sampled as a block in this manner, we only use two as an example.

The return value of get_cond_dists(vnt) should be one of the following:

• A single Distribution, if only one variable is being sampled.
• A VarNamedTuple of Distributions, which represents a mapping from variable names to their conditional posteriors. Please see the documentation linked above for information on how to construct VarNamedTuples.

For convenience, we also allow the following return values (which are internally converted into a VarNamedTuple):

• A NamedTuple of Distributions, which is like the AbstractDict case but can be used if all the variable names are single Symbols, e.g.: (; var1=dist1, var2=dist2).
• An AbstractDict{<:VarName,<:Distribution} that maps the variables being sampled to their conditional posteriors E.g. Dict(@varname(var1) => dist1, @varname(var2) => dist2).

Note that the AbstractDict case is likely to incur a performance penalty; we recommend using VarNamedTuples directly.

Examples

using Turing

# Define a model
@model function inverse_gdemo(x)
    precision ~ Gamma(2, inv(3))
    std = sqrt(1 / precision)
    m ~ Normal(0, std)
    for i in eachindex(x)
        x[i] ~ Normal(m, std)
    end
end

# Define analytical conditionals. See
# https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution
function cond_precision(vnt)
    a = 2.0
    b = 3.0
    m = vnt[@varname(m)]
    x = vnt[@varname(x)]
    n = length(x)
    a_new = a + (n + 1) / 2
    b_new = b + sum(abs2, x .- m) / 2 + m^2 / 2
    return Gamma(a_new, 1 / b_new)
end

function cond_m(vnt)
    precision = vnt[@varname(precision)]
    x = vnt[@varname(x)]
    n = length(x)
    m_mean = sum(x) / (n + 1)
    m_var = 1 / (precision * (n + 1))
    return Normal(m_mean, sqrt(m_var))
end

# Sample using GibbsConditional
model = inverse_gdemo([1.0, 2.0, 3.0])
chain = sample(model, Gibbs(
    :precision => GibbsConditional(cond_precision),
    :m => GibbsConditional(cond_m)
), 1000)

GibbsInitStrategy(varnames, strategies)

Initialise each variable with the strategy of the component sampler that samples it, choosing by ownership; a variable no component claims falls back to the prior.

This produces the single initial draw the first sweep starts from, and nothing else. Each component's own initialisation does not go through it: gibbs_initialstep_recursive asks init_strategy(sampler) directly, because ownership is ambiguous where components overlap. So an HMC component gets its InitFromUniform starting point from there, not from here.

HMC(ϵ::Float64, n_leapfrog::Int; adtype::ADTypes.AbstractADType = AutoForwardDiff())

Hamiltonian Monte Carlo sampler with static trajectory.

Arguments

• ϵ: The leapfrog step size to use.
• n_leapfrog: The number of leapfrog steps to use.
• adtype: The automatic differentiation (AD) backend. If not specified, ForwardDiff is used, with its chunksize automatically determined.

Usage

HMC(0.05, 10)

Tips

If you are receiving gradient errors when using HMC, try reducing the leapfrog step size ϵ, e.g.

# Original step size
sample(gdemo([1.5, 2]), HMC(0.1, 10), 1000)

# Reduced step size
sample(gdemo([1.5, 2]), HMC(0.01, 10), 1000)

HMCDA(
    n_adapts::Int, δ::Float64, λ::Float64; ϵ::Float64 = 0.0;
    adtype::ADTypes.AbstractADType = AutoForwardDiff(),
)

Hamiltonian Monte Carlo sampler with Dual Averaging algorithm.

Usage

HMCDA(200, 0.65, 0.3)

Arguments

• n_adapts: Numbers of samples to use for adaptation.
• δ: Target acceptance rate. 65% is often recommended.
• λ: Target leapfrog length.
• ϵ: Initial step size; 0 means automatically search by Turing.
• adtype: The automatic differentiation (AD) backend. If not specified, ForwardDiff is used, with its chunksize automatically determined.

Reference

For more information, please view the following paper (arXiv link):

Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning Research 15, no. 1 (2014): 1593-1623.

MH()

Construct a Metropolis-Hastings sampler that draws proposals from the model prior.

MH(cov_matrix)

Construct a Metropolis-Hastings sampler that performs random-walk sampling in linked space, with proposals drawn from a multivariate normal distribution with the given covariance matrix. Its dimension must match the complete linked parameter vector.

MHState(vi, accepted)

MH's state: the varinfo it returns, and whether its last step accepted.

The flag rides on the transition as well, but Gibbs steps its components with discard_sample=true and keeps only the state, so a statistic reachable only from the transition never reaches the chain.

Multinomial resampling: n independent draws from the categorical over weights.

NUTS(n_adapts::Int, δ::Float64; max_depth::Int=10, Δ_max::Float64=1000.0, init_ϵ::Float64=0.0, adtype::ADTypes.AbstractADType=AutoForwardDiff())

No-U-Turn Sampler (NUTS) sampler.

Usage

NUTS()            # Use default NUTS configuration.
NUTS(1000, 0.65)  # Use 1000 adaption steps, and target accept ratio 0.65.

Arguments

• n_adapts::Int: The number of samples to use with adaptation.
• δ::Float64: Target acceptance rate for dual averaging.
• max_depth::Int: Maximum doubling tree depth.
• Δ_max::Float64: Maximum divergence during doubling tree.
• init_ϵ::Float64: Initial step size; 0 means automatically searching using a heuristic procedure.
• adtype::ADTypes.AbstractADType: The automatic differentiation (AD) backend. If not specified, ForwardDiff is used, with its chunksize automatically determined.

struct PG{R<:Turing.Inference.AbstractResampler} <: Turing.Inference.ParticleInference

Particle Gibbs (conditional SMC) sampler.

Fields

• nparticles::Int64: number of particles

• resampler::Turing.Inference.AbstractResampler: resampling scheme

• multithreaded::Bool: reweight the particles across threads within each sweep

PG(n, [resampler = ESSThresholdResampler(0.5)]; multithreaded = false)
PG(n, [scheme = StratifiedResampler(), ]threshold; multithreaded = false)

Particle Gibbs sampler with n particles. By default resampling is triggered whenever the effective sample size drops below half the number of particles. The selected scheme applies to the unconditional first sweep only; conditional sweeps draw their ancestors from the categorical over the weights, for the reason given in the resampling-schemes section of this file.

Set multithreaded = true to evaluate the particles across threads within each sweep; results are unchanged (start Julia with multiple threads, e.g. julia -t auto, for this to have effect). Threads are the only option here: a suspended particle is a live Libtask task and cannot be serialised, so a single sweep cannot be spread across processes. Passing MCMCThreads() or MCMCDistributed() to sample parallelises whole chains instead, which is a separate axis and composes with this one.

PGState(trajectory)

Particle Gibbs sampler state: the retained trajectory's raw values, which the next sweep's reference particle reuses. Plain data, because sampler state has to survive save_state=true and MCMCDistributed(), whereas the Particle it is read off owns a live Libtask.TapedTask that cannot be serialised. Nothing else needs carrying over: the reference consumes no randomness of its own, and every other particle is seeded from the sampler's rng.

Particle(model, rng)
Particle(model, rng, reference::DynamicPPL.VarNamedTuple)

A single particle: a suspended model execution together with its varinfo, its own rng, and an accumulated logweight.

Without reference the particle draws from the prior. Given a retained trajectory's raw values it becomes a conditional-SMC reference pinned to that trajectory, erroring if its execution reaches an address that trajectory lacks or finishes without reaching one it has.

Prior()

Algorithm for sampling from the prior.

Every draw is an independent draw from the prior, so there is no starting point for initial_params to set; passing one warns and has no effect. To hold a variable at a value, condition or fix the model instead.

ProduceLogLikelihoodAccumulator{T} <: LogProbAccumulator{T}

A likelihood accumulator that Libtask.produces each increment as it accumulates it, which is what turns a model evaluation into a particle filter: one produce per likelihood term, so the sweep sees one filtering step per observe. Substituting it for LogLikelihoodAccumulator is the only thing that distinguishes a particle's varinfo.

RepeatSampler <: AbstractMCMC.AbstractSampler

A RepeatSampler is a container for a sampler and a number of times to repeat it.

Fields

• sampler: The sampler to repeat

• num_repeat: The number of times to repeat the sampler

Examples

repeated_sampler = RepeatSampler(sampler, 10)
# The initial step is a single step of `sampler`; it is the steps from a state that repeat.
_, state = AbstractMCMC.step(rng, model, repeated_sampler)
AbstractMCMC.step(rng, model, repeated_sampler, state) # take 10 steps of `sampler`

ReshapedBlock(variable::VarName, change::_BlockChange)

Say why a Gibbs component's parameter layout differs from the one it last saw:

• _BLOCK_JOINED: variable is now in the block;
• _BLOCK_LEFT: variable is no longer in the block;
• _BLOCK_REKEYED: variable is now written by a different tilde statement;
• _BLOCK_RESPECIFIED: variable now links to a different width.

Every case invalidates state sized for the block. change otherwise only shapes the error message.

Defined here rather than beside the rest of the Gibbs interface because the samplers that dispatch on it are loaded first. See the five-argument gibbs_update_state!! for what it is for.

struct SMC{R<:Turing.Inference.AbstractResampler} <: Turing.Inference.ParticleInference

Sequential Monte Carlo sampler.

Fields

• resampler::Turing.Inference.AbstractResampler: resampling scheme

• multithreaded::Bool: reweight the particles across threads within each sweep

SMC([resampler = ESSThresholdResampler(0.5)]; multithreaded = false)
SMC([scheme = StratifiedResampler(), ]threshold; multithreaded = false)

Sequential Monte Carlo sampler. By default stratified resampling is triggered whenever the effective sample size drops below half the number of particles.

Set multithreaded = true to evaluate the particles across threads within each sweep; results are unchanged (start Julia with multiple threads, e.g. julia -t auto, for this to have effect). Threads are the only option here: a suspended particle is a live Libtask task and cannot be serialised, so a single sweep cannot be spread across processes. Passing MCMCThreads() or MCMCDistributed() to sample parallelises whole chains instead, which is a separate axis and composes with this one.

The resampling scheme types (StratifiedResampler, SystematicResampler, MultinomialResampler, ESSThresholdResampler) are not exported; refer to them as e.g. Turing.Inference.SystematicResampler.

SMCContext

Leaf context marking a model evaluation as a particle-filter step: tilde_assume!! draws from the prior using the particle's own generator – or, for a conditional-SMC reference, reuses the retained trajectory's value at that address – and tilde_observe!! scores the observation and Libtask.produces the increment as the particle's weight.

Snapshot(values, layouts)

The values Gibbs threads through a sweep, with the linked layout of each tilde statement that produced them.

Both halves come from one model evaluation and reached_values is the only thing that builds one, so block_fingerprint can compare a block's layout either side of a step without evaluating the model again.

Stratified resampling: one independent uniform per stratum of width 1/n.

Systematic resampling: one shared uniform placed on a regular grid of n points.

_cond_dist_for(cond_dists, vn)

The conditional distribution governing the tilde statement at vn.

Keyed by subsumption, not by equality. The conditional may be stored at a coarser VarName than the statement uses – theta[:] ~ MvNormal(...) is one distribution for the whole of theta, and the single-distribution form of get_cond_dists never sees the tilde key – so an exact lookup raised a MethodError from inside the distribution.

_convert_initial_params(initial_params)

Convert initial_params to a DynamicPPl.AbstractInitStrategy if it is not already one, or throw a useful error message.

Turing.Inference._default_parameter_values(state)

The answer for a state with no gibbs_get_parameter_values method of its own.

An AbstractVarInfo carries its ~ values in a RawValueAccumulator, so there is one; any other state has to say for itself.

_require_owned(spl, sampler, varnames, leaf, what)

Throw unless leaf, whose existence just changed, is sampled by the component that changed it.

A variable coming or going under one component's proposal while another samples it means that other component conditions on a variable absent from the state being proposed, and the chain comes back biased. what says which direction it moved.

Turing.Inference.gibbs_update_state!!(
    sampler::AbstractSampler, state, model::Model, global_vals::VarNamedTuple,
    reshaped::ReshapedBlock
)

Update the state of a Gibbs component sampler whose block now holds a different set of variables from the one it last stepped.

That happens when another component samples the variable deciding whether one of this block's variables exists:

@model function f(y)
    b ~ Bernoulli(0.3)
    θ = zeros(2)
    θ[1] ~ Normal()
    b == 1 && (θ[2] ~ Normal())
    return y ~ Normal(sum(θ), 0.5)
end

Gibbs((@varname(b), @varname(θ)) => PG(20), @varname(θ) => HMC(0.1, 5))

PG samples b and θ together, then HMC samples θ again. When PG's step flips b, θ[2] appears or vanishes, so by the time HMC next steps its block has a different shape from the one it last saw. Each step still sees one fixed shape, so the scheme is valid; HMC only has to rebuild the parameter layout and phasepoint it cached for the old shape. Without a flip, θ never changes shape and the case never arises.

Defaults to throwing, because most components carry something so sized: an adapted mass matrix, a flat parameter layout, a set of prior means. Implementing this method is how a sampler says that it copes – the implementation is the declaration, so there is no separate trait that could fall out of step with what the sampler actually does, and a sampler written before this method existed keeps the safe answer.

Turing's own answers: MH, PG and GibbsConditional delegate to the four-argument form, caching nothing shaped like the block. A Hamiltonian that is not adapting rebuilds both the parameter layout and the phasepoint, which covers HMC and also NUTS(0, δ) and HMCDA(0, δ, λ), whose states carry NoAdaptation. An adapting NUTS or HMCDA does not implement it, since gen_metric renews an AdaptiveHamiltonian's metric from state.adaptor, whose mass matrix is sized for the shape it adapted to; resetting the adaptation instead would discard it mid-chain, which is a decision for the user rather than a detail of conditioning. Nor does ESS, whose prior means are gathered for the block it was built on, nor externalsampler, whose wrapped state is opaque.

_same_layout(before, now, linked)

Whether two (dist, val) records occupy the same layout, measured in linked space when the component holds one (see keeps_linked_layout) and at the values' own shape otherwise.

An unchanged distribution settles it either way without measuring: the layout is a function of the distribution, and measuring in linked space means deriving the linking transform, which under fix_transforms the caller has asked not to pay for once per tilde per step.

When they differ the width is measured by linking the value, not inferred from the family or the bijector's type. Those are proxies that diverge from the layout: keying on the family refused an adapting NUTS for a Normal()/TDist(3) branch whose block had not moved.

_shape_template(dist)

A container shaped and typed like dist, for templated_setindex!! to size an array by, or NoTemplate() when the distribution says nothing about a shape. Only its shape and element type are read, so it is left uninitialised.

The shape has to come from the distribution being written, not from the values the state holds: those describe the previous step, and another component may have changed the block's dimension since, in which case templating a two-element result onto a one-element hint throws a BoundsError from inside the setindex.

NoTemplate() alone is not enough either. A VarName carrying a Colon – m[:] ~ MvNormal(...) – gives the setindex nothing to infer a size from and it refuses outright. A multivariate or matrix-variate distribution knows its own size, which is exactly the missing hint; a univariate one at an indexed VarName does not, and does not need it.

advance!(particle) -> Union{Real,Nothing}

Run the particle to its next observe, returning the incremental log-likelihood, or nothing once the model finishes.

Turing.Inference.allow_discrete_variables(spl::AbstractSampler)

Whether spl can sample a model with discrete variables.

Defaults to true. Gradient-based samplers override it to false, and sample checks the model against it before starting.

Turing.Inference.allow_varying_dimension(spl::AbstractSampler)

Whether spl can move between supports within one of its own steps, that is, sample a block whose set of variables its own proposal changes.

Defaults to false, because proposing between two supports takes a construction built for it. A LogDensityFunction's layout is fixed for the whole of a step, so it has no slot for a variable the proposal reaches part-way through: a leapfrog step that crosses into another branch raises KeyError from DynamicPPL rather than reaching this check at all. PG and CSMC rebuild their trace each sweep, drawing whatever the model reaches, so they can.

MH answers true, but not because every crossing is safe for it. Whether one is depends on the variable that moved: drawn from its prior, the proposal density cancels against the prior and the ratio collapses to a likelihood ratio, which is defined across dimensions; proposed from, it does not cancel. That is a fact about one variable in one evaluation, which a trait asked once cannot supply, so MH answers true here and refuses the invalid crossings itself. A component whose answer really is uniform should give it here.

This is a different question from being handed a block at a new shape between one's own steps, which another component's step can do and which a component answers by implementing gibbs_update_state!! for a ReshapedBlock.

Returning true is a claim about the algorithm, and it carries an obligation: the component must have coherent semantics for a variable it samples going away and later coming back. PG and CSMC do. The reference particle replays the retained trajectory's values, so it reproduces that execution exactly and cannot reach an address the trajectory lacks, which it would otherwise refuse; the remaining particles draw from the prior and are free to reach a different set of addresses, which is what lets the block move between supports at all.

Returning true is not by itself enough. For an array whose length varies (for i in 1:n; x[i] ~ ...; end under a random n), the block has to hold n as well: with n in another component, the conditioned x[i] become observations whose number depends on it, and PG refuses with "the number of observations must not be random".

block_fingerprint(snapshot, varnames)

The shape of a component's block: the leaves it holds, and the linked width of each tilde statement that produced them.

Leaves alone are not enough. A tilde statement that changes which distribution it draws from can keep its leaf count and still move the block's linked dimension – x ~ Dirichlet(3) occupies two numbers and x ~ MvNormal(zeros(3), I) three, both with three leaves – and a component holding a parameter layout has to be told.

Only the width, so a block whose distributions change without moving it is not a reshape: x ~ truncated(Normal(); lower=a) with a in another block is relinked every sweep and still occupies one number, and refusing that would rule out a partition that samples perfectly well.

build_values_vnt(model::DynamicPPL.Model)

Build a VarNamedTuple of the values of every variable this component conditions on: those supplied as model arguments, and those Gibbs, the user, or fix conditioned.

merge is right-biased and replaces a whole key, so an array argument with a missing element – whose other elements are observations and whose missing one Gibbs conditions on the current draw – loses its observations to the partially-set conditioned value. Those elements are put back afterwards, rather than merging leaf by leaf throughout, so that a value stored under one key stays under one key: get_cond_dists sees these keys.

A missing in the context is put back the same way. condition(model; y=missing) leaves y absent as far as model execution is concerned, so the model uses the argument it was given, and get_cond_dists must see that argument rather than the missing the merge would otherwise leave in its place.

check_all_variables_handled(vnt, spl::Gibbs)

Check that every variable in vnt belongs to a component.

A key of vnt no declared varname subsumes is examined leaf by leaf, because ownership of a value stored as one key is a property of its leaves. One component owning every leaf can be handed the whole value, so that partition is fine; several owning parts of it is one Gibbs cannot express, since freeing one part of a single stored value is exactly what conditioned_values cannot do; and a leaf no component owns is one the user left out.

check_block_nonempty(varnames, block)

Throw if the model currently reaches none of the variables a component samples.

A component with nothing to sample is refused rather than skipped. MH would in fact sample such a partition correctly, and skipping the component for that sweep is the sensible semantics, but a LogDensityFunction over no variables is ill-formed and fails deep inside with VectorEvaluator requires a vector of floating-point values. Refusing uniformly is loud and easy to revisit; letting it through for some samplers and not others is neither.

check_no_missing_arguments(model)

Refuse a model with an argument bound to missing, or holding one.

Gibbs conditions every component on the variables it does not sample, and conditioning cannot reach a variable that is a model argument: the compiler reads the argument directly for such a tilde and never consults the condition context, so the missing arrives at the likelihood and throws MethodError: no method matching loglikelihood(::Normal{Float64}, ::Missing). An element of an array argument goes the same way as a whole one.

This is a capability GibbsContext had and conditioning does not: it made the variable an assumption and supplied the value itself. Restoring it needs condition to take precedence over a model argument (DynamicPPL.jl#1462, unmerged), and missing as a latent marker is due for deprecation anyway (DynamicPPL.jl#1464), so Gibbs refuses instead of failing inside the likelihood. Other samplers are unaffected.

Deliberately coarse: an argument that is missing but never the left-hand side of a tilde is refused too, though it would sample. Telling the two apart needs an evaluation, and a false refusal that names the argument is easier to act on than a MethodError from inside DynamicPPL.

check_reported_variables(varnames, conditioned, report)

Throw if report contains a value for a variable the component did not sample and was conditioned on.

gibbs_get_parameter_values' contract is that a component reports only what it samples, and merge gives the report priority, so a foreign value silently replaces another block's – either frozen for the rest of the run, or re-sampled afterwards from a corrupted conditional, which looks perfectly plausible in the chain.

A component may legitimately report a variable it does not sample when that variable did not exist to be conditioned on and its own step brought it into being: it is assumed rather than observed, so it lands in the component's accumulator. That case is not an error here – check_variable_set's appearance branch diagnoses it, and names the component that should have sampled it. Only a variable that was conditioned can have been overwritten.

check_variable_set(spl, sampler, varnames, old, new, proposed)

Check that a change in the set of variables is one the component that made it may make.

This is a best-effort check, not an enforcement of the rule in the Gibbs docstring. It compares the snapshots either side of a component's step, so it sees only what that component accepted: a proposal that crossed and was rejected leaves both snapshots equal and passes. In practice it is reliable, because the deciding variable usually does move – on the bad_dimension model in the Gibbs docstring it refused the split partition on all 40 seeds tried, at 50, 300 and 2000 draws alike – but that is a fact about those chains, not a guarantee. Returning normally therefore says nothing about whether the partition is valid; meeting the requirement is the caller's job.

Only a variable appearing or leaving is checked. A component moving the support of another block's variable is permitted, and deliberately unchecked; see the warning in the Gibbs docstring for what that costs and who owns it.

Two conditions have to hold together, and anything else throws:

• the component taking the step declares allow_varying_dimension; and
• the variable that came or went is one that same component samples.

The first is needed because during a component's step every variable it does not sample is conditioned, and so cannot move: if the set of variables the model reaches changed, it was that component's own proposal that moved between two different supports. The second is needed because a variable coming and going under one component's proposal while another component samples it means that other component conditions on a variable absent from the state being proposed – the two blocks disagree about the support, and the chain comes back biased even though both samplers can handle varying dimension on their own.

Together they say: the variable and whatever decides whether it exists belong in one block, sampled by a component built for it. The cost of getting this wrong is not a crash but a wrong answer: on the dynamic_bernoulli_normal model in the tests, where b decides whether θ[2] exists, Gibbs(@varname(b) => MH(), @varname(θ) => PG(20)) used to sample and return P(b=1) between 0.0 and 0.18 against an exact 0.394, while Gibbs((@varname(b), @varname(θ)) => PG(20)) gives 0.38.

Both directions are checked, and symmetrically: old and new are both the set of variables a model evaluation reached (see reached_values), so a leaf in one and not the other came or went during this step, and nothing else can have caused it – every variable outside the component's block was conditioned.

Taking both sets from the model, rather than from what the component reported, is what makes the verdict independent of the draw Gibbs initialises with and of the component's own bookkeeping. A component is free to keep a value for a variable it is not currently sampling; were that value taken into the snapshot it would be conditioned into the step of whichever component decides the variable's existence, that step would never assume it again, and the appearance branch would be dead for it – so a split partition would be rejected only for those seeds whose initialising draw happened to miss the variable.

check_walkers_same_layout(linked_vis)

Throw unless every walker occupies the same parameter layout.

The stretch move interpolates between two walkers' position vectors, and the single LogDensityFunction is built from the first walker, so every walker's vector is decoded against that one layout. All of them therefore have to hold the same variables in the same order at the same widths. On a model whose set of variables, or their sizes, depends on its own draws, walkers initialised from the prior do not.

Each variable's name, order, and width are compared, rather than only the names or the total width. Comparing names alone passed walkers agreeing on names while a variable's dimension differed, and names with the total width passed two variables trading dimensions – x of 2 and y of 3 against x of 3 and y of 2 – both of which then failed inside the decode or proposal.

Widths and not transforms: the layout fixes each variable's range, while its link is re-derived at every evaluation, so walkers holding one variable under different transforms still sample correctly, and comparing transforms would refuse them.

This is necessary rather than sufficient, and only the initial walkers are examined. A proposal can still cross into a branch none of them started in, and the decode then fails on a variable the layout has no range for – m ~ Normal(); m > 0 ? (x ~ Normal()) : (y ~ Normal()), started entirely in m > 0, raises KeyError: key y not found once a proposal reaches m < 0. Catching that would mean validating every evaluation. A model whose layout varies at all is best not sampled with Emcee.

component_stats(spl::Gibbs, states)

Collect the component samplers' statistics into one NamedTuple, prefixing each with the symbols of the variables that component samples, so that two components reporting e.g. acceptance_rate do not collide. Components sampling the same symbols are further distinguished by their index.

The prefix uses each variable's symbol rather than its whole VarName because chain packages read variable structure back out of these names: MCMCChains.namesingroup(chn, :x) matches anything beginning x[, so a statistic named x[1]_acceptance_rate would be served up as one of x's draws. A symbol carries no optic, so it cannot be parsed as part of another variable.

conditioned_values(global_vnt, target_variables)

Return the values in global_vnt for every variable not sampled by this Gibbs component, i.e. the ones it conditions on.

A component may own part of a value the model stores as a unit, which cannot be expressed: under x ~ MvNormal(zeros(2), I) the values hold x as one key, so freeing it for a component that samples x[1] would free x[2] too. This throws rather than hand a component a larger block than it owns. Whether a partition is expressible is a property of the model's tilde statements and not of the samplers – element-wise x[1] ~ and x[2] ~ give a key each and split cleanly, because the snapshot comes from a model evaluation rather than from what a component reports. Sampling a genuinely unsplittable block means writing a sampler against the AbstractMCMC interface directly.

Conditioned variables reach tilde_observe!!, so particle samplers reweight on them. That is what makes the component's target distribution correct: a conditioned variable the target depends on must reweight the sweep, and one it does not contributes the same increment to every particle, which ESS-gated resampling ignores.

externalsampler(
    sampler::AbstractSampler;
    adtype=AutoForwardDiff(),
    unconstrained=AbstractMCMC.requires_unconstrained_space(sampler),
)

Wrap a sampler so it can be used as an inference algorithm.

Arguments

• sampler::AbstractSampler: The sampler to wrap.

Keyword Arguments

• adtype::ADTypes.AbstractADType=ADTypes.AutoForwardDiff(): The automatic differentiation (AD) backend to use.
• unconstrained::Bool=AbstractMCMC.requires_unconstrained_space(sampler): Whether the sampler requires unconstrained space.

find_initial_params_ldf(rng, ldf, init_strategy; max_attempts=1000)

Given a LogDensityFunction and an initialization strategy, attempt to find valid initial parameters by sampling from the initialization strategy and checking that the log density (and gradient, if available) are finite. If valid parameters are not found after max_attempts, throw an error.

fork(particle, rng)

Copy particle into an independent, reseeded continuation. deepcopy forks the underlying TapedTask (Libtask defines copy as deepcopy) and preserves the task↔particle back-reference; reseed! then gives it its own random stream.

Turing.Inference.gibbs_get_parameter_values(state)

Return a VarNamedTuple containing the parameter values of all variables in the sampler state.

Turing's Gibbs sampler maintains, at all points during the sampling process, a single global VarNamedTuple that contains the raw values for all variables in the model. During the sampling process, it calls each component sampler in turn and rebuilds that VarNamedTuple around the new raw values returned by each sampler.

This function is used to pass that information from a component sampler to the Gibbs sampler. Note that this means that the VarNamedTuple returned by this function should only contain raw values for the variables that the component sampler is responsible for sampling, and should not contain any values for other variables. In particular it must leave out := quantities: they are not variables, so Gibbs would take one appearing inside a branch for a variable that appeared mid-run. DynamicPPL.get_parameter_values returns exactly the ~ values of a state whose accumulator holds both.

A step need not reach every variable the component samples: a target inside a branch the model did not take has no value that sweep, and leaving it out is the right answer. A component that wants to remember such a value – to reuse it if the branch comes back – should keep it in its own state and not report it. Gibbs does not read the report to decide which variables exist: it evaluates the model (see reached_values), because existence is a property of the model at the current values and not of a component's bookkeeping. Reporting a value for a variable the model no longer reaches has no effect on the snapshot, and reporting one for a variable the component does not sample is an error.

Turing.Inference.gibbs_get_raw_values(state)

Deprecated name for gibbs_get_parameter_values, still honoured so that a sampler written against it keeps working. Calls through to the new name, so a state that defines only the new one can still be asked by the old.

Turing.Inference.gibbs_get_stats(state)

Return a NamedTuple of sampler statistics (acceptance rates, step sizes, and so on) for the last step taken from state.

Gibbs discards its component samplers' transitions – reading parameters off them would cost a model re-evaluation – so a component that wants its statistics to reach the chain has to carry them on its state. Defaults to no statistics.

gibbs_initial_values(rng, model, spl, initial_params)

Return the values the sweep starts from, and check that every one of them has a component.

Take the first step of MCMC for the first component sampler, and call the same function recursively on the remaining samplers, until no samplers remain. Return the global VNT and a tuple of initial states for all component samplers.

The step_function argument should always be either AbstractMCMC.step or AbstractMCMC.step_warmup.

gibbs_recompute_ldf_and_params(
    old_ldf::LogDensityFunction,
    model::Model,
    global_vals::VarNamedTuple,
    extra_accs=()
)

Shared helper that is used in gibbs_update_state!! for any sampler that uses a LogDensityFunction.

Creates a new LogDensityFunction from the newly conditioned model, then reevaluates the model to obtain the correct vectorised parameters corresponding to the raw values in global_vals.

If extra information is needed (e.g. log-probabilities), extra_accs can be used to pass in other accumulators to be used in the same model evaluation, to avoid having to recompute them later.

Returns (new_ldf, new_params, accs) where accs is the set of accumulators after evaluation, from which extra accumulators (e.g. LogLikelihoodAccumulator) can be read.

The flat layout is built from global_vals rather than reused from old_ldf. A layout is keyed by tilde statement and carries each variable's transform, so a reused one is stale the moment the model takes a different branch – not only when a variable joins or leaves the block, but also when the same variables arrive under different keys, or with a different transform and so a different linked dimension. Evaluating against a stale layout fails inside DynamicPPL naming no variable. Built from the values, the layout is identical whenever the execution path is, which is exactly when the component's own state depends on the ordering, and differs only when that state was going to be wrong anyway.

Run a Gibbs step for the first varname/sampler/state tuple, and recursively call the same function on the tail, until there are no more samplers left.

The step_function argument should always be either AbstractMCMC.step or AbstractMCMC.step_warmup.

Turing.Inference.gibbs_update_state!!(
    sampler::AbstractSampler, state, model::Model, global_vals::VarNamedTuple
)

Update the state of a Gibbs component sampler to be consistent with the new values in global_vals. Each sampler should implement a method for its respective state type.

Note that the model argument passed in here will be 'conditioned' on the new values inside global_vals. Thus, evaluating it will reflect the log-probability associated with the new values.

Exactly what this function should do will depend on what the sampler state contains, but for example, it will often mean:

• Updating any raw or vectorised values stored in the sampler state to be consistent with global_vals.
• Reevaluating the (new) model to update any cached log-probabilities.
• Updating any log-density callables (such as a DynamicPPL.LogDensityFunction) stored in the sampler state, to be consistent with the new model.

For examples of this, please see the implementations of this function for the samplers in Turing.jl. In particular, the HMC and ExternalSampler implementations work with LogDensityFunction and demonstrate how information such as that can be updated based on the new model.

See also the five-argument form, gibbs_update_state!!(sampler, state, model, global_vals, reshaped::ReshapedBlock), which Gibbs calls in place of this one when the block has a different set of variables from the one the sampler last stepped.

Turing.Inference.init_strategy(spl::AbstractSampler)

Get the default initialization strategy for a given sampler spl, i.e. how initial parameters for sampling are chosen if not specified by the user. By default, this is InitFromPrior(), which samples initial parameters from the prior distribution.

isgibbscomponent(spl::AbstractSampler)

Deprecated name for supports_gibbs, still honoured so that a sampler written against it keeps working.

Defined in terms of supports_gibbs rather than as a constant true, so that a wrapper delegating through this name – isgibbscomponent(w::MyWrapper) = isgibbscomponent(w.inner) – still gets the right answer for a sampler Turing declares unusable.

Whether particle is a conditional-SMC reference, i.e. pinned to a retained trajectory. Carried on the particle rather than inferred from its slot, so forking and resampling cannot get it wrong.

keeps_linked_layout(spl::AbstractSampler)

Whether spl carries a parameter vector in linked space, so that a change in a variable's linked width reshapes its block.

true by default, since a sampler that does carry one is broken by a reshape it is not told about. A sampler that keeps nothing between steps answers false and has its block compared at the values' native shape, which costs no Bijectors transform: deriving one is not free, and for a distribution that defines no link it is not possible.

An approximation of the question that decides anything, which is whether the sampler's ReshapedBlock method throws. That cannot be asked – the throwing fallback is a method like any other – so a sampler whose method merely delegates still answers true and still pays for a measurement whose verdict cannot change what it does. HMC is the case in point.

Turing.loadstate(chain::FlexiChain{<:VarName})

Extracts the last sampler state from a FlexiChain. This is the same function as FlexiChains.last_sampler_state.

model_argument_values(model)

The values model's arguments supply, keyword as well as positional.

model.args holds only the positional ones; a keyword argument, whether defaulted or passed, lands in model.defaults. Reading one and not the other is how a keyword argument became invisible to two separate callers, so both ask here.

A fresh counter-based generator for one particle, seeded from rng.

post_sample_hook(chain, sampler::AbstractSampler; kwargs...)

A post-sampling hook that can e.g. print info about the results of sampling.

Implementations of this should be careful to take kwargs... as keyword arguments instead of restricting the signature to specific keyword arguments. Right now, the only keyword argument that is passed to this function is verbose, but in the future additional keyword arguments may be passed here.

post_sample_hook(chain::FlexiChains.VNChain, sampler::Union{HMC,NUTS,HMCDA}; kwargs...)

Emit a warning message if there are divergent transitions in the chain.

reached_values(rng, model, proposed)

Return the values of the variables model reaches when evaluated at proposed.

This is how the snapshot after a component's step is built, rather than from the values the component reports. Which variables exist is a property of the model at the current values, so only an evaluation settles it, and a component's report cannot be trusted for it: a component is free to keep a value for a variable it is not currently sampling, which is useful to it and none of Gibbs's business, but were that value to enter the snapshot it would be conditioned into the step of whichever component decides the variable's existence, and check_variable_set would never see the variable come back.

The evaluation costs one model evaluation per component step, and draws from the prior only for a variable proposed has no value for. That happens when a component under-reports a variable it samples, or when a variable appeared during this step, and check_variable_set throws in both cases – so nothing that returns normally has consumed rng.

Draw n ancestor indices from 1:length(weights) with probabilities weights.

reseed!(particle, rng)

Restart particle as a fresh continuation seeded from rng, so that a particle descended from the reference stops reusing retained values and samples afresh. Mutates and returns particle.

Whether to resample given the normalized weights. Bare schemes always resample.

supports_gibbs(spl::AbstractSampler)

Return a boolean indicating whether spl is a valid component for a Gibbs sampler.

Defaults to true if no method has been defined for a particular sampler.

Effective sample size of a normalised weight vector, 1 / Σ wᵢ². Named for the weights to keep it distinct from MCMCDiagnosticTools.ess, which Turing re-exports and which measures a chain's autocorrelation rather than a population's weight degeneracy.