API: Turing.Inference

CSMC(...)

Equivalent to 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 │

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{S<:AbstractSampler,AD<:ADTypes.AbstractADType,Unconstrained}

Represents a sampler that is not an implementation of InferenceAlgorithm.

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

Gibbs(algs...)

Compositional MCMC interface. Gibbs sampling combines one or more sampling algorithms, each of which samples from a different set of variables in a model.

Example:

@model function gibbs_example(x)
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
end

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(HMC(0.2, 3, :v1), PG(20, :v2))

One can also pass the number of iterations for each Gibbs component using the following syntax:

• alg = Gibbs((HMC(0.2, 3, :v1), n_hmc), (PG(20, :v2), n_pg))

where n_hmc and n_pg are the number of HMC and PG iterations for each Gibbs iteration.

Tips:

• HMC and NUTS are fast samplers and can throw off particle-based

methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including more particles in the particle sampler.

GibbsConditional(sym, conditional)

A "pseudo-sampler" to manually provide analytical Gibbs conditionals to Gibbs. GibbsConditional(:x, cond) will sample the variable x according to the conditional cond, which must therefore be a function from a NamedTuple of the conditioned variables to a Distribution.

The NamedTuple that is passed in contains all random variables from the model in an unspecified order, taken from the VarInfo object over which the model is run. Scalars and vectors are stored in their respective shapes. The tuple also contains the value of the conditioned variable itself, which can be useful, but using it creates something that is not a Gibbs sampler anymore (see here).

Examples

α_0 = 2.0
θ_0 = inv(3.0)
x = [1.5, 2.0]
N = length(x)

@model function inverse_gdemo(x)
    λ ~ Gamma(α_0, θ_0)
    σ = sqrt(1 / λ)
    m ~ Normal(0, σ)
    @. x ~ $(Normal(m, σ))
end

# The conditionals can be formulated in terms of the following statistics:
x_bar = mean(x) # sample mean
s2 = var(x; mean=x_bar, corrected=false) # sample variance
m_n = N * x_bar / (N + 1)

function cond_m(c)
    λ_n = c.λ * (N + 1)
    σ_n = sqrt(1 / λ_n)
    return Normal(m_n, σ_n)
end

function cond_λ(c)
    α_n = α_0 + (N - 1) / 2 + 1
    β_n = s2 * N / 2 + c.m^2 / 2 + inv(θ_0)
    return Gamma(α_n, inv(β_n))
end

m = inverse_gdemo(x)

sample(m, Gibbs(GibbsConditional(:λ, cond_λ), GibbsConditional(:m, cond_m)), 10)

GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}}

Stores a VarInfo for use in sampling, and a Tuple of Samplers that the Gibbs sampler iterates through for each step!.

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.

IS()

Importance sampling algorithm.

Usage:

IS()

Example:

# Define a simple Normal model with unknown mean and variance.
@model function gdemo(x)
    s² ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt.(s))
    x[1] ~ Normal(m, sqrt.(s))
    x[2] ~ Normal(m, sqrt.(s))
    return s², m
end

sample(gdemo([1.5, 2]), IS(), 1000)

MH(space...)

Construct a Metropolis-Hastings algorithm.

The arguments space can be

• Blank (i.e. MH()), in which case MH defaults to using the prior for each parameter as the proposal distribution.
• A set of one or more symbols to sample with MH in conjunction with Gibbs, i.e. Gibbs(MH(:m), PG(10, :s))
• An iterable of pairs or tuples mapping a Symbol to a AdvancedMH.Proposal, Distribution, or Function that generates returns a conditional proposal distribution.
• A covariance matrix to use as for mean-zero multivariate normal proposals.

Examples

The default MH will draw proposal samples from the prior distribution using AdvancedMH.StaticProposal.

@model function gdemo(x, y)
    s² ~ InverseGamma(2,3)
    m ~ Normal(0, sqrt(s²))
    x ~ Normal(m, sqrt(s²))
    y ~ Normal(m, sqrt(s²))
