API

@model(expr[, warn = false])

Macro to specify a probabilistic model.

If warn is true, a warning is displayed if internal variable names are used in the model definition.

Examples

Model definition:

@model function model(x, y = 42)
    ...
end

To generate a Model, call model(xvalue) or model(xvalue, yvalue).

struct Model{F,argnames,defaultnames,missings,Targs,Tdefaults,Ctx<:AbstractContext,Threaded}
    f::F
    args::NamedTuple{argnames,Targs}
    defaults::NamedTuple{defaultnames,Tdefaults}
    context::Ctx=DefaultContext()
end

A Model struct with model evaluation function of type F, arguments of names argnames types Targs, default arguments of names defaultnames with types Tdefaults, missing arguments missings, and evaluation context of type Ctx.

Here argnames, defaultargnames, and missings are tuples of symbols, e.g. (:a, :b). context is by default DefaultContext().

An argument with a type of Missing will be in missings by default. However, in non-traditional use-cases missings can be defined differently. All variables in missings are treated as random variables rather than observations.

The Threaded type parameter indicates whether the model requires threadsafe evaluation (i.e., whether the model contains statements which modify the internal VarInfo that are executed in parallel). By default, this is set to false.

The default arguments are used internally when constructing instances of the same model with different arguments.

Examples

julia> Model(f, (x = 1.0, y = 2.0))
Model{typeof(f),(:x, :y),(),(),Tuple{Float64,Float64},Tuple{}}(f, (x = 1.0, y = 2.0), NamedTuple())

julia> Model(f, (x = 1.0, y = 2.0), (x = 42,))
Model{typeof(f),(:x, :y),(:x,),(),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))

julia> Model{(:y,)}(f, (x = 1.0, y = 2.0), (x = 42,)) # with special definition of missings
Model{typeof(f),(:x, :y),(:x,),(:y,),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))

(model::Model)([rng, varinfo])

Sample from the prior of the model with random number generator rng.

Returns the model's return value.

Note that calling this with an existing varinfo object will mutate it.

nameof(model::Model)

Get the name of the model as Symbol.

getargnames(model::Model)

Get a tuple of the argument names of the model.

getmissings(model::Model)

Get a tuple of the names of the missing arguments of the model.

contextualize(model::Model, context::AbstractContext)

Return a new Model with the same evaluation function and other arguments, but with its underlying context set to context.

setthreadsafe(model::Model, threadsafe::Bool)

Returns a new Model with its threadsafe flag set to threadsafe.

Threadsafe evaluation ensures correctness when executing model statements that mutate the internal VarInfo object in parallel. For example, this is needed if tilde-statements are nested inside Threads.@threads or similar constructs.

It is not needed for generic multithreaded operations that don't involve VarInfo. For example, calculating a log-likelihood term in parallel and then calling @addlogprob! outside of the parallel region is safe without needing to set threadsafe=true.

It is also not needed for multithreaded sampling with AbstractMCMC's MCMCThreads().

Setting threadsafe to true increases the overhead in evaluating the model. Please see the Turing.jl docs for more details.

requires_threadsafe(model::Model)

Return whether model has been marked as needing threadsafe evaluation (using setthreadsafe).

rand([rng=Random.default_rng()], [T=NamedTuple], model::Model)

Generate a sample of type T from the prior distribution of the model.

logprior(model::Model, params)
logprior(model::Model, varinfo::AbstractVarInfo)

Return the log prior probability of variables params for the probabilistic model, or the log prior of the data in varinfo if provided.

Note that this probability always refers to the parameters in unlinked space, i.e., the return value of logprior does not depend on whether VarInfo has been linked or not.

See also logjoint and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logprior(demo([1.0]), (m = 100.0, ))
-5000.918938533205