end

chain = sample(gdemo(1.5, 2.0), MH(), 1_000)
mean(chain)

Alternatively, you can specify particular parameters to sample if you want to combine sampling from multiple samplers:

# Samples s² with MH and m with PG
chain = sample(gdemo(1.5, 2.0), Gibbs(MH(:s²), PG(10, :m)), 1_000)
mean(chain)

Specifying a single distribution implies the use of static MH:

# Use a static proposal for s² (which happens to be the same
# as the prior) and a static proposal for m (note that this 
# isn't a random walk proposal).
chain = sample(
    gdemo(1.5, 2.0),
    MH(
        :s² => InverseGamma(2, 3),
        :m => Normal(0, 1)
    ),
    1_000
)
mean(chain)

Specifying explicit proposals using the AdvancedMH interface:

# Use a static proposal for s² and random walk with proposal
# standard deviation of 0.25 for m.
chain = sample(
    gdemo(1.5, 2.0),
    MH(
        :s² => AdvancedMH.StaticProposal(InverseGamma(2,3)),
        :m => AdvancedMH.RandomWalkProposal(Normal(0, 0.25))
    ),
    1_000
)
mean(chain)

Using a custom function to specify a conditional distribution:

# Use a static proposal for s and and a conditional proposal for m,
# where the proposal is centered around the current sample.
chain = sample(
    gdemo(1.5, 2.0),
    MH(
        :s² => InverseGamma(2, 3),
        :m => x -> Normal(x, 1)
    ),
    1_000
)
mean(chain)

Providing a covariance matrix will cause MH to perform random-walk sampling in the transformed space with proposals drawn from a multivariate normal distribution. The provided matrix must be positive semi-definite and square:

# Providing a custom variance-covariance matrix
chain = sample(
    gdemo(1.5, 2.0),
    MH(
        [0.25 0.05;
         0.05 0.50]
    ),
    1_000
)
mean(chain)

MHLogDensityFunction

A log density function for the MH sampler.

This variant uses the set_namedtuple! function to update the VarInfo.

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.

PG(n, space...)
PG(n, [resampler = AdvancedPS.ResampleWithESSThreshold(), space = ()])
PG(n, [resampler = AdvancedPS.resample_systematic, ]threshold[, space = ()])

Create a Particle Gibbs sampler of type PG with n particles for the variables in space.

If the algorithm for the resampling step is not specified explicitly, systematic resampling is performed if the estimated effective sample size per particle drops below 0.5.

struct PG{space, R} <: Turing.Inference.ParticleInference

Particle Gibbs sampler.

Fields

• nparticles::Int64: Number of particles.

• resampler::Any: Resampling algorithm.

PolynomialStepsize(a[, b=0, γ=0.55])

Create a polynomially decaying stepsize function.

At iteration t, the step size is

\[a (b + t)^{-γ}.\]

Prior()

Algorithm for sampling from the prior.

SGHMC{AD,space}

Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) sampler.e

Fields

• learning_rate::Real

• momentum_decay::Real

• adtype::Any

Reference

Tianqi Chen, Emily Fox, & Carlos Guestrin (2014). Stochastic Gradient Hamiltonian Monte Carlo. In: Proceedings of the 31st International Conference on Machine Learning (pp. 1683–1691).

SGHMC(
    space::Symbol...;
    learning_rate::Real,
    momentum_decay::Real,
    adtype::ADTypes.AbstractADType = AutoForwardDiff(),
)

Create a Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) sampler.

If the automatic differentiation (AD) backend adtype is not provided, ForwardDiff with automatically determined chunksize is used.

Reference

Tianqi Chen, Emily Fox, & Carlos Guestrin (2014). Stochastic Gradient Hamiltonian Monte Carlo. In: Proceedings of the 31st International Conference on Machine Learning (pp. 1683–1691).

SGLD

Stochastic gradient Langevin dynamics (SGLD) sampler.

Fields

• stepsize::Any: Step size function.

• adtype::Any

Reference

Max Welling & Yee Whye Teh (2011). Bayesian Learning via Stochastic Gradient Langevin Dynamics. In: Proceedings of the 28th International Conference on Machine Learning (pp. 681–688).

SGLD(
    space::Symbol...;
    stepsize = PolynomialStepsize(0.01),
    adtype::ADTypes.AbstractADType = AutoForwardDiff(),
)

Stochastic gradient Langevin dynamics (SGLD) sampler.

By default, a polynomially decaying stepsize is used.

If the automatic differentiation (AD) backend adtype is not provided, ForwardDiff with automatically determined chunksize is used.

Reference

Max Welling & Yee Whye Teh (2011). Bayesian Learning via Stochastic Gradient Langevin Dynamics. In: Proceedings of the 28th International Conference on Machine Learning (pp. 681–688).

See also: PolynomialStepsize

SMC(space...)
SMC([resampler = AdvancedPS.ResampleWithESSThreshold(), space = ()])
SMC([resampler = AdvancedPS.resample_systematic, ]threshold[, space = ()])

Create a sequential Monte Carlo sampler of type SMC for the variables in space.

If the algorithm for the resampling step is not specified explicitly, systematic resampling is performed if the estimated effective sample size per particle drops below 0.5.

struct SMC{space, R} <: Turing.Inference.ParticleInference

Sequential Monte Carlo sampler.

Fields

• resampler::Any

predict([rng::AbstractRNG,] model::Model, chain::MCMCChains.Chains; include_all=false)

Execute model conditioned on each sample in chain, and return the resulting Chains.

If include_all is false, the returned Chains will contain only those variables sampled/not present in chain.

Details

Internally calls Turing.Inference.transitions_from_chain to obtained the samples and then converts these into a Chains object using AbstractMCMC.bundle_samples.

Example