julia> # Using a `OrderedDict`.
       logprior(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-5000.918938533205

julia> # Truth.
       logpdf(Normal(), 100.0)
-5000.918938533205

logprior(model::Model, values::Union{NamedTuple,AbstractDict})

Return the log prior probability of variables values for the probabilistic model.

See also logjoint and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logprior(demo([1.0]), (m = 100.0, ))
-5000.918938533205

julia> # Using a `OrderedDict`.
       logprior(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-5000.918938533205

julia> # Truth.
       logpdf(Normal(), 100.0)
-5000.918938533205

logprior(model::DynamicPPL.Model, chain::MCMCChains.Chains)

Return an array of log prior probabilities evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> logprior(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -3.2956988239086447
 -3.2956988239086447
 -3.2956988239086447

loglikelihood(model::Model, params)
loglikelihood(model::Model, varinfo::AbstractVarInfo)

Return the log likelihood of variables params for the probabilistic model, or the log likelihood of the data in varinfo if provided.

See also logjoint and logprior.

Examples

```jldoctest; setup=:(using Distributions) julia> @model function demo(x) m ~ Normal() for i in eachindex(x) x[i] ~ Normal(m, 1.0) end end demo (generic function with 2 methods)

julia> # Using a NamedTuple. loglikelihood(demo([1.0]), (m = 100.0, )) -4901.418938533205

julia> # Using a OrderedDict. loglikelihood(demo([1.0]), OrderedDict(@varname(m) => 100.0)) -4901.418938533205

julia> # Truth. logpdf(Normal(100.0, 1.0), 1.0) -4901.418938533205

loglikelihood(model::Model, values::Union{NamedTuple,AbstractDict})

Return the log likelihood of variables values for the probabilistic model.

See also logjoint and logprior.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       loglikelihood(demo([1.0]), (m = 100.0, ))
-4901.418938533205

julia> # Using a `OrderedDict`.
       loglikelihood(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-4901.418938533205

julia> # Truth.
       logpdf(Normal(100.0, 1.0), 1.0)
-4901.418938533205

loglikelihood(model::DynamicPPL.Model, chain::MCMCChains.Chains)

Return an array of log likelihoods evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> loglikelihood(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -2.1447298858494
 -2.1447298858494
 -2.1447298858494

logjoint(model::Model, params)
logjoint(model::Model, varinfo::AbstractVarInfo)

Return the log joint probability of variables params for the probabilistic model, or the log joint of the data in varinfo if provided.

Note that this probability always refers to the parameters in unlinked space, i.e., the return value of logjoint does not depend on whether VarInfo has been linked or not.

See also logprior and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logjoint(demo([1.0]), (m = 100.0, ))
-9902.33787706641

julia> # Using a `OrderedDict`.
       logjoint(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-9902.33787706641

julia> # Truth.
       logpdf(Normal(100.0, 1.0), 1.0) + logpdf(Normal(), 100.0)
-9902.33787706641

logjoint(model::Model, values::Union{NamedTuple,AbstractDict})

Return the log joint probability of variables values for the probabilistic model.

See logprior and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logjoint(demo([1.0]), (m = 100.0, ))
-9902.33787706641

julia> # Using a `OrderedDict`.
       logjoint(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-9902.33787706641

julia> # Truth.
       logpdf(Normal(100.0, 1.0), 1.0) + logpdf(Normal(), 100.0)
-9902.33787706641

logjoint(model::Model, chain::MCMCChains.Chains)

Return an array of log joint probabilities evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> logjoint(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -5.440428709758045
 -5.440428709758045
 -5.440428709758045

DynamicPPL.LogDensityFunction(
    model::Model,
    getlogdensity::Any=getlogjoint_internal,
    varinfo::AbstractVarInfo=VarInfo(model)
    accs::Union{NTuple{<:Any,AbstractAccumulator},AccumulatorTuple}=DynamicPPL.ldf_accs(getlogdensity);
    adtype::Union{ADTypes.AbstractADType,Nothing}=nothing,
)

A struct which contains a model, along with all the information necessary to:

• calculate its log density at a given point;
• and if adtype is provided, calculate the gradient of the log density at that point.

This information can be extracted using the LogDensityProblems.jl interface, specifically, using LogDensityProblems.logdensity and LogDensityProblems.logdensity_and_gradient. If adtype is nothing, then only logdensity is implemented. If adtype is a concrete AD backend type, then logdensity_and_gradient is also implemented.

getlogdensity should be a callable which takes a single argument: a VarInfo, and returns a Real corresponding to the log density of interest. There are several functions in DynamicPPL that are 'supported' out of the box:

• getlogjoint_internal: calculate the log joint, including the log-Jacobian term for any variables that have been linked in the provided VarInfo.
• getlogprior_internal: calculate the log prior, including the log-Jacobian term for any variables that have been linked in the provided VarInfo.
• getlogjoint: calculate the log joint in the model space, ignoring any effects of linking
• getlogprior: calculate the log prior in the model space, ignoring any effects of linking
• getloglikelihood: calculate the log likelihood (this is unaffected by linking, since transforms are only applied to random variables)

If you provide one of these functions, a VarInfo will be automatically created for you. If you provide a different function, you have to manually create a VarInfo and pass it as the third argument.

accs allows you to specify an AccumulatorTuple or a tuple of AbstractAccumulators which will be used when evaluating the log density. (Note that the accumulators from theVarInfoargument are discarded.) By default, this uses an internal function,DynamicPPL.ldf_accs`, which attempts to choose an appropriate set of accumulators based on which kind of log-density is being calculated.

If the adtype keyword argument is provided, then this struct will also store the adtype along with other information for efficient calculation of the gradient of the log density. Note that preparing a LogDensityFunction with an AD type AutoBackend() requires the AD backend itself to have been loaded (e.g. with import Backend).

Fields

Note that it is undefined behaviour to access any of a LogDensityFunction's fields, apart from:

• ldf.model: The original model from which this LogDensityFunction was constructed.
• ldf.adtype: The AD type used for gradient calculations, or nothing if no AD type was provided.

Extended help

Up until DynamicPPL v0.38, there have been two ways of evaluating a DynamicPPL model at a given set of parameters:

1. With unflatten!! + evaluate!! with DefaultContext: this stores a vector of parameters inside a VarInfo's metadata, then reads parameter values from the VarInfo during evaluation.

2. With InitFromParams: this reads parameter values from a NamedTuple or a Dict, and stores them inside a VarInfo's metadata.

In general, both of these approaches work fine, but the fact that they modify the VarInfo's metadata can often be quite wasteful. In particular, it is very common that the only outputs we care about from model evaluation are those which are stored in accumulators, such as log probability densities, or ValuesAsInModel.

To avoid this issue, we use OnlyAccsVarInfo, which is a VarInfo that only contains accumulators. It implements enough of the AbstractVarInfo interface to not error during model evaluation.

Because OnlyAccsVarInfo does not store any parameter values, when evaluating a model with it, it is mandatory that parameters are provided from outside the VarInfo, namely via InitContext.

The main problem that we face is that it is not possible to directly implement DynamicPPL.init(rng, vn, dist, strategy) for strategy::InitFromParams{<:AbstractVector}. In particular, it is not clear:

• which parts of the vector correspond to which random variables, and
• whether the variables are linked or unlinked.

Traditionally, this problem has been solved by unflatten!!, because that function would place values into the VarInfo's metadata alongside the information about ranges and linking. That way, when we evaluate with DefaultContext, we can read this information out again. However, we want to avoid using a metadata. Thus, here, we extract this information from the VarInfo a single time when constructing a LogDensityFunction object. Inside the LogDensityFunction, we store a mapping from VarNames to ranges in that vector, along with link status.

When evaluating the model, this allows us to combine the parameter vector together with those ranges to create an InitFromParams{VectorWithRanges}, which lets us very quickly read parameter values from the vector.

Note that this assumes that the ranges and link status are static throughout the lifetime of the LogDensityFunction object. Therefore, a LogDensityFunction object cannot handle models which have variable numbers of parameters, or models which may visit random variables in different orders depending on stochastic control flow. Indeed, silent errors may occur with such models. This is a general limitation of vectorised parameters: the original unflatten!! + evaluate!! approach also fails with such models.

OnlyAccsVarInfo

This is a wrapper around an AccumulatorTuple that implements the minimal AbstractVarInfo interface to work with the tilde_assume!! and tilde_observe!! functions for InitContext.

Note that this does not implement almost every other AbstractVarInfo interface function, and so using this with a different leaf context such as DefaultContext will result in errors.

Conceptually, one can also think of this as a VarInfo that doesn't contain a metadata field. This is also why it only works with InitContext: in this case, the parameters used for evaluation are supplied by the context instead of the metadata.

model | (x = 1.0, ...)

Return a Model which now treats variables on the right-hand side as observations.

See condition for more information and examples.

condition(model::Model; values...)
condition(model::Model, values::NamedTuple)

Return a Model which now treats the variables in values as observations.

See also: decondition, conditioned

Limitations

This does currently not work with variables that are provided to the model as arguments, e.g. @model function demo(x) ... end means that condition will not affect the variable x.

Therefore if one wants to make use of condition and decondition one should not be specifying any random variables as arguments.

This is done for the sake of backwards compatibility.

Examples

Simple univariate model

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> m, x = model(); (m ≠ 1.0 && x ≠ 100.0)
true

julia> # Create a new instance which treats `x` as observed
       # with value `100.0`, and similarly for `m=1.0`.
       conditioned_model = condition(model, x=100.0, m=1.0);

julia> m, x = conditioned_model(); (m == 1.0 && x == 100.0)
true

julia> # Let's only condition on `x = 100.0`.
       conditioned_model = condition(model, x = 100.0);

julia> m, x =conditioned_model(); (m ≠ 1.0 && x == 100.0)
true

julia> # We can also use the nicer `|` syntax.
       conditioned_model = model | (x = 100.0, );

julia> m, x = conditioned_model(); (m ≠ 1.0 && x == 100.0)
true

The above uses a NamedTuple to hold the conditioning variables, which allows us to perform some additional optimizations; in many cases, the above has zero runtime-overhead.

But we can also use a Dict, which offers more flexibility in the conditioning (see examples further below) but generally has worse performance than the NamedTuple approach:

julia> conditioned_model_dict = condition(model, Dict(@varname(x) => 100.0));

julia> m, x = conditioned_model_dict(); (m ≠ 1.0 && x == 100.0)
true

julia> # There's also an option using `|` by letting the right-hand side be a tuple
       # with elements of type `Pair{<:VarName}`, i.e. `vn => value` with `vn isa VarName`.
       conditioned_model_dict = model | (@varname(x) => 100.0, );

julia> m, x = conditioned_model_dict(); (m ≠ 1.0 && x == 100.0)
true

Condition only a part of a multivariate variable

Not only can be condition on multivariate random variables, but we can also use the standard mechanism of setting something to missing in the call to condition to only condition on a part of the variable.

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> conditioned_model = condition(model, m = [missing, 1.0]);

julia> # (✓) `m[1]` sampled while `m[2]` is fixed
       m = conditioned_model(); (m[1] ≠ 1.0 && m[2] == 1.0)
true

Intuitively one might also expect to be able to write model | (m[1] = 1.0, ). Unfortunately this is not supported as it has the potential of increasing compilation times but without offering any benefit with respect to runtime:

julia> # (×) `m[2]` is not set to 1.0.
       m = condition(model, var"m[2]" = 1.0)(); m[2] == 1.0
false

But you can do this if you use a Dict as the underlying storage instead:

julia> # Alternatives:
       # - `model | (@varname(m[2]) => 1.0,)`
       # - `condition(model, Dict(@varname(m[2] => 1.0)))`
       # (✓) `m[2]` is set to 1.0.
       m = condition(model, @varname(m[2]) => 1.0)(); (m[1] ≠ 1.0 && m[2] == 1.0)
true

Nested models

condition of course also supports the use of nested models through the use of to_submodel.

julia> @model demo_inner() = m ~ Normal()
demo_inner (generic function with 2 methods)

julia> @model function demo_outer()
           # By default, `to_submodel` prefixes the variables using the left-hand side of `~`.
           inner ~ to_submodel(demo_inner())
           return inner
       end
demo_outer (generic function with 2 methods)

julia> model = demo_outer();

julia> model() ≠ 1.0
true

julia> # To condition the variable inside `demo_inner` we need to refer to it as `inner.m`.
       conditioned_model = model | (@varname(inner.m) => 1.0, );

julia> conditioned_model()
1.0

julia> # However, it's not possible to condition `inner` directly.
       conditioned_model_fail = model | (inner = 1.0, );

julia> conditioned_model_fail()
ERROR: ArgumentError: `x ~ to_submodel(...)` is not supported when `x` is observed
[...]

conditioned(context::AbstractContext)

Return NamedTuple of values that are conditioned on under context`.

Note that this will recursively traverse the context stack and return a merged version of the condition values.

conditioned(model::Model)

Return the conditioned values in model.

Examples

julia> using Distributions

julia> using DynamicPPL: conditioned, contextualize, PrefixContext, ConditionContext

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

julia> m = demo();

julia> # Returns all the variables we have conditioned on + their values.
       conditioned(condition(m, x=100.0, m=1.0))
(x = 100.0, m = 1.0)

julia> # Nested ones also work.
       # (Note that `PrefixContext` also prefixes the variables of any
       # ConditionContext that is _inside_ it; because of this, the type of the
       # container has to be broadened to a `Dict`.)
       cm = condition(contextualize(m, PrefixContext(@varname(a), ConditionContext((m=1.0,)))), x=100.0);

julia> Set(keys(conditioned(cm))) == Set([@varname(a.m), @varname(x)])
true

julia> # Since we conditioned on `a.m`, it is not treated as a random variable.
       # However, `a.x` will still be a random variable.
       keys(VarInfo(cm))
1-element Vector{VarName}:
 a.x

julia> # We can also condition on `a.m` _outside_ of the PrefixContext:
       cm = condition(contextualize(m, PrefixContext(@varname(a))), (@varname(a.m) => 1.0));

julia> conditioned(cm)
Dict{VarName{:a, AbstractPPL.Property{:m, AbstractPPL.Iden}}, Float64} with 1 entry:
  a.m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(cm))
1-element Vector{VarName}:
 a.x

decondition(model::Model)
decondition(model::Model, variables...)

Return a Model for which variables... are not considered observations. If no variables are provided, then all variables currently considered observations will no longer be.

This is essentially the inverse of condition. This also means that it suffers from the same limitiations.

Note that currently we only support variables to take on explicit values provided to condition.

Examples

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> conditioned_model = condition(demo(), m = 1.0, x = 10.0);

julia> conditioned_model()
(m = 1.0, x = 10.0)

julia> # By specifying the `VarName` to `decondition`.
       model = decondition(conditioned_model, @varname(m));

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # When `NamedTuple` is used as the underlying, you can also provide
       # the symbol directly (though the `@varname` approach is preferable if
       # if the variable is known at compile-time).
       model = decondition(conditioned_model, :m);

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # `decondition` multiple at once:
       (m, x) = decondition(model, :m, :x)(); (m ≠ 1.0 && x ≠ 10.0)
true

julia> # `decondition` without any symbols will `decondition` all variables.
       (m, x) = decondition(model)(); (m ≠ 1.0 && x ≠ 10.0)
true

julia> # Usage of `Val` to perform `decondition` at compile-time if possible
       # is also supported.
       model = decondition(conditioned_model, Val{:m}());

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

Similarly when using a Dict:

julia> conditioned_model_dict = condition(demo(), @varname(m) => 1.0, @varname(x) => 10.0);

julia> conditioned_model_dict()
(m = 1.0, x = 10.0)

julia> deconditioned_model_dict = decondition(conditioned_model_dict, @varname(m));

julia> (m, x) = deconditioned_model_dict(); m ≠ 1.0 && x == 10.0
true

But, as mentioned, decondition is only supported for variables explicitly provided to condition earlier;

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> conditioned_model = condition(model, @varname(m) => [1.0, 2.0]);

julia> conditioned_model()
2-element Vector{Float64}:
 1.0
 2.0

julia> deconditioned_model = decondition(conditioned_model, @varname(m[1]));

julia> deconditioned_model()  # (×) `m[1]` is still conditioned
2-element Vector{Float64}:
 1.0
 2.0

julia> # (✓) this works though
       deconditioned_model_2 = deconditioned_model | (@varname(m[1]) => missing);

julia> m = deconditioned_model_2(); (m[1] ≠ 1.0 && m[2] == 2.0)
true

fix([context::AbstractContext,] values::NamedTuple)
fix([context::AbstractContext]; values...)

Return FixedContext with values and context if values is non-empty, otherwise return context which is DefaultContext by default.

See also: unfix

fix(model::Model; values...)
fix(model::Model, values::NamedTuple)

Return a Model which now treats the variables in values as fixed.

See also: unfix, fixed

Examples

Simple univariate model

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> m, x = model(); (m ≠ 1.0 && x ≠ 100.0)
true

julia> # Create a new instance which treats `x` as observed
       # with value `100.0`, and similarly for `m=1.0`.
       fixed_model = fix(model, x=100.0, m=1.0);

julia> m, x = fixed_model(); (m == 1.0 && x == 100.0)
true

julia> # Let's only fix on `x = 100.0`.
       fixed_model = fix(model, x = 100.0);

julia> m, x = fixed_model(); (m ≠ 1.0 && x == 100.0)
true

The above uses a NamedTuple to hold the fixed variables, which allows us to perform some additional optimizations; in many cases, the above has zero runtime-overhead.

But we can also use a Dict, which offers more flexibility in the fixing (see examples further below) but generally has worse performance than the NamedTuple approach:

julia> fixed_model_dict = fix(model, Dict(@varname(x) => 100.0));

julia> m, x = fixed_model_dict(); (m ≠ 1.0 && x == 100.0)
true

julia> # Alternative: pass `Pair{<:VarName}` as positional argument.
       fixed_model_dict = fix(model, @varname(x) => 100.0, );

julia> m, x = fixed_model_dict(); (m ≠ 1.0 && x == 100.0)
true

Fix only a part of a multivariate variable

We can not only fix multivariate random variables, but we can also use the standard mechanism of setting something to missing in the call to fix to only fix a part of the variable.

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> fixed_model = fix(model, m = [missing, 1.0]);

julia> # (✓) `m[1]` sampled while `m[2]` is fixed
       m = fixed_model(); (m[1] ≠ 1.0 && m[2] == 1.0)
true

Intuitively one might also expect to be able to write something like fix(model, var"m[1]" = 1.0, ). Unfortunately this is not supported as it has the potential of increasing compilation times but without offering any benefit with respect to runtime:

julia> # (×) `m[2]` is not set to 1.0.
       m = fix(model, var"m[2]" = 1.0)(); m[2] == 1.0
false

But you can do this if you use a Dict as the underlying storage instead:

julia> # Alternative: `fix(model, Dict(@varname(m[2] => 1.0)))`
       # (✓) `m[2]` is set to 1.0.
       m = fix(model, @varname(m[2]) => 1.0)(); (m[1] ≠ 1.0 && m[2] == 1.0)
true

Nested models

fix of course also supports the use of nested models through the use of to_submodel, similar to condition.

julia> @model demo_inner() = m ~ Normal()
demo_inner (generic function with 2 methods)

julia> @model function demo_outer()
           inner ~ to_submodel(demo_inner())
           return inner
       end
demo_outer (generic function with 2 methods)

julia> model = demo_outer();

julia> model() ≠ 1.0
true

julia> fixed_model = fix(model, (@varname(inner.m) => 1.0, ));

julia> fixed_model()
1.0

However, unlike condition, fix can also be used to fix the return-value of the submodel:

julia> fixed_model = fix(model, inner = 2.0,);

julia> fixed_model()
2.0

Difference from condition

A very similar functionality is also provided by condition. The only difference between fixing and conditioning is as follows:

• conditioned variables are considered to be observations, and are thus included in the computation logjoint and loglikelihood, but not in logprior.
• fixed variables are considered to be constant, and are thus not included in any log-probability computations.

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> model_fixed = fix(model, m = 1.0);

julia> model_conditioned = condition(model, m = 1.0);

julia> logjoint(model_fixed, (x=1.0,))
-0.9189385332046728

julia> # Different!
       logjoint(model_conditioned, (x=1.0,))
-2.3378770664093453

julia> # And the difference is the missing log-probability of `m`:
       logjoint(model_fixed, (x=1.0,)) + logpdf(Normal(), 1.0) == logjoint(model_conditioned, (x=1.0,))
true

fixed(context::AbstractContext)

Return the values that are fixed under context.

Note that this will recursively traverse the context stack and return a merged version of the fix values.

fixed(model::Model)

Return the fixed values in model.

Examples

julia> using Distributions

julia> using DynamicPPL: fixed, contextualize, PrefixContext

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

julia> m = demo();

julia> # Returns all the variables we have fixed on + their values.
       fixed(fix(m, x=100.0, m=1.0))
(x = 100.0, m = 1.0)

julia> # The rest of this is the same as the `condition` example above.
       fm = fix(contextualize(m, PrefixContext(@varname(a), fix(m=1.0))), x=100.0);

julia> Set(keys(fixed(fm))) == Set([@varname(a.m), @varname(x)])
true

julia> keys(VarInfo(fm))
1-element Vector{VarName}:
 a.x

julia> # We can also fix `a.m` _outside_ of the PrefixContext:
       fm = fix(contextualize(m, PrefixContext(@varname(a))), (@varname(a.m) => 1.0));

julia> fixed(fm)
Dict{VarName{:a, AbstractPPL.Property{:m, AbstractPPL.Iden}}, Float64} with 1 entry:
  a.m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(fm))
1-element Vector{VarName}:
 a.x

unfix(context::AbstractContext, syms...)

Return context but with syms no longer fixed.

Note that this recursively traverses contexts, unfixing all along the way.

See also: fix

unfix(model::Model)
unfix(model::Model, variables...)

Return a Model for which variables... are not considered fixed. If no variables are provided, then all variables currently considered fixed will no longer be.

This is essentially the inverse of fix. This also means that it suffers from the same limitiations.

Note that currently we only support variables to take on explicit values provided to fix.

Examples

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> fixed_model = fix(demo(), m = 1.0, x = 10.0);

julia> fixed_model()
(m = 1.0, x = 10.0)

julia> # By specifying the `VarName` to `unfix`.
       model = unfix(fixed_model, @varname(m));

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # When `NamedTuple` is used as the underlying, you can also provide
       # the symbol directly (though the `@varname` approach is preferable if
       # if the variable is known at compile-time).
       model = unfix(fixed_model, :m);

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # `unfix` multiple at once:
       (m, x) = unfix(model, :m, :x)(); (m ≠ 1.0 && x ≠ 10.0)
true

julia> # `unfix` without any symbols will `unfix` all variables.
       (m, x) = unfix(model)(); (m ≠ 1.0 && x ≠ 10.0)
true

julia> # Usage of `Val` to perform `unfix` at compile-time if possible
       # is also supported.
       model = unfix(fixed_model, Val{:m}());

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

Similarly when using a Dict:

julia> fixed_model_dict = fix(demo(), @varname(m) => 1.0, @varname(x) => 10.0);

julia> fixed_model_dict()
(m = 1.0, x = 10.0)

julia> unfixed_model_dict = unfix(fixed_model_dict, @varname(m));

julia> (m, x) = unfixed_model_dict(); m ≠ 1.0 && x == 10.0
true

But, as mentioned, unfix is only supported for variables explicitly provided to fix earlier:

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> fixed_model = fix(model, @varname(m) => [1.0, 2.0]);

julia> fixed_model()
2-element Vector{Float64}:
 1.0
 2.0

julia> unfixed_model = unfix(fixed_model, @varname(m[1]));

julia> unfixed_model()  # (×) `m[1]` is still fixed
2-element Vector{Float64}:
 1.0
 2.0

julia> # (✓) this works though
       unfixed_model_2 = fix(unfixed_model, @varname(m[1]) => missing);

julia> m = unfixed_model_2(); (m[1] ≠ 1.0 && m[2] == 2.0)
true

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

Sample from the posterior predictive distribution by executing model with parameters fixed to each sample in chain, and return the resulting Chains.

The model passed to predict is often different from the one used to generate chain. Typically, the model from which chain originated treats certain variables as observed (i.e., data points), while the model you pass to predict may mark these same variables as missing or unobserved. Calling predict then leverages the previously inferred parameter values to simulate what new, unobserved data might look like, given your posterior beliefs.

For each parameter configuration in chain:

1. All random variables present in chain are fixed to their sampled values.
2. Any variables not included in chain are sampled from their prior distributions.

If include_all is false, the returned Chains will contain only those variables that were not fixed by the samples in chain. This is useful when you want to sample only new variables from the posterior predictive distribution.

Examples

using AbstractMCMC, Distributions, DynamicPPL, Random

@model function linear_reg(x, y, σ = 0.1)
    β ~ Normal(0, 1)
    for i in eachindex(y)
        y[i] ~ Normal(β * x[i], σ)
    end
end

# Generate synthetic chain using known ground truth parameter
ground_truth_β = 2.0

# Create chain of samples from a normal distribution centered on ground truth
β_chain = MCMCChains.Chains(
    rand(Normal(ground_truth_β, 0.002), 1000), [:β,]
)

# Generate predictions for two test points
xs_test = [10.1, 10.2]

m_train = linear_reg(xs_test, fill(missing, length(xs_test)))

predictions = DynamicPPL.AbstractPPL.predict(
    Random.default_rng(), m_train, β_chain
)

ys_pred = vec(mean(Array(predictions); dims=1))

# Check if predictions match expected values within tolerance
(
    isapprox(ys_pred[1], ground_truth_β * xs_test[1], atol = 0.01),
    isapprox(ys_pred[2], ground_truth_β * xs_test[2], atol = 0.01)
)

# output

(true, true)

marginalize(
    model::DynamicPPL.Model,
    marginalized_varnames::AbstractVector{<:VarName};
    varinfo::DynamicPPL.AbstractVarInfo=link(VarInfo(model), model),
    getlogprob=DynamicPPL.getlogjoint,
    method::MarginalLogDensities.AbstractMarginalizer=MarginalLogDensities.LaplaceApprox();
    kwargs...,
)

Construct a MarginalLogDensities.MarginalLogDensity object that represents the marginal log-density of the given model, after marginalizing out the variables specified in varnames.

The resulting object can be called with a vector of parameter values to compute the marginal log-density.

Keyword arguments

• varinfo: The varinfo to use for the model. By default we use a linked VarInfo, meaning that the resulting log-density function accepts parameters that have been transformed to unconstrained space.

• getlogprob: A function which specifies which kind of marginal log-density to compute. Its default value is DynamicPPL.getlogjoint which returns the marginal log-joint probability.

• method: The marginalization method; defaults to a Laplace approximation. Please see the MarginalLogDensities.jl package for other options.

• Other keyword arguments are passed to the MarginalLogDensities.MarginalLogDensity constructor.

Example

julia> using DynamicPPL, Distributions, MarginalLogDensities

julia> @model function demo()
           x ~ Normal(1.0)
           y ~ Normal(2.0)
       end
demo (generic function with 2 methods)

julia> marginalized = marginalize(demo(), [:x]);

julia> # The resulting callable computes the marginal log-density of `y`.
       marginalized([1.0])
-1.4189385332046727

julia> logpdf(Normal(2.0), 1.0)
-1.4189385332046727

@model function f()
    x ~ Normal()
    y ~ truncated(Normal(); lower=x)
end

VarInfo(
    mld::MarginalLogDensities.MarginalLogDensity{<:LogDensityFunctionWrapper},
    unmarginalized_params::Union{AbstractVector,Nothing}=nothing
)

Retrieve the VarInfo object used in the marginalisation process.

If a Laplace approximation was used for the marginalisation, the values of the marginalized parameters are also set to their mode (note that this only happens if the mld object has been used to compute the marginal log-density at least once, so that the mode has been computed).

If a vector of unmarginalized_params is specified, the values for the corresponding parameters will also be updated in the returned VarInfo. This vector may be obtained e.g. by performing an optimization of the marginal log-density.

All other aspects of the VarInfo, such as link status, are preserved from the original VarInfo used in the marginalisation.

Example

julia> using DynamicPPL, Distributions, MarginalLogDensities

julia> @model function demo()
           x ~ Normal()
           y ~ Beta(2, 2)
       end
demo (generic function with 2 methods)

julia> # Note that by default `marginalize` uses a linked VarInfo.
       mld = marginalize(demo(), [@varname(x)]);

julia> using MarginalLogDensities: Optimization, OptimizationOptimJL

julia> # Find the mode of the marginal log-density of `y`, with an initial point of `y0`.
       y0 = 2.0; opt_problem = Optimization.OptimizationProblem(mld, [y0])
OptimizationProblem. In-place: true
u0: 1-element Vector{Float64}:
 2.0

julia> # This tells us the optimal (linked) value of `y` is around 0.
       opt_solution = Optimization.solve(opt_problem, OptimizationOptimJL.NelderMead())
retcode: Success
u: 1-element Vector{Float64}:
 4.88281250001733e-5

julia> # Get the VarInfo corresponding to the mode of `y`.
       vi = VarInfo(mld, opt_solution.u);

julia> # `x` is set to its mode (which for `Normal()` is zero).
       vi[@varname(x)]
0.0

julia> # `y` is set to the optimal value we found above.
       DynamicPPL.getindex_internal(vi, @varname(y))
1-element Vector{Float64}:
 4.88281250001733e-5

julia> # To obtain values in the original constrained space, we can either
       # use `getindex`:
       vi[@varname(y)]
0.5000122070312476

julia> # Or invlink the entire VarInfo object using the model:
       vi_unlinked = DynamicPPL.invlink(vi, demo()); vi_unlinked[:]
2-element Vector{Float64}:
 0.0
 0.5000122070312476

to_submodel(model::Model[, auto_prefix::Bool])

Return a model wrapper indicating that it is a sampleable model over the return-values.

This is mainly meant to be used on the right-hand side of a ~ operator to indicate that the model can be sampled from but not necessarily evaluated for its log density.

Arguments

• model::Model: the model to wrap.
• auto_prefix::Bool: whether to automatically prefix the variables in the model using the left-hand side of the ~ statement. Default: true.

Examples

Simple example

julia> @model function demo1(x)
           x ~ Normal()
           return 1 + abs(x)
       end;

julia> @model function demo2(x, y)
            a ~ to_submodel(demo1(x))
            return y ~ Uniform(0, a)
       end;

When we sample from the model demo2(missing, 0.4) random variable x will be sampled:

julia> vi = VarInfo(demo2(missing, 0.4));

julia> @varname(a.x) in keys(vi)
true

The variable a is not tracked. However, it will be assigned the return value of demo1, and can be used in subsequent lines of the model, as shown above.

julia> @varname(a) in keys(vi)
false

We can check that the log joint probability of the model accumulated in vi is correct:

julia> x = vi[@varname(a.x)];

julia> getlogjoint(vi) ≈ logpdf(Normal(), x) + logpdf(Uniform(0, 1 + abs(x)), 0.4)
true

Without automatic prefixing

As mentioned earlier, by default, the auto_prefix argument specifies whether to automatically prefix the variables in the submodel. If auto_prefix=false, then the variables in the submodel will not be prefixed.

julia> @model function demo1(x)
           x ~ Normal()
           return 1 + abs(x)
       end;

julia> @model function demo2_no_prefix(x, z)
            a ~ to_submodel(demo1(x), false)
            return z ~ Uniform(-a, 1)
       end;

julia> vi = VarInfo(demo2_no_prefix(missing, 0.4));

julia> @varname(x) in keys(vi)  # here we just use `x` instead of `a.x`
true

However, not using prefixing is generally not recommended as it can lead to variable name clashes unless one is careful. For example, if we're re-using the same model twice in a model, not using prefixing will lead to variable name clashes: However, one can manually prefix using the prefix(::Model, input):

julia> @model function demo2(x, y, z)
            a ~ to_submodel(prefix(demo1(x), :sub1), false)
            b ~ to_submodel(prefix(demo1(y), :sub2), false)
            return z ~ Uniform(-a, b)
       end;

julia> vi = VarInfo(demo2(missing, missing, 0.4));

julia> @varname(sub1.x) in keys(vi)
true

julia> @varname(sub2.x) in keys(vi)
true

Variables a and b are not tracked, but are assigned the return values of the respective calls to demo1:

julia> @varname(a) in keys(vi)
false

julia> @varname(b) in keys(vi)
false

We can check that the log joint probability of the model accumulated in vi is correct:

julia> sub1_x = vi[@varname(sub1.x)];

julia> sub2_x = vi[@varname(sub2.x)];

julia> logprior = logpdf(Normal(), sub1_x) + logpdf(Normal(), sub2_x);

julia> loglikelihood = logpdf(Uniform(-1 - abs(sub1_x), 1 + abs(sub2_x)), 0.4);

julia> getlogjoint(vi) ≈ logprior + loglikelihood
true

Usage as likelihood is illegal

Note that it is illegal to use a to_submodel model as a likelihood in another model:

julia> @model inner() = x ~ Normal()
inner (generic function with 2 methods)

julia> @model illegal_likelihood() = a ~ to_submodel(inner())
illegal_likelihood (generic function with 2 methods)

julia> model = illegal_likelihood() | (a = 1.0,);

julia> model()
ERROR: ArgumentError: `x ~ to_submodel(...)` is not supported when `x` is observed
[...]

prefix(ctx::AbstractContext, vn::VarName)

Apply the prefixes in the context ctx to the variable name vn.

prefix(model::Model, x::VarName)
prefix(model::Model, x::Val{sym})
prefix(model::Model, x::Any)

Return model but with all random variables prefixed by x, where x is either:

• a VarName (e.g. @varname(a)),
• a Val{sym} (e.g. Val(:a)), or
• for any other type, x is converted to a Symbol and then to a VarName. Note that this will introduce runtime overheads so is not recommended unless absolutely necessary.

Examples

julia> using DynamicPPL: prefix

julia> @model demo() = x ~ Dirac(1)
demo (generic function with 2 methods)

julia> rand(prefix(demo(), @varname(my_prefix)))
(var"my_prefix.x" = 1,)

julia> rand(prefix(demo(), Val(:my_prefix)))
(var"my_prefix.x" = 1,)

typed_identity(x)

Identity function, but with an overload for with_logabsdet_jacobian to ensure that it returns a sensible zero logjac.

The problem with plain old identity is that the default definition of with_logabsdet_jacobian for identity returns zero(eltype(x)): https://github.com/JuliaMath/ChangesOfVariables.jl/blob/d6a8115fc9b9419decbdb48e2c56ec9675b4c6a4/src/with_ladj.jl#L154

This is fine for most samples x, but if eltype(x) doesn't return a sensible type (e.g. if it's Any), then using identity will error with zero(Any). This can happen with, for example, ProductNamedTupleDistribution:

julia> using Distributions; d = product_distribution((a = Normal(), b = LKJCholesky(3, 0.5)));

julia> eltype(rand(d))
Any

The same problem precludes us from eventually broadening the scope of DynamicPPL.jl to support distributions with non-numeric samples.

Furthermore, in principle, the type of the log-probability should be separate from the type of the sample. Thus, instead of using zero(LogProbType), we should use the eltype of the LogJacobianAccumulator. There's no easy way to thread that through here, but if a way to do this is discovered, then typed_identity is what will allow us to obtain that custom behaviour.

@addlogprob!(ex)

Add a term to the log joint.

If ex evaluates to a NamedTuple with keys :loglikelihood and/or :logprior, the values are added to the log likelihood and log prior respectively.

If ex evaluates to a number it is added to the log likelihood.

Examples

julia> mylogjoint(x, μ) = (; loglikelihood=loglikelihood(Normal(μ, 1), x), logprior=1.0);

julia> @model function demo(x)
           μ ~ Normal()
           @addlogprob! mylogjoint(x, μ)
       end;

julia> x = [1.3, -2.1];

julia> loglikelihood(demo(x), (μ=0.2,)) ≈ mylogjoint(x, 0.2).loglikelihood
true

julia> logprior(demo(x), (μ=0.2,)) ≈ logpdf(Normal(), 0.2) + mylogjoint(x, 0.2).logprior
true

and to reject samples:

julia> @model function demo(x)
           m ~ MvNormal(zero(x), I)
           if dot(m, x) < 0
               @addlogprob! (; loglikelihood=-Inf)
               # Exit the model evaluation early
               return
           end
           x ~ MvNormal(m, I)
           return
       end;

julia> logjoint(demo([-2.1]), (m=[0.2],)) == -Inf
true

returned(model::Model, chain::MCMCChains.Chains)

Execute model for each of the samples in chain and return an array of the values returned by the model for each sample.

Examples

General

Often you might have additional quantities computed inside the model that you want to inspect, e.g.

@model function demo(x)
    # sample and observe
    θ ~ Prior()
    x ~ Likelihood()
    return interesting_quantity(θ, x)
end
m = demo(data)
chain = sample(m, alg, n)
# To inspect the `interesting_quantity(θ, x)` where `θ` is replaced by samples
# from the posterior/`chain`:
returned(m, chain) # <= results in a `Vector` of returned values
                               #    from `interesting_quantity(θ, x)`

Concrete (and simple)

julia> using Turing

julia> @model function demo(xs)
           s ~ InverseGamma(2, 3)
           m_shifted ~ Normal(10, √s)
           m = m_shifted - 10

           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end

           return (m, )
       end
demo (generic function with 1 method)

julia> model = demo(randn(10));

julia> chain = sample(model, MH(), 10);

julia> returned(model, chain)
10×1 Array{Tuple{Float64},2}:
 (2.1964758025119338,)
 (2.1964758025119338,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.043088571494005024,)
 (-0.16489786710222099,)
 (-0.16489786710222099,)

returned(model::Model, parameters...)

Initialise a model using the given parameters and return the model's return value. The parameters must be provided in a format that can be wrapped in an InitFromParams, i.e., InitFromParams(parameters..., nothing) must be a valid AbstractInitStrategy (where nothing is the fallback strategy to use if parameters are not provided).

As far as DynamicPPL is concerned, parameters can be either a singular NamedTuple or an AbstractDict{<:VarName}; however this method is left flexible to allow for other packages that wish to extend InitFromParams.

Example

julia> using DynamicPPL, Distributions

julia> @model function demo()
           m ~ Normal()
           return (mp1 = m + 1,)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> returned(model, (; m = 1.0))
(mp1 = 2.0,)

julia> returned(model, Dict{VarName,Float64}(@varname(m) => 2.0))
(mp1 = 3.0,)

DynamicPPL.pointwise_logdensities(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains,
    ::Type{Tout}=MCMCChains.Chains
    ::Val{whichlogprob}=Val(:both),
)

Runs model on each sample in chain, returning a new MCMCChains.Chains object where the log-density of each variable at each sample is stored (rather than its value).

whichlogprob specifies which log-probabilities to compute. It can be :both, :prior, or :likelihood.

You can pass Tout=OrderedDict to get the result as an OrderedDict{VarName, Matrix{Float64}} instead.

See also: DynamicPPL.pointwise_loglikelihoods, DynamicPPL.pointwise_prior_logdensities.

Examples

julia> using MCMCChains

julia> @model function demo(xs, y)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, √s)
           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end
           y ~ Normal(m, √s)
       end
demo (generic function with 2 methods)

julia> # Example observations.
       model = demo([1.0, 2.0, 3.0], [4.0]);

julia> # A chain with 3 iterations.
       chain = Chains(
           reshape(1.:6., 3, 2),
           [:s, :m];
           info=(varname_to_symbol=Dict(
               @varname(s) => :s,
               @varname(m) => :m,
           ),),
       );

julia> plds = pointwise_logdensities(model, chain)
Chains MCMC chain (3×6×1 Array{Float64, 3}):

Iterations        = 1:1:3
Number of chains  = 1
Samples per chain = 3
parameters        = s, m, xs[1], xs[2], xs[3], y
[...]

julia> plds[:s]
2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 1:1:3
    :chain, 1:1
And data, a 3×1 Matrix{Float64}:
 -0.8027754226637804
 -1.3822169643436162
 -2.0986122886681096

julia> # The above is the same as:
       logpdf.(InverseGamma(2, 3), chain[:s])
3×1 Matrix{Float64}:
 -0.8027754226637804
 -1.3822169643436162
 -2.0986122886681096

julia> # Alternatively: pldsdict = pointwiselogdensities(model, chain, OrderedDict) OrderedDict{VarName, Matrix{Float64}} with 6 entries: s => [-0.802775; -1.38222; -2.09861;;] m => [-8.91894; -7.51551; -7.46824;;] xs[1] => [-5.41894; -5.26551; -5.63491;;] xs[2] => [-2.91894; -3.51551; -4.13491;;] xs[3] => [-1.41894; -2.26551; -2.96824;;] y => [-0.918939; -1.51551; -2.13491;;]

DynamicPPL.pointwise_loglikelihoods(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains,
    ::Type{Tout}=MCMCChains.Chains
)

Compute the pointwise log-likelihoods of the model given the chain. This is the same as pointwise_logdensities(model, chain), but only including the likelihood terms.

See also: DynamicPPL.pointwise_logdensities, DynamicPPL.pointwise_prior_logdensities.

DynamicPPL.pointwise_prior_logdensities(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains
)

Compute the pointwise log-prior-densities of the model given the chain. This is the same as pointwise_logdensities(model, chain), but only including the prior terms.

See also: DynamicPPL.pointwise_logdensities, DynamicPPL.pointwise_loglikelihoods.

value_iterator_from_chain(model::Model, chain)
value_iterator_from_chain(varinfo::AbstractVarInfo, chain)

Return an iterator over the values in chain for each variable in model/varinfo.

Example

julia> using MCMCChains, DynamicPPL, Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end

           return s, m
       end
demo_model (generic function with 2 methods)

julia> model = demo_model([1.0, 2.0]);

julia> chain = Chains(rand(rng, 10, 2, 3), [:s, :m]);

julia> iter = value_iterator_from_chain(model, chain);

julia> first(iter)
OrderedDict{VarName, Any} with 2 entries:
  s => 0.580515
  m => 0.739328

julia> collect(iter)
10×3 Matrix{OrderedDict{VarName, Any}}:
 OrderedDict(s=>0.580515, m=>0.739328)  …  OrderedDict(s=>0.186047, m=>0.402423)
 OrderedDict(s=>0.191241, m=>0.627342)     OrderedDict(s=>0.776277, m=>0.166342)
 OrderedDict(s=>0.971133, m=>0.637584)     OrderedDict(s=>0.651655, m=>0.712044)
 OrderedDict(s=>0.74345, m=>0.110359)      OrderedDict(s=>0.469214, m=>0.104502)
 OrderedDict(s=>0.170969, m=>0.598514)     OrderedDict(s=>0.853546, m=>0.185399)
 OrderedDict(s=>0.704776, m=>0.322111)  …  OrderedDict(s=>0.638301, m=>0.853802)
 OrderedDict(s=>0.441044, m=>0.162285)     OrderedDict(s=>0.852959, m=>0.0956922)
 OrderedDict(s=>0.803972, m=>0.643369)     OrderedDict(s=>0.245049, m=>0.871985)
 OrderedDict(s=>0.772384, m=>0.646323)     OrderedDict(s=>0.906603, m=>0.385502)
 OrderedDict(s=>0.70882, m=>0.253105)      OrderedDict(s=>0.413222, m=>0.953288)

julia> # This can be used to `condition` a `Model`.
       conditioned_model = model | first(iter);

julia> conditioned_model()  # <= results in same values as the `first(iter)` above
(0.5805148626851955, 0.7393275279160691)

extract_priors([rng::Random.AbstractRNG, ]model::Model)

Extract the priors from a model.

This is done by sampling from the model and recording the distributions that are used to generate the samples.

Examples

Changing parameters

julia> using Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function model_dynamic_parameters()
           x ~ Normal(0, 1)
           y ~ Normal(x, 1)
       end;

julia> model = model_dynamic_parameters();

julia> extract_priors(rng, model)[@varname(y)]
Normal{Float64}(μ=-0.6702516921145671, σ=1.0)

julia> extract_priors(rng, model)[@varname(y)]
Normal{Float64}(μ=1.3736306979834252, σ=1.0)

Changing support

julia> using LinearAlgebra, Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function model_dynamic_support()
           n ~ Categorical(ones(10) ./ 10)
           x ~ MvNormal(zeros(n), I)
       end;

julia> model = model_dynamic_support();

julia> length(extract_priors(rng, model)[@varname(x)])
6

julia> length(extract_priors(rng, model)[@varname(x)])
9

extract_priors(model::Model, varinfo::AbstractVarInfo)

Extract the priors from a model.

This is done by evaluating the model at the values present in varinfo and recording the distributions that are present at each tilde statement.

values_as_in_model(model::Model, include_colon_eq::Bool, varinfo::AbstractVarInfo)

Get the values of varinfo as they would be seen in the model.

More specifically, this method attempts to extract the realization as seen in the model. For example, x[1] ~ truncated(Normal(); lower=0) will result in a realization that is compatible with truncated(Normal(); lower=0) – i.e. one where the value of x[1] is positive – regardless of whether varinfo is working in unconstrained space.

Hence this method is a "safe" way of obtaining realizations in constrained space at the cost of additional model evaluations.

Returns a VarNamedTuple.

Arguments

• model::Model: model to extract realizations from.
• include_colon_eq::Bool: whether to also include variables on the LHS of :=.
• varinfo::AbstractVarInfo: variable information to use for the extraction.

Examples

When VarInfo fails

The following demonstrates a common pitfall when working with VarInfo and constrained variables.

julia> using Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function model_changing_support()
           x ~ Bernoulli(0.5)
           y ~ x == 1 ? Uniform(0, 1) : Uniform(11, 12)
       end;

julia> model = model_changing_support();

julia> # Construct initial `VarInfo`.
       varinfo = VarInfo(rng, model);

julia> # Link it so it works in unconstrained space.
       varinfo_linked = DynamicPPL.link(varinfo, model);

julia> # Perform computations in unconstrained space, e.g. changing the values of `vals`.
       # Flip `x` so we hit the other support of `y`.
       vals = [!varinfo[@varname(x)], rand(rng)];

julia> # Update the `VarInfo` with the new values.
       varinfo_linked = DynamicPPL.unflatten!!(varinfo_linked, vals);

julia> # Determine the expected support of `y`.
       lb, ub = vals[1] == 1 ? (0, 1) : (11, 12)
(0, 1)

julia> # Approach 1: Convert back to constrained space using `invlink` and extract.
       varinfo_invlinked = DynamicPPL.invlink(varinfo_linked, model);

julia> lb ≤ first(varinfo_invlinked[@varname(y)]) ≤ ub
true

julia> # Approach 2: Extract realizations using `values_as_in_model`.
       lb ≤ values_as_in_model(model, true, varinfo_linked)[@varname(y)] ≤ ub
true

A named distribution that carries the name of the random variable with it.

run_ad(
    model::Model,
    adtype::ADTypes.AbstractADType;
    test::Union{AbstractADCorrectnessTestSetting,Bool}=WithBackend(),
    benchmark=false,
    atol::AbstractFloat=1e-8,
    rtol::AbstractFloat=sqrt(eps()),
    getlogdensity::Function=getlogjoint_internal,
    rng::Random.AbstractRNG=Random.default_rng(),
    varinfo::AbstractVarInfo=link(VarInfo(model), model),
    params::Union{Nothing,Vector{<:AbstractFloat}}=nothing,
    verbose=true,
)::ADResult

Description

Test the correctness and/or benchmark the AD backend adtype for the model model.

Whether to test and benchmark is controlled by the test and benchmark keyword arguments. By default, test is true and benchmark is false.

Note that to run AD successfully you will need to import the AD backend itself. For example, to test with AutoReverseDiff() you will need to run import ReverseDiff.

Arguments

There are two positional arguments, which absolutely must be provided:

1. model - The model being tested.
2. adtype - The AD backend being tested.

Everything else is optional, and can be categorised into several groups:

1. How to specify the VarInfo.

   DynamicPPL contains several different types of VarInfo objects which change the way model evaluation occurs. If you want to use a specific type of VarInfo, pass it as the varinfo argument. Otherwise, it will default to using a linked TypedVarInfo generated from the model. Here, linked means that the parameters in the VarInfo have been transformed to unconstrained Euclidean space if they aren't already in that space.

2. How to specify the parameters.

   For maximum control over this, generate a vector of parameters yourself and pass this as the params argument. If you don't specify this, it will be taken from the contents of the VarInfo.

   Note that if the VarInfo is not specified (and thus automatically generated) the parameters in it will have been sampled from the prior of the model. If you want to seed the parameter generation for the VarInfo, you can pass the rng keyword argument, which will then be used to create the VarInfo.

   Finally, note that these only reflect the parameters used for evaluating the gradient. If you also want to control the parameters used for preparing the gradient, then you need to manually set these parameters in the VarInfo object, for example using vi = DynamicPPL.unflatten!!(vi, prep_params). You could then evaluate the gradient at a different set of parameters using the params keyword argument.

3. Which type of logp is being calculated.

   By default, run_ad evaluates the 'internal log joint density' of the model, i.e., the log joint density in the unconstrained space. Thus, for example, in

   @model f() = x ~ LogNormal()

   the internal log joint density is logpdf(Normal(), log(x)). This is the relevant log density for e.g. Hamiltonian Monte Carlo samplers and is therefore the most useful to test.

   If you want the log joint density in the original model parameterisation, you can use getlogjoint. Likewise, if you want only the prior or likelihood, you can use getlogprior or getloglikelihood, respectively.

4. How to specify the results to compare against.

   Once logp and its gradient has been calculated with the specified adtype, it can optionally be tested for correctness. The exact way this is tested is specified in the test parameter.

   There are several options for this:

   • You can explicitly specify the correct value using WithExpectedResult().
   • You can compare against the result obtained with a different AD backend using WithBackend(adtype).
   • You can disable testing by passing NoTest().
   • The default is to compare against the result obtained with ForwardDiff, i.e. WithBackend(AutoForwardDiff()).
   • test=false and test=true are synonyms for NoTest() and WithBackend(AutoForwardDiff()), respectively.

5. How to specify the tolerances. (Only if testing is enabled.)

   Both absolute and relative tolerances can be specified using the atol and rtol keyword arguments respectively. The behaviour of these is similar to isapprox(), i.e. the value and gradient are considered correct if either atol or rtol is satisfied. The default values are 100*eps() for atol and sqrt(eps()) for rtol.

   For the most part, it is the rtol check that is more meaningful, because we cannot know the magnitude of logp and its gradient a priori. The atol value is supplied to handle the case where gradients are equal to zero.

6. Whether to benchmark.

   By default, benchmarking is disabled. To enable it, set benchmark=true. When enabled, the time taken to evaluate logp as well as its gradient is measured using Chairmarks.jl, and the ADResult object returned will contain grad_time and primal_time fields with the median times (in seconds).

7. Whether to output extra logging information.

   By default, this function prints messages when it runs. To silence it, set verbose=false.

Returns / Throws

Returns an ADResult object, which contains the results of the test and/or benchmark.

If test is true and the AD backend returns an incorrect value or gradient, an ADIncorrectException is thrown. If a different error occurs, it will be thrown as-is.

AbstractADCorrectnessTestSetting

Different ways of testing the correctness of an AD backend.

WithBackend(adtype::AbstractADType=AutoForwardDiff()) <: AbstractADCorrectnessTestSetting

Test correctness by comparing it against the result obtained with adtype.

adtype defaults to ForwardDiff.jl, since it's the default AD backend used in Turing.jl.

WithExpectedResult(
    value::T,
    grad::AbstractVector{T}
) where {T <: AbstractFloat}
<: AbstractADCorrectnessTestSetting

Test correctness by comparing it against a known result (e.g. one obtained analytically, or one obtained with a different backend previously). Both the value of the primal (i.e. the log-density) as well as its gradient must be supplied.

NoTest() <: AbstractADCorrectnessTestSetting

Disable correctness testing.

ADResult{Tparams<:AbstractFloat,Tresult<:AbstractFloat,Ttol<:AbstractFloat}

Data structure to store the results of the AD correctness test.

The type parameter Tparams is the numeric type of the parameters passed in; Tresult is the type of the value and the gradient; and Ttol is the type of the absolute and relative tolerances used for correctness testing.

Fields

• model::Model: The DynamicPPL model that was tested

• getlogdensity::Function: The function used to extract the log density from the model

• varinfo::AbstractVarInfo: The VarInfo that was used

• params::Vector{Tparams} where Tparams<:AbstractFloat: The values at which the model was evaluated

• adtype::ADTypes.AbstractADType: The AD backend that was tested

• atol::AbstractFloat: Absolute tolerance used for correctness test

• rtol::AbstractFloat: Relative tolerance used for correctness test

• value_expected::Union{Nothing, Tresult} where Tresult<:AbstractFloat: The expected value of logp

• grad_expected::Union{Nothing, Vector{Tresult}} where Tresult<:AbstractFloat: The expected gradient of logp

• value_actual::AbstractFloat: The value of logp (calculated using adtype)

• grad_actual::Vector{Tresult} where Tresult<:AbstractFloat: The gradient of logp (calculated using adtype)

• grad_time::Union{Nothing, Tresult} where Tresult<:AbstractFloat: If benchmarking was requested, the time taken by the AD backend to evaluate the gradient of logp (in seconds)

• primal_time::Union{Nothing, Tresult} where Tresult<:AbstractFloat: If benchmarking was requested, the time taken by the AD backend to evaluate logp (in seconds)

ADIncorrectException{T<:AbstractFloat}

Exception thrown when an AD backend returns an incorrect value or gradient.

The type parameter T is the numeric type of the value and gradient.

Fields

• value_expected::AbstractFloat

• value_actual::AbstractFloat

• grad_expected::Vector{T} where T<:AbstractFloat

• grad_actual::Vector{T} where T<:AbstractFloat

• atol::AbstractFloat

• rtol::AbstractFloat

A collection of models corresponding to the posterior distribution defined by the generative process

s ~ InverseGamma(2, 3)
m ~ Normal(0, √s)
1.5 ~ Normal(m, √s)
2.0 ~ Normal(m, √s)

or by

s[1] ~ InverseGamma(2, 3)
s[2] ~ InverseGamma(2, 3)
m[1] ~ Normal(0, √s)
m[2] ~ Normal(0, √s)
1.5 ~ Normal(m[1], √s[1])
2.0 ~ Normal(m[2], √s[2])

These are examples of a Normal-InverseGamma conjugate prior with Normal likelihood, for which the posterior is known in closed form.

In particular, for the univariate model (the former one):

mean(s) == 49 / 24
mean(m) == 7 / 6

And for the multivariate one (the latter one):

mean(s[1]) == 19 / 8
mean(m[1]) == 3 / 4
mean(s[2]) == 8 / 3
mean(m[2]) == 1

A tuple of all models defined in DynamicPPL.TestUtils.

logprior_true(model, args...)

Return the logprior of model for args.

This should generally be implemented by hand for every specific model.

See also: logjoint_true, loglikelihood_true.

loglikelihood_true(model, args...)

Return the loglikelihood of model for args.

This should generally be implemented by hand for every specific model.

See also: logjoint_true, logprior_true.

logjoint_true(model, args...)

Return the logjoint of model for args.

Defaults to logprior_true(model, args...) + loglikelihood_true(model, args..).

This should generally be implemented by hand for every specific model so that the returned value can be used as a ground-truth for testing things like:

1. Validity of evaluation of model using a particular implementation of AbstractVarInfo.
2. Validity of a sampler when combined with DynamicPPL by running the sampler twice: once targeting ground-truth functions, e.g. logjoint_true, and once targeting model.

And more.

See also: logprior_true, loglikelihood_true.

logprior_true_with_logabsdet_jacobian(model::Model, args...)

Return a tuple (args_unconstrained, logprior_unconstrained) of model for args....

Unlike logprior_true, the returned logprior computation includes the log-absdet-jacobian adjustment, thus computing logprior for the unconstrained variables.

Note that args are assumed be in the support of model, while args_unconstrained will be unconstrained.

See also: logprior_true.

logjoint_true_with_logabsdet_jacobian(model::Model, args...)

Return a tuple (args_unconstrained, logjoint) of model for args.

Unlike logjoint_true, the returned logjoint computation includes the log-absdet-jacobian adjustment, thus computing logjoint for the unconstrained variables.

Note that args are assumed be in the support of model, while args_unconstrained will be unconstrained.

This should generally not be implemented directly, instead one should implement logprior_true_with_logabsdet_jacobian for a given model.

See also: logjoint_true, logprior_true_with_logabsdet_jacobian.

varnames(model::Model)

Return the VarNames defined in model, as a Vector.

posterior_mean(model::Model)

Return a NamedTuple compatible with varnames(model) where the values represent the posterior mean under model.

"Compatible" means that a varname from varnames(model) can be used to extract the corresponding value using e.g. AbstractPPL.getvalue(posterior_mean(model), varname).

setup_varinfos(model::Model, example_values::NamedTuple, varnames; include_threadsafe::Bool=false)

Return a tuple of instances for different implementations of AbstractVarInfo with each vi, supposedly, satisfying vi[vn] == get(example_values, vn) for vn in varnames.

If include_threadsafe is true, then the returned tuple will also include thread-safe versions of the varinfo instances.

update_values!!(vi::AbstractVarInfo, vals::NamedTuple, vns)

Return instance similar to vi but with vns set to values from vals.

test_values(vi::AbstractVarInfo, vals::NamedTuple, vns)

Test that vi[vn] corresponds to the correct value in vals for every vn in vns.

check_model(model::Model, varinfo::AbstractVarInfo; error_on_failure=false)

Check that model is valid, warning about any potential issues (or erroring if error_on_failure is true).

Returns

• issuccess::Bool: Whether the model check succeeded.

check_model_and_trace(model::Model, varinfo::AbstractVarInfo; error_on_failure=false)

Check that evaluating model with the given varinfo is valid, warning about any potential issues.

This will check the model for the following issues:

1. Repeated usage of the same varname in a model.
2. NaN on the left-hand side of observe statements.

Arguments

• model::Model: The model to check.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model.

Keyword Argument

• error_on_failure::Bool: Whether to throw an error if the model check fails. Default: false.

Returns

• issuccess::Bool: Whether the model check succeeded.
• trace::Vector{Stmt}: The trace of statements executed during the model check.

Examples

Correct model

julia> using StableRNGs

julia> rng = StableRNG(42);

julia> @model demo_correct() = x ~ Normal()
demo_correct (generic function with 2 methods)

julia> model = demo_correct(); varinfo = VarInfo(rng, model);

julia> issuccess, trace = check_model_and_trace(model, varinfo);

julia> issuccess
true

julia> print(trace)
 assume: x ~ Normal{Float64}(μ=0.0, σ=1.0) ⟼ -0.670252

julia> cond_model = model | (x = 1.0,);

julia> issuccess, trace = check_model_and_trace(cond_model, VarInfo(cond_model));
┌ Warning: The model does not contain any parameters.
└ @ DynamicPPL.DebugUtils DynamicPPL.jl/src/debug_utils.jl:342

julia> issuccess
true

julia> print(trace)
 observe: x (= 1.0) ~ Normal{Float64}(μ=0.0, σ=1.0)

Incorrect model

julia> @model function demo_incorrect()
           # (×) Sampling `x` twice will lead to incorrect log-probabilities!
           x ~ Normal()
           x ~ Exponential()
       end
demo_incorrect (generic function with 2 methods)

julia> # Notice that VarInfo(model_incorrect) evaluates the model, but doesn't actually
       # alert us to the issue of `x` being sampled twice.
       model = demo_incorrect(); varinfo = VarInfo(model);

julia> issuccess, trace = check_model_and_trace(model, varinfo; error_on_failure=true);
ERROR: varname x used multiple times in model

has_static_constraints([rng, ]model::Model; num_evals=5, error_on_failure=false)

Return true if the model has static constraints, false otherwise.

Note that this is a heuristic check based on sampling from the model multiple times and checking if the model is consistent across runs.

Arguments

• rng::Random.AbstractRNG: The random number generator to use when evaluating the model.
• model::Model: The model to check.

Keyword Arguments

• num_evals::Int: The number of evaluations to perform. Default: 5.
• error_on_failure::Bool: Whether to throw an error if any of the num_evals model checks fail. Default: false.

model_warntype(model[, varinfo]; optimize=true)

Check the type stability of the model's evaluator, warning about any potential issues.

This simply calls @code_warntype on the model's evaluator, filling in internal arguments where needed.

Arguments

• model::Model: The model to check.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Keyword Arguments

• optimize::Bool: Whether to generate optimized code. Default: false.

model_typed(model[, varinfo]; optimize=true)

Return the type inference for the model's evaluator.

This simply calls @code_typed on the model's evaluator, filling in internal arguments where needed.

Arguments

• model::Model: The model to check.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Keyword Arguments

• optimize::Bool: Whether to generate optimized code. Default: true.

gen_evaluator_call_with_types(model[, varinfo])

Generate the evaluator call and the types of the arguments.

Arguments

• model::Model: The model whose evaluator is of interest.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Returns

A 2-tuple with the following elements:

• f: This is either model.f or Core.kwcall, depending on whether the model has keyword arguments.
• argtypes::Type{<:Tuple}: The types of the arguments for the evaluator.

AbstractVarInfo

Abstract supertype for data structures that capture random variables when executing a probabilistic model and accumulate log densities such as the log likelihood or the log joint probability of the model.

See also: VarInfo

VarInfo{Linked,T<:VarNamedTuple,Accs<:AccumulatorTuple} <: AbstractVarInfo

The default implementation of AbstractVarInfo, storing variable values and accumulators.

The Linked type parameter is either true or false to mark that all variables in this VarInfo are linked, or nothing to indicate that some variables may be linked and some not, and a runtime check is needed.

VarInfo is quite a thin wrapper around a VarNamedTuple storing the variable values, and a tuple of accumulators. The only really noteworthy thing about it is that it stores the values of variables vectorised as instances of AbstractTransformedValue. That is, it stores each value as a special vector with a flag indicating whether it is just a vectorised value (VectorValue), or whether it is also linked (LinkedVectorValue). It also stores the size of the actual post-transformation value. These are all accessible via AbstractTransformedValue.

Note that setindex!! and getindex on VarInfo take and return values in the support of the original distribution. To get access to the internal vectorised values, use getindex_internal, setindex_internal!!, and unflatten!!.

There's also a VarInfo-specific function setindex_with_dist!!, which sets a variable's value with a transformation based on the statistical distribution this value is a sample for.

For more details on the internal storage, see documentation of AbstractTransformedValue and VarNamedTuple.

Fields

• values::VarNamedTuple

• accs::DynamicPPL.AccumulatorTuple

setindex_with_dist!!(vi::VarInfo, val, dist::Distribution, vn::VarName)

Set the value of vn in vi to val, applying a transformation based on dist.

val is taken to be the actual value of the variable, and is transformed into the internal (vectorised) representation using a transformation based on dist. If the variable is currently linked in vi, or doesn't exist in vi but all other variables in vi are linked, the linking transformation is used; otherwise, the standard vector transformation is used.

Returns three things:

• the modified vi,
• the log absolute determinant of the Jacobian of the transformation applied,
• the AbstractTransformedValue used to store the value.

is_transformed(vi::AbstractVarInfo[, vns::Union{VarName, AbstractVector{<:Varname}}])

Return true if vi is working in unconstrained space, and false if vi is assuming realizations to be in support of the corresponding distributions.

If vns is provided, then only check if this/these varname(s) are transformed.

set_transformed!!(vi::AbstractVarInfo, trans::Bool[, vn::VarName])

Return vi with is_transformed(vi, vn) evaluating to true.

If vn is not specified, then is_transformed(vi) evaluates to true for all variables.

VarNamedTuple{names,Values}

A NamedTuple-like structure with VarName keys.

VarNamedTuple is a data structure for storing arbitrary data, keyed by VarNames, in an efficient and type stable manner. It is mainly used through getindex, setindex!!, templated_setindex!!, and haskey, all of which only accept VarNames as keys. Other notable methods are merge and subset.

VarNamedTuple has an ordering to its elements, and two VarNamedTuples with the same keys and values but in different orders are considered different for equality and hashing. Iterations such as keys and values respect this ordering. The ordering is dependent on the order in which elements were inserted into the VarNamedTuple, though isn't always equal to it. More specifically

• Any new keys that have a joint parent VarName with an existing key are inserted after that key. For instance, if one first inserts, in order, @varname(a.x), @varname(b), and @varname(a.y), the resulting order will be (@varname(a.x), @varname(a.y), @varname(b)).
• Index keys, like@varname(a[3])or@varname(b[2,3,4:5]), are always iterated in the same order anArraywith the same indices would be iterated. For instance, if one first inserts, in order,@varname(a[2]),@varname(b), and@varname(a[1]), the resulting order will be(@varname(a[1]), @varname(a[2]), @varname(b))`.

Otherwise insertion order is respected.

setindex!! and getindex on VarNamedTuple are type stable as long as one does not store heterogeneous data under different indices of the same symbol. That is, if either

• one sets a[1] and a[2] to be of different types, or
• if a[1] and a[2] both exist, one sets a[1].b without setting a[2].b,

then getting values for a[1] or a[2] will not be type stable.

vnt_size(x)

Get the size of an object x for use in VarNamedTuple and PartialArray.

By default, this falls back onto Base.size, but can be overloaded for custom types. This notion of type is used to determine whether a value can be set into a PartialArray as a block, see the docstring of PartialArray and ArrayLikeBlock for details.

apply!!(func, vnt::VarNamedTuple, name::VarName)

Apply func to the subdata at name in vnt, and set the result back at name.

Like map_values!!, but only for a single VarName.

julia> using DynamicPPL: VarNamedTuple, setindex!!

julia> using DynamicPPL.VarNamedTuples: apply!!

julia> vnt = VarNamedTuple()
VarNamedTuple()

julia> vnt = setindex!!(vnt, [1, 2, 3], @varname(a))
VarNamedTuple
└─ a => [1, 2, 3]

julia> apply!!(x -> x .+ 1, vnt, @varname(a))
VarNamedTuple
└─ a => [2, 3, 4]

map_pairs!!(func, vnt::VarNamedTuple)

Apply func to all key => value pairs of vnt, in place if possible.

func should accept a pair of VarName and value, and return the new value to be set.

map_values!!(func, vnt::VarNamedTuple)

Apply func to elements of vnt, in place if possible.

PartialArray{
    ElType,
    num_dims,
    D<:AbstractArray{ElType,num_dims},
    M<:AbstractArray{Bool,num_dims}
}

An array-like like structure that may only have some of its elements defined.

One can set values in a PartialArray either element-by-element, or with ranges like arr[1:3,2] = [5,10,15]. When setting values over a range of indices, the value being set must either be an AbstractArray or otherwise something for which vnt_size(value) or Base.size(value) (which vnt_size falls back onto) is defined, and the size matches the range. If the value is an AbstractArray, the elements are copied individually, but if it is not, the value is stored as a block, that takes up the whole range, e.g. [1:3,2], but is only a single object. Getting such a block-value must be done with the exact same range of indices, otherwise an error is thrown.

If the element type of a PartialArray is not concrete, any call to setindex!! will check if, after the new value has been set, the element type can be made more concrete. If so, a new PartialArray with a more concrete element type is returned. Thus the element type of any PartialArray should always be as concrete as is allowed by the elements in it.

The internal implementation of an PartialArray consists of two arrays: one holding the data and the other one being a boolean mask indicating which elements are defined. These internal arrays may need resizing when new elements are set that have index ranges larger than the current internal arrays. To avoid resizing too often, the internal arrays are resized in exponentially increasing steps. This means that most setindex!! calls are very fast, but some may incur substantial overhead due to resizing and copying data. It also means that the largest index set so far determines the memory usage of the PartialArray. PartialArrays are thus well-suited when most values in it will eventually be set. If only a few scattered values are set, a structure like SparseArray may be more appropriate.

DynamicPPL.VarNamedTuples.templated_setindex!!(vnt, value, vn, template; allow_new=Val(true))

Assign value to the location in vnt specified by vn.

The argument template must be provided in order to guide the creation of PartialArrays, as well as to concretise any dynamic indices in vn. It must be an object that has the shape of the top-level symbol in vn. For example:

vnt = VarNamedTuple()
templated_setindex!!(vnt, 10, @varname(x[1]), rand(2, 2))

Here, rand(2, 2) is the template for the top-level symbol x, which tells setindex!! that x should be a PartialArray that is backed by a matrix.

The actual data inside template is not needed, and template is never mutated by this call.

NoTemplate

A singleton struct representing the fact that there is no template for a top-level variable. When NoTemplate is used, several things happen:

• If you attempt to call make_leaf with an Index optic, this will error.
• When recursing into substructures, NoTemplate will be propagated.

Collectively this means that you can only set values for variables that only have Property optics, unless a template is provided.

It might seem more idiomatic to use nothing or missing for this. However, this causes a bug with BangBang.setindex!!: https://github.com/JuliaFolds2/BangBang.jl/issues/43 so we use a dedicated struct instead.

SkipTemplate{N}(value)

A struct representing the fact that value is the template for the variable N levels down from the top-level variable. In other words, SkipTemplate{0}(value) is equivalent to just value, and SkipTemplate{1}(value) means that value is the template for a when setting the variable @varname(x.a).

AbstractAccumulator

An abstract type for accumulators.

An accumulator is an object that may change its value at every tildeassume!! or tildeobserve!! call based on the random variable in question. The obvious examples of accumulators are the log prior and log likelihood. Other examples might be a variable that counts the number of observations in a trace, or a list of the names of random variables seen so far.

An accumulator type T <: AbstractAccumulator must implement the following methods:

• accumulator_name(acc::T) or accumulator_name(::Type{T})
• accumulate_observe!!(acc::T, dist, val, vn)
• accumulate_assume!!(acc::T, val, tval, logjac, vn, dist, template)
• reset(acc::T)
• Base.copy(acc::T)

In these functions:

• val is the new value of the random variable sampled from a distribution (always in the original unlinked space), or the value on the left-hand side of an observe statement.
• tval is the original AbstractTransformedValue that was obtained from the initialisation strategy. This is passed through unchanged to accumulate_assume!! since it can be reused for some accumulators (e.g. when storing linked values, if the linked value was already provided, it is faster to reuse it than to re-link val).
• dist is the distribution on the RHS of the tilde statement.
• vn is the VarName that is on the left-hand side of the tilde-statement. If the tilde-statement is a literal observation like 0.0 ~ Normal(), then vn is nothing.
• logjac is the log determinant of the Jacobian of the link transformation, if the variable is stored as a linked value in the VarInfo. If the variable is stored in its original, unlinked form, then logjac is zero.
• template is a value that conveys the shape of the top-level symbol in vn, and is used specifically for accumulators that carry VarNamedTuples.

To be able to work with multi-threading, it should also implement:

• split(acc::T)
• combine(acc::T, acc2::T)

See the documentation for each of these functions for more details.

LogPriorAccumulator{T<:Real} <: LogProbAccumulator{T}

An accumulator that tracks the cumulative log prior during model execution.

Note that the log prior stored in here is always calculated based on unlinked parameters, i.e., the value of logp is independent of whether tha VarInfo is linked or not.

Fields

• logp::Real: the scalar log prior value

LogJacobianAccumulator{T<:Real} <: LogProbAccumulator{T}

An accumulator that tracks the cumulative log Jacobian (technically, log(abs(det(J)))) during model execution. Specifically, J refers to the Jacobian of the link transform, i.e., from the space of the original distribution to unconstrained space.

getlogjoint_internal(vi) = getlogjoint(vi) - getlogjac(vi)

Fields

• logjac::Real: the logabsdet of the link transform Jacobian

LogLikelihoodAccumulator{T<:Real} <: LogProbAccumulator{T}

An accumulator that tracks the cumulative log likelihood during model execution.

Fields

• logp::Real: the scalar log likelihood value

getlogp(vi::AbstractVarInfo)

Return a NamedTuple of the log prior, log Jacobian, and log likelihood probabilities.

The keys are called logprior, logjac, and loglikelihood. If any of them are not present in vi an error will be thrown.

setlogp!!(vi::AbstractVarInfo, logp::NamedTuple)

Set both the log prior and the log likelihood probabilities in vi.

logp should have fields logprior and loglikelihood and no other fields.

See also: setlogprior!!, setloglikelihood!!, getlogp.

acclogp!!(vi::AbstractVarInfo, logp::NamedTuple; ignore_missing_accumulator::Bool=false)

Add to both the log prior and the log likelihood probabilities in vi.

logp should have fields logprior and/or loglikelihood, and no other fields.

By default if the necessary accumulators are not in vi an error is thrown. If ignore_missing_accumulator is set to true then this is silently ignored instead.

getlogjoint(vi::AbstractVarInfo)

Return the log of the joint probability of the observed data and parameters in vi.

See also: getlogprior, getloglikelihood.

getlogjoint_internal(vi::AbstractVarInfo)

Return the log of the joint probability of the observed data and parameters as they are stored internally in vi, including the log-Jacobian for any linked parameters.

In general, we have that:

getlogjoint_internal(vi) == getlogjoint(vi) - getlogjac(vi)

getlogjac(vi::AbstractVarInfo)

Return the accumulated log-Jacobian term for any linked parameters in vi. The Jacobian here is taken with respect to the forward (link) transform.

See also: setlogjac!!.

setlogjac!!(vi::AbstractVarInfo, logjac)

Set the accumulated log-Jacobian term for any linked parameters in vi. The Jacobian here is taken with respect to the forward (link) transform.

See also: getlogjac, acclogjac!!.

acclogjac!!(vi::AbstractVarInfo, logjac)

Add logjac to the value of the log Jacobian in vi.

See also: getlogjac, setlogjac!!.

getlogprior(vi::AbstractVarInfo)

Return the log of the prior probability of the parameters in vi.

See also: getlogjoint, getloglikelihood, setlogprior!!.

getlogprior_internal(vi::AbstractVarInfo)

Return the log of the prior probability of the parameters as stored internally in vi. This includes the log-Jacobian for any linked parameters.

In general, we have that:

getlogprior_internal(vi) == getlogprior(vi) - getlogjac(vi)

setlogprior!!(vi::AbstractVarInfo, logp)

Set the log of the prior probability of the parameters sampled in vi to logp.

See also: setloglikelihood!!, setlogp!!, getlogprior.

acclogprior!!(vi::AbstractVarInfo, logp)

Add logp to the value of the log of the prior probability in vi.

See also: accloglikelihood!!, acclogp!!, getlogprior, setlogprior!!.

getloglikelihood(vi::AbstractVarInfo)

Return the log of the likelihood probability of the observed data in vi.

See also: getlogjoint, getlogprior, setloglikelihood!!.

setloglikelihood!!(vi::AbstractVarInfo, logp)

Set the log of the likelihood probability of the observed data sampled in vi to logp.

See also: setlogprior!!, setlogp!!, getloglikelihood.

accloglikelihood!!(vi::AbstractVarInfo, logp)

Add logp to the value of the log of the likelihood in vi.

See also: accloglikelihood!!, acclogp!!, getloglikelihood, setloglikelihood!!.

keys(vi::AbstractVarInfo)

Return an iterator over all vns in vi.

getindex(vi::AbstractVarInfo, vn::VarName[, dist::Distribution])
getindex(vi::AbstractVarInfo, vns::Vector{<:VarName}[, dist::Distribution])

Return the current value(s) of vn (vns) in vi in the support of its (their) distribution(s).

If dist is specified, the value(s) will be massaged into the representation expected by dist.

empty!!(vi::AbstractVarInfo)

Empty vi of variables and reset all accumulators.

This is useful when using a sampling algorithm that assumes an empty vi, e.g. SMC.

isempty(vi::AbstractVarInfo)

Return true if vi is empty and false otherwise.

getindex_internal(vi::AbstractVarInfo, vn::VarName)
getindex_internal(vi::AbstractVarInfo, vns::Vector{<:VarName})
getindex_internal(vi::AbstractVarInfo, ::Colon)

Return the internal value of the varname vn, varnames vns, or all varnames in vi respectively. The internal value is the value of the variables that is stored in the varinfo object; this may be the actual realisation of the random variable (i.e. the value sampled from the distribution), or it may have been transformed to Euclidean space, depending on whether the varinfo was linked.

See https://turinglang.org/docs/developers/transforms/dynamicppl/ for more information on how transformed variables are stored in DynamicPPL.

See also: getindex(vi::AbstractVarInfo, vn::VarName, dist::Distribution)

setindex_internal!!(vi::VarInfo, val, vn::VarName)

Set the internal (vectorised) value of variable vn in vi to val.

This does not change the transformation or linked status of the variable.

values_as(varinfo[, Type])

Return the values/realizations in varinfo as Type, if implemented.

If no Type is provided, return values as stored in varinfo.

Examples

julia> # Just use an example model to construct the `VarInfo` because we're lazy.
       vi = DynamicPPL.VarInfo(DynamicPPL.TestUtils.demo_assume_dot_observe());

julia> vi = DynamicPPL.setindex!!(vi, 1.0, @varname(s));

julia> vi = DynamicPPL.setindex!!(vi, 2.0, @varname(m));

julia> values_as(vi, NamedTuple)
(s = 1.0, m = 2.0)

julia> values_as(vi, OrderedDict)
OrderedDict{Any, Any} with 2 entries:
  s => 1.0
  m => 2.0

julia> values_as(vi, Vector)
2-element Vector{Float64}:
 1.0
 2.0

abstract type AbstractTransformation

Represents a transformation to be used in link!! and invlink!!, amongst others.

A concrete implementation of this should implement the following methods:

• link!!: transforms the AbstractVarInfo to the unconstrained space.
• invlink!!: transforms the AbstractVarInfo to the constrained space.

And potentially:

• maybe_invlink_before_eval!!: hook to decide whether to transform before evaluating the model.

See also: link!!, invlink!!, maybe_invlink_before_eval!!.

struct NoTransformation <: DynamicPPL.AbstractTransformation

Transformation which applies the identity function.

struct DynamicTransformation <: DynamicPPL.AbstractTransformation

Transformation which transforms the variables on a per-need-basis in the execution of a given Model.

This is in constrast to StaticTransformation which transforms all variables before the execution of a given Model.

Different VarInfo types should implement their own methods for link!! and invlink!! for DynamicTransformation.

See also: StaticTransformation.

struct StaticTransformation{F} <: DynamicPPL.AbstractTransformation

Transformation which transforms all variables before the execution of a given Model.

This is done through the maybe_invlink_before_eval!! method.

See also: DynamicTransformation, maybe_invlink_before_eval!!.

Fields

• bijector::Any: The function, assumed to implement the Bijectors interface, to be applied to the variables

transformation(vi::AbstractVarInfo)

Return the AbstractTransformation related to vi.

LinkAll() <: AbstractLinkStrategy

Indicate that all variables should be linked.

UnlinkAll() <: AbstractLinkStrategy

Indicate that all variables should be unlinked.

LinkSome(vns) <: AbstractLinkStrategy

Indicate that the variables in vns must be linked. The link statuses of other variables are preserved. vns should be some iterable collection of VarNames, although there is no strict type requirement.

UnlinkSome(vns}) <: AbstractLinkStrategy

Indicate that the variables in vns must not Be linked. The link statuses of other variables are preserved. vns should be some iterable collection of VarNames, although there is no strict type requirement.

link([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
link([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their linked space without mutating vi.

Either transform all variables, or only ones specified in vns.

Use the transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, invlink.

invlink([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
invlink([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their constrained space without mutating vi.

Either transform all variables, or only ones specified in vns.

Use the (inverse of) transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, link.

link!!([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
link!!([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their linked space, mutating vi if possible.

Either transform all variables, or only ones specified in vns.

Use the transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, invlink!!.

invlink!!([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
invlink!!([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their constrained space, mutating vi if possible.

Either transform all variables, or only ones specified in vns.

Use the (inverse of) transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, link!!.

default_transformation(model::Model[, vi::AbstractVarInfo])

Return the AbstractTransformation currently related to model and, potentially, vi.

link_transform(dist)

Return the constrained-to-unconstrained bijector for distribution dist.

By default, this is just Bijectors.bijector(dist).

invlink_transform(dist)

Return the unconstrained-to-constrained bijector for distribution dist.

By default, this is just inverse(link_transform(dist)).

maybe_invlink_before_eval!!([t::Transformation,] vi, model)

Return a possibly invlinked version of vi.

This will be called prior to model evaluation, allowing one to perform a single invlink!! before evaluation rather than lazyily evaluating the transforms on as-we-need basis as is done with DynamicTransformation.

See also: StaticTransformation, DynamicTransformation.

Examples

julia> using DynamicPPL, Distributions, Bijectors

julia> @model demo() = x ~ Normal()
demo (generic function with 2 methods)

julia> # By subtyping `Transform`, we inherit the `(inv)link!!`.
       struct MyBijector <: Bijectors.Transform end

julia> # Define some dummy `inverse` which will be used in the `link!!` call.
       Bijectors.inverse(f::MyBijector) = identity

julia> # We need to define `with_logabsdet_jacobian` for `MyBijector`
       # (`identity` already has `with_logabsdet_jacobian` defined)
       function Bijectors.with_logabsdet_jacobian(::MyBijector, x)
           # Just using a large number of the logabsdet-jacobian term
           # for demonstration purposes.
           return (x, 1000)
       end

julia> # Change the `default_transformation` for our model to be a
       # `StaticTransformation` using `MyBijector`.
       function DynamicPPL.default_transformation(::Model{typeof(demo)})
           return DynamicPPL.StaticTransformation(MyBijector())
       end

julia> model = demo();

julia> vi = setindex!!(VarInfo(), 1.0, @varname(x));

julia> vi[@varname(x)]
1.0

julia> vi_linked = link!!(vi, model);

julia> # Now performs a single `invlink!!` before model evaluation.
       logjoint(model, vi_linked)
-1001.4189385332047

merge(varinfo, other_varinfos...)

Merge varinfos into one, giving precedence to the right-most varinfo when sensible.

This is particularly useful when combined with subset(varinfo, vns).

See docstring of subset(varinfo, vns) for examples.

subset(varinfo::AbstractVarInfo, vns::AbstractVector{<:VarName})

Subset a varinfo to only contain the variables vns.

The ordering of variables in the return value will be the same as in varinfo.

Examples

julia> @model function demo()
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           x = Vector{Float64}(undef, 2)
           x[1] ~ Normal(m, sqrt(s))
           x[2] ~ Normal(m, sqrt(s))
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> vi = VarInfo(model);

julia> keys(vi)
4-element Vector{VarName}:
 s
 m
 x[1]
 x[2]

julia> for (i, vn) in enumerate(keys(vi))
           vi = DynamicPPL.setindex!!(vi, Float64(i), vn)
       end

julia> vi[[@varname(s), @varname(m), @varname(x[1]), @varname(x[2])]]
4-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0

julia> # Extract one with only `m`.
       vi_subset1 = subset(vi, [@varname(m),]);

julia> keys(vi_subset1)
1-element Vector{VarName}:
 m

julia> vi_subset1[@varname(m)]
2.0

julia> # Extract one with both `s` and `x[2]`.
       vi_subset2 = subset(vi, [@varname(s), @varname(x[2])]);

julia> keys(vi_subset2)
2-element Vector{VarName}:
 s
 x[2]

julia> vi_subset2[[@varname(s), @varname(x[2])]]
2-element Vector{Float64}:
 1.0
 4.0

subset is particularly useful when combined with merge(vi::AbstractVarInfo)

julia> # Merge the two.
       vi_subset_merged = merge(vi_subset1, vi_subset2);

julia> keys(vi_subset_merged)
3-element Vector{VarName}:
 m
 s
 x[2]

julia> vi_subset_merged[[@varname(s), @varname(m), @varname(x[2])]]
3-element Vector{Float64}:
 1.0
 2.0
 4.0

julia> # Merge the two with the original.
       vi_merged = merge(vi, vi_subset_merged);

julia> keys(vi_merged)
4-element Vector{VarName}:
 s
 m
 x[1]
 x[2]

julia> vi_merged[[@varname(s), @varname(m), @varname(x[1]), @varname(x[2])]]
4-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0

Notes

Type-stability

unflatten!!(vi::AbstractVarInfo, x::AbstractVector)

Return a new instance of vi with the values of x assigned to the variables.

evaluate!!(model::Model, varinfo)

Evaluate the model with the given varinfo.

If the model has been marked as requiring threadsafe evaluation, are available, the varinfo provided will be wrapped in a ThreadSafeVarInfo before evaluation.

Returns a tuple of the model's return value, plus the updated varinfo (unwrapped if necessary).

struct DefaultContext <: AbstractContext end

DefaultContext, as the name suggests, is the default context used when instantiating a model.

julia> @model f() = x ~ Normal();

julia> model = f(); model.context
DefaultContext()

As an evaluation context, the behaviour of DefaultContext is to require all variables to be present in the AbstractVarInfo used for evaluation. Thus, semantically, evaluating a model with DefaultContext means 'calculating the log-probability associated with the variables in the AbstractVarInfo'.

InitContext(
        [rng::Random.AbstractRNG=Random.default_rng()],
        [strategy::AbstractInitStrategy=InitFromPrior()],
)

A leaf context that indicates that new values for random variables are currently being obtained through sampling. Used e.g. when initialising a fresh VarInfo. Note that, if leafcontext(model.context) isa InitContext, then evaluate!!(model, varinfo) will override all values in the VarInfo.

DynamicPPL.tilde_assume!!(
    context::AbstractContext,
    right::Distribution,
    vn::VarName,
    template::Any,
    vi::AbstractVarInfo
)::Tuple{Any,AbstractVarInfo}

Handle assumed variables, i.e. anything which is not observed (see tilde_observe!!). Accumulate the associated log probability, and return the sampled value and updated vi.

vn is the VarName on the left-hand side of the tilde statement.

template is the value of the top-level symbol in vn.

This function should return a tuple (x, vi), where x is the sampled value (which must be untransformed, i.e., insupport(right, x) must be true!) and vi is the updated VarInfo.

DynamicPPL.tilde_assume!!(
    ::DefaultContext,
    right::Distribution,
    vn::VarName,
    template::Any,
    vi::AbstractVarInfo
)

Handle assumed variables. For DefaultContext, this function extracts the value associated with vn from vi, If vi does not contain an appropriate value then this will error.

DynamicPPL.tilde_assume!!(
    context::AbstractContext,
    right::DynamicPPL.Submodel,
    vn::VarName,
    ::Any,
    vi::AbstractVarInfo
)

Evaluate the submodel with the given context.

DynamicPPL.tilde_observe!!(
    context::AbstractContext,
    right::Distribution,
    left,
    vn::Union{VarName, Nothing},
    vi::AbstractVarInfo
)::Tuple{Any,AbstractVarInfo}

This function handles observed variables, which may be:

• literals on the left-hand side, e.g., 3.0 ~ Normal()
• a model input, e.g. x ~ Normal() in a model @model f(x) ... end
• a conditioned or fixed variable, e.g. x ~ Normal() in a model model | (; x = 3.0).

The relevant log-probability associated with the observation is computed and accumulated in the VarInfo object vi (except for fixed variables, which do not contribute to the log-probability).

left is the actual value that the left-hand side evaluates to. vn is the VarName on the left-hand side, or nothing if the left-hand side is a literal value.

Observations of submodels are not yet supported in DynamicPPL.

This function should return a tuple (left, vi), where left is the same as the input, and vi is the updated VarInfo.

DynamicPPL.tilde_observe!!(
    ::DefaultContext,
    right::Distribution,
    left,
    vn::Union{VarName,Nothing},
    vi::AbstractVarInfo,
)

Handle observed variables. This just accumulates the log-likelihood for left.

AbstractParentContext

An abstract context that has a child context.

Subtypes of AbstractParentContext must implement the following interface:

• DynamicPPL.childcontext(context::AbstractParentContext): Return the child context.
• DynamicPPL.setchildcontext(parent::AbstractParentContext, child::AbstractContext): Reconstruct parent but now using child as its child context.

childcontext(context::AbstractParentContext)

Return the descendant context of context.

setchildcontext(parent::AbstractParentContext, child::AbstractContext)

Reconstruct parent but now using child is its childcontext, effectively updating the child context.

Examples

julia> using DynamicPPL: InitContext, ConditionContext

julia> ctx = ConditionContext((; a = 1));

julia> DynamicPPL.childcontext(ctx)
DefaultContext()

julia> ctx_prior = DynamicPPL.setchildcontext(ctx, InitContext(MersenneTwister(23), InitFromPrior()));

julia> DynamicPPL.childcontext(ctx_prior)
InitContext{MersenneTwister, InitFromPrior}(MersenneTwister(23), InitFromPrior())

leafcontext(context::AbstractContext)

Return the leaf of context, i.e. the first descendant context that is not an AbstractParentContext.

setleafcontext(left::AbstractContext, right::AbstractContext)

Return left but now with its leaf context replaced by right.

Note that this also works even if right is not a leaf context, in which case effectively append right to left, dropping the original leaf context of left.

Examples

julia> using DynamicPPL: leafcontext, setleafcontext, childcontext, setchildcontext, AbstractContext, InitContext

julia> struct ParentContext{C} <: AbstractParentContext
           context::C
       end

julia> DynamicPPL.childcontext(context::ParentContext) = context.context

julia> DynamicPPL.setchildcontext(::ParentContext, child) = ParentContext(child)

julia> Base.show(io::IO, c::ParentContext) = print(io, "ParentContext(", childcontext(c), ")")

julia> ctx = ParentContext(ParentContext(DefaultContext()))
ParentContext(ParentContext(DefaultContext()))

julia> # Replace the leaf context with another leaf.
       leafcontext(setleafcontext(ctx, InitContext(MersenneTwister(23), InitFromPrior())))
InitContext{MersenneTwister, InitFromPrior}(MersenneTwister(23), InitFromPrior())

julia> # Append another parent context.
       setleafcontext(ctx, ParentContext(DefaultContext()))
ParentContext(ParentContext(ParentContext(DefaultContext())))

setleafcontext(model::Model, context::AbstractContext)

Return a new Model with its leaf context set to context. This is a convenience shortcut for contextualize(model, setleafcontext(model.context, context)).

init!!(
    [rng::Random.AbstractRNG,]
    model::Model,
    varinfo::AbstractVarInfo,
    [init_strategy::AbstractInitStrategy=InitFromPrior()]
)

Evaluate the model and replace the values of the model's random variables in the given varinfo with new values, using a specified initialisation strategy. If the values in varinfo are not set, they will be added using a specified initialisation strategy.

If init_strategy is not provided, defaults to InitFromPrior().

Returns a tuple of the model's return value, plus the updated varinfo object.

InitFromPrior()

Obtain new values by sampling from the prior distribution.

InitFromUniform()
InitFromUniform(lower, upper)

Obtain new values by first transforming the distribution of the random variable to unconstrained space, then sampling a value uniformly between lower and upper, and transforming that value back to the original space.

If lower and upper are unspecified, they default to (-2, 2), which mimics Stan's default initialisation strategy.

Requires that lower <= upper.

References

Stan reference manual page on initialization

InitFromParams(
    params::Any
    fallback::Union{AbstractInitStrategy,Nothing}=InitFromPrior()
)

Obtain new values by extracting them from the given set of params.

The most common use case is to provide a NamedTuple or AbstractDict{<:VarName}, which provides a mapping from variable names to values. However, we leave the type of params open in order to allow for custom parameter storage types.

Custom parameter storage types

For InitFromParams to work correctly with a custom params::P, you need to implement

DynamicPPL.init(rng, vn::VarName, dist::Distribution, p::InitFromParams{P}) where {P}

This tells you how to obtain values for the random variable vn from p.params. Note that the last argument is InitFromParams(params), not just params itself. Please see the docstring of DynamicPPL.init for more information on the expected behaviour.

If you only use InitFromParams with DynamicPPL.OnlyAccsVarInfo, as is usually the case, then you will not need to implement anything else. So far, this is the same as you would do for creating any new AbstractInitStrategy subtype.

However, to use InitFromParams with a full DynamicPPL.VarInfo, you may also need to implement

DynamicPPL.get_param_eltype(p::InitFromParams{P}) where {P}

See the docstring of DynamicPPL.get_param_eltype for more information on when this is needed.

The argument fallback specifies how new values are to be obtained if they cannot be found in params, or they are specified as missing. fallback can either be an initialisation strategy itself, in which case it will be used to obtain new values, or it can be nothing, in which case an error will be thrown. The default for fallback is InitFromPrior().

AbstractInitStrategy

Abstract type representing the possible ways of initialising new values for the random variables in a model (e.g., when creating a new VarInfo).

Any subtype of AbstractInitStrategy must implement the DynamicPPL.init method, and in some cases, DynamicPPL.get_param_eltype (see its docstring for details).

init(rng::Random.AbstractRNG, vn::VarName, dist::Distribution, strategy::AbstractInitStrategy)

Generate a new value for a random variable with the given distribution.

This function must return an AbstractTransformedValue.

If strategy provides values that are already untransformed (e.g., a Float64 within (0, 1) for dist::Beta, then you should return an UntransformedValue.

Otherwise, often there are cases where this will return either a VectorValue or a LinkedVectorValue, for example, if the strategy is reading from an existing VarInfo.

DynamicPPL.get_param_eltype(strategy::AbstractInitStrategy)

Return the element type of the parameters generated from the given initialisation strategy.

The default implementation returns Any. However, for InitFromParams which provides known parameters for evaluating the model, methods are implemented in order to return more specific types.

In general, if you are implementing a custom AbstractInitStrategy, correct behaviour can only be guaranteed if you implement this method as well. However, quite often, the default return value of Any will actually suffice. The cases where this does not suffice, and where you do have to manually implement get_param_eltype, are explained in the extended help (see ??DynamicPPL.get_param_eltype in the REPL).

Extended help

There are a few edge cases in DynamicPPL where the element type is needed. These largely relate to determining the element type of accumulators ahead of time (before evaluation), as well as promoting type parameters in model arguments. The classic case is when evaluating a model with ForwardDiff: the accumulators must be set to Duals, and any Vector{Float64} arguments must be promoted to Vector{Dual}. Other tracer types, for example those in SparseConnectivityTracer.jl, also require similar treatment.

If the AbstractInitStrategy is never used in combination with tracer types, then it is perfectly safe to return Any. This does not lead to type instability downstream because the actual accumulators will still be created with concrete Float types (the Any is just used to determine whether the float type needs to be modified).

In case that wasn't enough: in fact, even the above is not always true. Firstly, the accumulator argument is only true when evaluating with ThreadSafeVarInfo. See the comments in DynamicPPL.unflatten!! for more details. For non-threadsafe evaluation, Julia is capable of automatically promoting the types on its own. Secondly, the promotion only matters if you are trying to directly assign into a Vector{Float64} with a ForwardDiff.Dual or similar tracer type, for example using xs[i] = MyDual. This doesn't actually apply to tilde-statements like xs[i] ~ ... because those use Accessors.set under the hood, which also does the promotion for you. For the gory details, see the following issues:

• https://github.com/TuringLang/DynamicPPL.jl/issues/906 for accumulator types
• https://github.com/TuringLang/DynamicPPL.jl/issues/823 for type argument promotion

get_param_eltype(varinfo::AbstractVarInfo, context::AbstractContext)

Get the element type of the parameters being used to evaluate a model, using a varinfo under the given context. For example, when evaluating a model with ForwardDiff AD, this should return ForwardDiff.Dual.

By default, this uses eltype(varinfo) which is slightly cursed. This relies on the fact that typically, before evaluation, the parameters will have been inserted into the VarInfo's metadata field.

For InitContext, it's quite different: because InitContext is responsible for supplying the parameters, we can avoid using eltype(varinfo) and instead query the parameters inside it. See the docstring of get_param_eltype(strategy::AbstractInitStrategy) for more explanation.

AbstractTransformedValue

An abstract type for values that enter the DynamicPPL tilde-pipeline.

These values are generated by an AbstractInitStrategy: the function DynamicPPL.init should return an AbstractTransformedValue.

Each AbstractTransformedValue contains some version of the actual variable's value, together with a transformation that can be used to convert the internal value back to the original space.

Current subtypes are VectorValue, LinkedVectorValue, and UntransformedValue. DynamicPPL's VarInfo type stores either VectorValues or LinkedVectorValues internally, depending on the link status of the VarInfo.

Subtypes of this should implement the following functions:

• DynamicPPL.get_transform(tv::AbstractTransformedValue): Get the transformation that converts the internal value back to the original space.

• DynamicPPL.get_internal_value(tv::AbstractTransformedValue): Get the internal value stored in tv.

• DynamicPPL.set_internal_value(tv::AbstractTransformedValue, new_val): Create a new AbstractTransformedValue with the same transformation as tv, but with internal value new_val.

• DynamicPPL.VarNamedTuples.vnt_size(tv::AbstractTransformedValue): Get the size of the original value before transformation.

VectorValue{V<:AbstractVector,T,S}

A transformed value that stores its internal value as a vectorised form. This is what VarInfo sees as an "unlinked value".

These values can be generated when using InitFromParams with a VarInfo's internal values.

LinkedVectorValue{V<:AbstractVector,T,S}

A transformed value that stores its internal value as a linked andvectorised form. This is what VarInfo sees as a "linked value".

These values can be generated when using InitFromParams with a VarInfo's internal values.

UntransformedValue{V}

A raw, untransformed, value.

These values can be generated from initialisation strategies such as InitFromPrior, InitFromUniform, and InitFromParams on a standard container type.

get_transform(tv::AbstractTransformedValue)

Get the transformation that converts the internal value back to the raw value.

get_internal_value(tv::AbstractTransformedValue)

Get the internal value stored in tv.

set_internal_value(tv::AbstractTransformedValue, new_val)

Create a new AbstractTransformedValue with the same transformation as tv, but with internal value new_val.

ParamsWithStats

A struct which contains parameter values extracted from a VarInfo, along with any statistics associated with the VarInfo. The statistics are provided as a NamedTuple and are optional.

AbstractMCMC.from_samples(
    ::Type{MCMCChains.Chains},
    params_and_stats::AbstractMatrix{<:DynamicPPL.ParamsWithStats}
)

Convert an array of DynamicPPL.ParamsWithStats to an MCMCChains.Chains object.

AbstractMCMC.to_samples(
    ::Type{DynamicPPL.ParamsWithStats},
    chain::MCMCChains.Chains,
)

Convert an MCMCChains.Chains object to an array of DynamicPPL.ParamsWithStats.

For this to work, chain must contain the varname_to_symbol mapping in its info field.