julia> using Turing; Turing.setprogress!(false);
[ Info: [Turing]: progress logging is disabled globally

julia> @model function linear_reg(x, y, σ = 0.1)
           β ~ Normal(0, 1)

           for i ∈ eachindex(y)
               y[i] ~ Normal(β * x[i], σ)
           end
       end;

julia> σ = 0.1; f(x) = 2 * x + 0.1 * randn();

julia> Δ = 0.1; xs_train = 0:Δ:10; ys_train = f.(xs_train);

julia> xs_test = [10 + Δ, 10 + 2 * Δ]; ys_test = f.(xs_test);

julia> m_train = linear_reg(xs_train, ys_train, σ);

julia> chain_lin_reg = sample(m_train, NUTS(100, 0.65), 200);
┌ Info: Found initial step size
└   ϵ = 0.003125

julia> m_test = linear_reg(xs_test, Vector{Union{Missing, Float64}}(undef, length(ys_test)), σ);

julia> predictions = predict(m_test, chain_lin_reg)
Object of type Chains, with data of type 100×2×1 Array{Float64,3}

Iterations        = 1:100
Thinning interval = 1
Chains            = 1
Samples per chain = 100
parameters        = y[1], y[2]

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters     mean     std  naive_se     mcse       ess   r_hat
  ──────────  ───────  ──────  ────────  ───────  ────────  ──────
        y[1]  20.1974  0.1007    0.0101  missing  101.0711  0.9922
        y[2]  20.3867  0.1062    0.0106  missing  101.4889  0.9903

Quantiles
  parameters     2.5%    25.0%    50.0%    75.0%    97.5%
  ──────────  ───────  ───────  ───────  ───────  ───────
        y[1]  20.0342  20.1188  20.2135  20.2588  20.4188
        y[2]  20.1870  20.3178  20.3839  20.4466  20.5895

julia> ys_pred = vec(mean(Array(group(predictions, :y)); dims = 1));

julia> sum(abs2, ys_test - ys_pred) ≤ 0.1
true

dist_val_tuple(spl::Sampler{<:MH}, vi::VarInfo)

Return two NamedTuples.

The first NamedTuple has symbols as keys and distributions as values. The second NamedTuple has model symbols as keys and their stored values as values.

externalsampler(sampler::AbstractSampler; adtype=AutoForwardDiff(), unconstrained=true)

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=true: Whether the sampler requires unconstrained space.

getparams(model, t)

Return a named tuple of parameters.

gibbs_rerun(prev_alg, alg)

Check if the model should be rerun to recompute the log density before sampling with the Gibbs component alg and after sampling from Gibbs component prev_alg.

By default, the function returns true.

gibbs_state(model, sampler, state, varinfo)

Return an updated state, taking into account the variables sampled by other Gibbs components.

Arguments

• model: model targeted by the Gibbs sampler.
• sampler: the sampler for this Gibbs component.
• state: the state of sampler computed in the previous iteration.
• varinfo: the variables, including the ones sampled by other Gibbs components.

gibbs_varinfo(model, sampler, state)

Return the variables corresponding to the current state of the Gibbs component sampler.

group_varnames_by_symbol(vns)

Group the varnames by their symbol.

Arguments

• vns: Iterable of VarName.

Returns

• OrderedDict{Symbol, Vector{VarName}}: A dictionary mapping symbol to a vector of varnames.

isgibbscomponent(alg)

Determine whether algorithm alg is allowed as a Gibbs component.

mh_accept(logp_current::Real, logp_proposal::Real, log_proposal_ratio::Real)

Decide if a proposal $x'$ with log probability $\log p(x') = logp_proposal$ and log proposal ratio $\log k(x', x) - \log k(x, x') = log_proposal_ratio$ in a Metropolis-Hastings algorithm with Markov kernel $k(x_t, x_{t+1})$ and current state $x$ with log probability $\log p(x) = logp_current$ is accepted by evaluating the Metropolis-Hastings acceptance criterion

\[\log U \leq \log p(x') - \log p(x) + \log k(x', x) - \log k(x, x')\]

for a uniform random number $U \in [0, 1)$.

recompute_logprob!!(rng, model, sampler, state)

Recompute the log-probability of the model based on the given state and return the resulting state.

requires_unconstrained_space(sampler::ExternalSampler)

Return true if the sampler requires unconstrained space, and false otherwise.

set_namedtuple!(vi::VarInfo, nt::NamedTuple)

Places the values of a NamedTuple into the relevant places of a VarInfo.

transitions_from_chain(
    [rng::AbstractRNG,]
    model::Model,
    chain::MCMCChains.Chains;
    sampler = DynamicPPL.SampleFromPrior()
)

Execute model conditioned on each sample in chain, and return resulting transitions.

The returned transitions are represented in a Vector{<:Turing.Inference.Transition}.

Details

In a bit more detail, the process is as follows:

1. For every sample in chain
   1. For every variable in sample
      1. Set variable in model to its value in sample
   2. Execute model with variables fixed as above, sampling variables NOT present in chain using SampleFromPrior
   3. Return sampled variables and log-joint

Example

julia> using Turing

julia> @model function demo()
           m ~ Normal(0, 1)
           x ~ Normal(m, 1)
       end;

julia> m = demo();

julia> chain = Chains(randn(2, 1, 1), ["m"]); # 2 samples of `m`

julia> transitions = Turing.Inference.transitions_from_chain(m, chain);

julia> [Turing.Inference.getlogp(t) for t in transitions] # extract the logjoints
2-element Array{Float64,1}:
 -3.6294991938628374
 -2.5697948166987845

julia> [first(t.θ.x) for t in transitions] # extract samples for `x`
2-element Array{Array{Float64,1},1}:
 [-2.0844148956440796]
 [-1.704630494695469]