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

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, varinfo::AbstractVarInfo)

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

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.

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

Return the log prior probability of variables θ 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, varinfo::AbstractVarInfo)

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

See also logjoint and logprior.

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

Return the log likelihood of variables θ 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, varinfo::AbstractVarInfo)

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

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 logprior and loglikelihood.

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

Return the log joint probability of variables θ 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

LogDensityFunction(
    model::Model,
    getlogdensity::Function=getlogjoint_internal,
    varinfo::AbstractVarInfo=ldf_default_varinfo(model, 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.

There are several options for getlogdensity 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.

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

• model: model used for evaluation

• getlogdensity: function to be called on varinfo to extract the log density. By default getlogjoint_internal.

• varinfo: varinfo used for evaluation. If not specified, generated with ldf_default_varinfo.

• adtype: AD type used for evaluation of log density gradient. If nothing, no gradient can be calculated

• prep: (internal use only) gradient preparation object for the model

Examples

julia> using Distributions

julia> using DynamicPPL: LogDensityFunction, setaccs!!

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

julia> model = demo(1.0);

julia> f = LogDensityFunction(model);

julia> # It implements the interface of LogDensityProblems.jl.
       using LogDensityProblems

julia> LogDensityProblems.logdensity(f, [0.0])
-2.3378770664093453

julia> LogDensityProblems.dimension(f)
1

julia> # By default it uses `VarInfo` under the hood, but this is not necessary.
       f = LogDensityFunction(model, getlogjoint_internal, SimpleVarInfo(model));

julia> LogDensityProblems.logdensity(f, [0.0])
-2.3378770664093453

julia> # One can also specify evaluating e.g. the log prior only:
       f_prior = LogDensityFunction(model, getlogprior);

julia> LogDensityProblems.logdensity(f_prior, [0.0]) == logpdf(Normal(), 0.0)
true

julia> # If we also need to calculate the gradient, we can specify an AD backend.
       import ForwardDiff, ADTypes

julia> f = LogDensityFunction(model, adtype=ADTypes.AutoForwardDiff());

julia> LogDensityProblems.logdensity_and_gradient(f, [0.0])
(-2.3378770664093453, [1.0])

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

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, Accessors.PropertyLens{:x}}}:
 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, Accessors.PropertyLens{:m}}, Float64} with 1 entry:
  a.m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(cm))
1-element Vector{VarName{:a, Accessors.PropertyLens{:x}}}:
 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

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.
       cm = fix(contextualize(m, PrefixContext(@varname(a), fix(m=1.0))), x=100.0);

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

julia> keys(VarInfo(cm))
1-element Vector{VarName{:a, Accessors.PropertyLens{:x}}}:
 a.x

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

julia> fixed(cm)
Dict{VarName{:a, Accessors.PropertyLens{:m}}, Float64} with 1 entry:
  a.m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(cm))
1-element Vector{VarName{:a, Accessors.PropertyLens{:x}}}:
 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,)

@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::NamedTuple)
returned(model::Model, parameters::AbstractDict{<:VarName})

Execute model with variables keys set to values and return the values returned by the model.

returned(model::Model, values, keys)

Execute model with variables keys set to values and return the values returned by the model. This method is deprecated; use the NamedTuple or AbstractDict version instead.

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.

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 type-stable `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 `θ`.
       # Flip `x` so we hit the other support of `y`.
       θ = [!varinfo[@varname(x)], rand(rng)];

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

julia> # Determine the expected support of `y`.
       lb, ub = θ[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> # (×) Fails! Because `VarInfo` _saves_ the original distributions
       # used in the very first model evaluation, hence the support of `y`
       # is not updated even though `x` has changed.
       lb ≤ first(varinfo_invlinked[@varname(y)]) ≤ ub
false

julia> # Approach 2: Extract realizations using `values_as_in_model`.
       # (✓) `values_as_in_model` will re-run the model and extract
       # the correct realization of `y` given the new values of `x`.
       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 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

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

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

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

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 a collection of VarName as they are expected to appear in the model.

Even though it is recommended to implement this by hand for a particular Model, a default implementation using SimpleVarInfo{<:Dict} is provided.

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 get, e.g. get(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, SimpleVarInfo.

struct VarInfo{Tmeta,Accs<:AccumulatorTuple} <: AbstractVarInfo
    metadata::Tmeta
    accs::Accs
end

A light wrapper over some kind of metadata.

The type of the metadata can be one of a number of options. It may either be a Metadata or a VarNamedVector, or, it may be a NamedTuple which maps symbols to Metadata or VarNamedVector instances. Here, a symbol refers to a Julia variable and may consist of one or more VarNames which appear on the left-hand side of tilde statements. For example, x[1] and x[2] both have the same symbol x.

Several type aliases are provided for these forms of VarInfos:

• VarInfo{<:Metadata} is UntypedVarInfo
• VarInfo{<:VarNamedVector} is UntypedVectorVarInfo
• VarInfo{<:NamedTuple} is NTVarInfo

The NamedTuple form, i.e. NTVarInfo, is useful for maintaining type stability of model evaluation. However, the element type of NamedTuples are not contained in its type itself: thus, there is no way to use the type system to determine whether the elements of the NamedTuple are Metadata or VarNamedVector.

Note that for NTVarInfo, it is the user's responsibility to ensure that each symbol is visited at least once during model evaluation, regardless of any stochastic branching.

untyped_varinfo([rng, ]model[, init_strategy])

Construct a VarInfo object for the given model, which has just a single Metadata as its metadata field.

Arguments

• rng::Random.AbstractRNG: The random number generator to use during model evaluation
• model::Model: The model for which to create the varinfo object
• init_strategy::AbstractInitStrategy: How the values are to be initialised. Defaults to InitFromPrior().

typed_varinfo(vi::UntypedVarInfo)

This function finds all the unique syms from the instances of VarName{sym} found in vi.metadata.vns. It then extracts the metadata associated with each symbol from the global vi.metadata field. Finally, a new VarInfo is created with a new metadata as a NamedTuple mapping from symbols to type-stable Metadata instances, one for each symbol.

typed_varinfo([rng, ]model[, init_strategy])

Return a VarInfo object for the given model, which has a NamedTuple of Metadata structs as its metadata field.

Arguments

• rng::Random.AbstractRNG: The random number generator to use during model evaluation
• model::Model: The model for which to create the varinfo object
• init_strategy::AbstractInitStrategy: How the values are to be initialised. Defaults to InitFromPrior().

untyped_vector_varinfo([rng, ]model[, init_strategy])

Return a VarInfo object for the given model, which has just a single VarNamedVector as its metadata field.

Arguments

• rng::Random.AbstractRNG: The random number generator to use during model evaluation
• model::Model: The model for which to create the varinfo object
• init_strategy::AbstractInitStrategy: How the values are to be initialised. Defaults to InitFromPrior().

typed_vector_varinfo([rng, ]model[, init_strategy])

Return a VarInfo object for the given model, which has a NamedTuple of VarNamedVectors as its metadata field.

Arguments

• rng::Random.AbstractRNG: The random number generator to use during model evaluation
• model::Model: The model for which to create the varinfo object
• init_strategy::AbstractInitStrategy: How the values are to be initialised. Defaults to InitFromPrior().

is_transformed(vnv::VarNamedVector, vn::VarName)

Return a boolean for whether vn is guaranteed to have been transformed so that its domain is all of Euclidean space.

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.

is_transformed(vi::VarInfo)

Check whether vi is in the transformed space.

Turing's Hamiltonian samplers use the link and invlink functions from Bijectors.jl to map a constrained variable (for example, one bounded to the space [0, 1]) from its constrained space to the set of real numbers. is_transformed checks if the number is in the constrained space or the real space.

If some but only some of the variables in vi are transformed, this function will return true. This behavior will likely change in the future.

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.

empty!(meta::Metadata)

Empty the fields of meta.

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

struct SimpleVarInfo{NT, Accs<:DynamicPPL.AccumulatorTuple, C<:DynamicPPL.AbstractTransformation} <: AbstractVarInfo

A simple wrapper of the parameters with a logp field for accumulation of the logdensity.

Currently only implemented for NT<:NamedTuple and NT<:AbstractDict.

Fields

• values: underlying representation of the realization represented

• accs: tuple of accumulators for things like log prior and log likelihood

• transformation: represents whether it assumes variables to be transformed

Notes

The major differences between this and NTVarInfo are:

1. SimpleVarInfo does not require linearization.
2. SimpleVarInfo can use more efficient bijectors.
3. SimpleVarInfo is only type-stable if NT<:NamedTuple and either a) no indexing is used in tilde-statements, or b) the values have been specified with the correct shapes.

Examples

General usage

julia> using StableRNGs

julia> @model function demo()
           m ~ Normal()
           x = Vector{Float64}(undef, 2)
           for i in eachindex(x)
               x[i] ~ Normal()
           end
           return x
       end
demo (generic function with 2 methods)

julia> m = demo();

julia> rng = StableRNG(42);

julia> # In the `NamedTuple` version we need to provide the place-holder values for
       # the variables which are using "containers", e.g. `Array`.
       # In this case, this means that we need to specify `x` but not `m`.
       _, vi = DynamicPPL.init!!(rng, m, SimpleVarInfo((x = ones(2), )));

julia> # (✓) Vroom, vroom! FAST!!!
       vi[@varname(x[1])]
0.4471218424633827

julia> # We can also access arbitrary varnames pointing to `x`, e.g.
       vi[@varname(x)]
2-element Vector{Float64}:
 0.4471218424633827
 1.3736306979834252

julia> vi[@varname(x[1:2])]
2-element Vector{Float64}:
 0.4471218424633827
 1.3736306979834252

julia> # (×) If we don't provide the container...
       _, vi = DynamicPPL.init!!(rng, m, SimpleVarInfo());
ERROR: FieldError: type NamedTuple has no field `x`, available fields: `m`
[...]

julia> # If one does not know the varnames, we can use a `OrderedDict` instead.
       _, vi = DynamicPPL.init!!(rng, m, SimpleVarInfo{Float64}(OrderedDict{VarName,Any}()));

julia> # (✓) Sort of fast, but only possible at runtime.
       vi[@varname(x[1])]
-1.019202452456547

julia> # In addtion, we can only access varnames as they appear in the model!
       vi[@varname(x)]
ERROR: x was not found in the dictionary provided
[...]

julia> vi[@varname(x[1:2])]
ERROR: x[1:2] was not found in the dictionary provided
[...]

Technically, it's possible to use any implementation of AbstractDict in place of OrderedDict, but OrderedDict ensures that certain operations, e.g. linearization/flattening of the values in the varinfo, are consistent between evaluations. Hence OrderedDict is the preferred implementation of AbstractDict to use here.

You can also sample in transformed space:

julia> @model demo_constrained() = x ~ Exponential()
demo_constrained (generic function with 2 methods)

julia> m = demo_constrained();

julia> _, vi = DynamicPPL.init!!(rng, m, SimpleVarInfo());

julia> vi[@varname(x)] # (✓) 0 ≤ x < ∞
1.8632965762164932

julia> _, vi = DynamicPPL.init!!(rng, m, DynamicPPL.set_transformed!!(SimpleVarInfo(), true));

julia> vi[@varname(x)] # (✓) -∞ < x < ∞
-0.21080155351918753

julia> xs = [last(DynamicPPL.init!!(rng, m, DynamicPPL.set_transformed!!(SimpleVarInfo(), true)))[@varname(x)] for i = 1:10];

julia> any(xs .< 0)  # (✓) Positive probability mass on negative numbers!
true

julia> # And with `OrderedDict` of course!
       _, vi = DynamicPPL.init!!(rng, m, DynamicPPL.set_transformed!!(SimpleVarInfo(OrderedDict{VarName,Any}()), true));

julia> vi[@varname(x)] # (✓) -∞ < x < ∞
0.6225185067787314

julia> xs = [last(DynamicPPL.init!!(rng, m, DynamicPPL.set_transformed!!(SimpleVarInfo(), true)))[@varname(x)] for i = 1:10];

julia> any(xs .< 0) # (✓) Positive probability mass on negative numbers!
true

Evaluation in transformed space of course also works:

julia> vi = DynamicPPL.set_transformed!!(SimpleVarInfo((x = -1.0,)), true)
Transformed SimpleVarInfo((x = -1.0,), (LogPrior = LogPriorAccumulator(0.0), LogJacobian = LogJacobianAccumulator(0.0), LogLikelihood = LogLikelihoodAccumulator(0.0)))

julia> # (✓) Positive probability mass on negative numbers!
       getlogjoint_internal(last(DynamicPPL.evaluate!!(m, vi)))
-1.3678794411714423

julia> # While if we forget to indicate that it's transformed:
       vi = DynamicPPL.set_transformed!!(SimpleVarInfo((x = -1.0,)), false)
SimpleVarInfo((x = -1.0,), (LogPrior = LogPriorAccumulator(0.0), LogJacobian = LogJacobianAccumulator(0.0), LogLikelihood = LogLikelihoodAccumulator(0.0)))

julia> # (✓) No probability mass on negative numbers!
       getlogjoint_internal(last(DynamicPPL.evaluate!!(m, vi)))
-Inf

Indexing

Using NamedTuple as underlying storage.

julia> svi_nt = SimpleVarInfo((m = (a = [1.0], ), ));

julia> svi_nt[@varname(m)]
(a = [1.0],)

julia> svi_nt[@varname(m.a)]
1-element Vector{Float64}:
 1.0

julia> svi_nt[@varname(m.a[1])]
1.0

julia> svi_nt[@varname(m.a[2])]
ERROR: BoundsError: attempt to access 1-element Vector{Float64} at index [2]
[...]

julia> svi_nt[@varname(m.b)]
ERROR: FieldError: type NamedTuple has no field `b`, available fields: `a`
[...]

Using OrderedDict as underlying storage.

julia> svi_dict = SimpleVarInfo(OrderedDict(@varname(m) => (a = [1.0], )));

julia> svi_dict[@varname(m)]
(a = [1.0],)

julia> svi_dict[@varname(m.a)]
1-element Vector{Float64}:
 1.0

julia> svi_dict[@varname(m.a[1])]
1.0

julia> svi_dict[@varname(m.a[2])]
ERROR: m.a[2] was not found in the dictionary provided
[...]

julia> svi_dict[@varname(m.b)]
ERROR: m.b was not found in the dictionary provided
[...]

DynamicPPL.tilde_assume!!(
    context::AbstractContext,
    right::Distribution,
    vn::VarName,
    vi::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.

This function should return a tuple (x, vi), where x is the sampled value (which must be in unlinked space!) and vi is the updated VarInfo.

DynamicPPL.tilde_assume!!(
    ::DefaultContext, right::Distribution, vn::VarName, 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,
    vi::AbstractVarInfo
)

Evaluate the submodel with the given context.

DynamicPPL.tilde_observe!!(
    context::AbstractContext,
    right::Distribution,
    left,
    vn::Union{VarName, Nothing},
    vi::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.

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, logjac, vn, dist)
• 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.
• 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.

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.

push!!(vi::VarInfo, vn::VarName, r, dist::Distribution)

Push a new random variable vn with a sampled value r from a distribution dist to the VarInfo vi, mutating if it makes sense.

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!(vnv::VarNamedVector, val, i::Int)

Sets the ith element of the internal storage vector, ignoring inactive entries.

setindex_internal!(vnv::VarNamedVector, val, vn::VarName[, transform])

Like setindex!, but sets the values as they are stored internally in vnv.

Optionally can set the transformation, such that transform(val) is the original value of the variable. By default, the transform is the identity if creating a new entry in vnv, or the existing transform if updating an existing entry.

update_internal!(vnv::VarNamedVector, vn::VarName, val::AbstractVector[, transform])

Update an existing entry for vn in vnv with the value val.

Like setindex_internal!, but errors if the key vn doesn't exist.

transform should be a function that converts val to the original representation. By default it's the same as the old transform for vn.

insert_internal!(vnv::VarNamedVector, val::AbstractVector, vn::VarName[, transform])

Add a variable with given value to vnv.

Like setindex_internal!, but errors if the key vn already exists.

transform should be a function that converts val to the original representation. By default it's identity.

length_internal(vnv::VarNamedVector)

Return the length of the internal storage vector of vnv, ignoring inactive entries.

reset!(vnv::VarNamedVector, val, vn::VarName)

Reset the value of vn in vnv to val.

This differs from setindex! in that it will always change the transform of the variable to be the default vectorisation transform. This undoes any possible linking.

Examples

julia> using DynamicPPL: VarNamedVector, @varname, reset!

julia> vnv = VarNamedVector{VarName,Any,Any}();

julia> vnv[@varname(x)] = reshape(1:9, (3, 3));

julia> setindex!(vnv, 2.0, @varname(x))
ERROR: An error occurred while assigning the value 2.0 to variable x. If you are changing the type or size of a variable you'll need to call reset!
[...]

julia> reset!(vnv, 2.0, @varname(x));

julia> vnv[@varname(x)]
2.0

update!(vnv::VarNamedVector, val, vn::VarName)

Update the value of vn in vnv to val.

Like setindex!, but errors if the key vn doesn't exist.

insert!(vnv::VarNamedVector, val, vn::VarName)

Add a variable with given value to vnv.

Like setindex!, but errors if the key vn already exists.

loosen_types!!(vnv::VarNamedVector, ::Type{KNew}, ::Type{VNew}, ::Type{TNew})

Loosen the types of vnv to allow varname type KNew and transformation type TransNew.

If KNew is a subtype of K and TransNew is a subtype of the element type of the TTrans then this is a no-op and vnv is returned as is. Otherwise a new VarNamedVector is returned with the same data but more abstract types, so that variables of type KNew and transformations of type TransNew can be pushed to it. Some of the underlying storage is shared between vnv and the return value, and thus mutating one may affect the other.

See also

tighten_types!!

Examples

julia> using DynamicPPL: VarNamedVector, @varname, loosen_types!!, setindex_internal!

julia> vnv = VarNamedVector(@varname(x) => [1.0]);

julia> y_trans(x) = reshape(x, (2, 2));

julia> setindex_internal!(vnv, collect(1:4), @varname(y), y_trans)
ERROR: MethodError: Cannot `convert` an object of type
[...]

julia> vnv_loose = DynamicPPL.loosen_types!!(
           vnv, typeof(@varname(y)), Float64, typeof(y_trans)
       );

julia> setindex_internal!(vnv_loose, collect(1:4), @varname(y), y_trans)

julia> vnv_loose[@varname(y)]
2×2 Matrix{Float64}:
 1.0  3.0
 2.0  4.0

tighten_types!!(vnv::VarNamedVector)

Return a VarNamedVector like vnv with the most concrete types possible.

This function either returns vnv itself or new VarNamedVector with the same values in it, but with the element types of various containers made as concrete as possible.

For instance, if vnv has its vector of transforms have eltype Any, but all the transforms are actually identity transformations, this function will return a new VarNamedVector with the transforms vector having eltype typeof(identity).

This is a lot like the reverse of loosen_types!!. Like with loosen_types!!, the return value may share some of its underlying storage with vnv, and thus mutating one may affect the other.

See also

loosen_types!!

Examples

julia> using DynamicPPL: VarNamedVector, @varname, loosen_types!!, setindex_internal!

julia> vnv = VarNamedVector(@varname(x) => Real[23], @varname(y) => randn(2,2));

julia> vnv = delete!(vnv, @varname(y));

julia> eltype(vnv)
Real

julia> vnv.transforms
1-element Vector{Any}:
 identity (generic function with 1 method)

julia> vnv_tight = DynamicPPL.tighten_types!!(vnv);

julia> eltype(vnv_tight) == Int
true

julia> vnv_tight.transforms
1-element Vector{typeof(identity)}:
 identity (generic function with 1 method)

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

SimpleVarInfo with NamedTuple:

julia> data = (x = 1.0, m = [2.0]);

julia> values_as(SimpleVarInfo(data))
(x = 1.0, m = [2.0])

julia> values_as(SimpleVarInfo(data), NamedTuple)
(x = 1.0, m = [2.0])

julia> values_as(SimpleVarInfo(data), OrderedDict)
OrderedDict{VarName{sym, typeof(identity)} where sym, Any} with 2 entries:
  x => 1.0
  m => [2.0]

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

SimpleVarInfo with OrderedDict:

julia> data = OrderedDict{Any,Any}(@varname(x) => 1.0, @varname(m) => [2.0]);

julia> values_as(SimpleVarInfo(data))
OrderedDict{Any, Any} with 2 entries:
  x => 1.0
  m => [2.0]

julia> values_as(SimpleVarInfo(data), NamedTuple)
(x = 1.0, m = [2.0])

julia> values_as(SimpleVarInfo(data), OrderedDict)
OrderedDict{Any, Any} with 2 entries:
  x => 1.0
  m => [2.0]

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

VarInfo with NamedTuple of Metadata:

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

julia> vi[@varname(s)] = 1.0; vi[@varname(m)] = 2.0;

julia> # For the sake of brevity, let's just check the type.
       md = values_as(vi); md.s isa Union{DynamicPPL.Metadata, DynamicPPL.VarNamedVector}
true

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

julia> values_as(vi, OrderedDict)
OrderedDict{VarName{sym, typeof(identity)} where sym, Float64} with 2 entries:
  s => 1.0
  m => 2.0

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

VarInfo with Metadata:

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

julia> vi[@varname(s)] = 1.0; vi[@varname(m)] = 2.0;

julia> # For the sake of brevity, let's just check the type.
       values_as(vi) isa Union{DynamicPPL.Metadata, Vector}
true

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

julia> values_as(vi, OrderedDict)
OrderedDict{VarName{sym, typeof(identity)} where sym, Float64} with 2 entries:
  s => 1.0
  m => 2.0

julia> values_as(vi, Vector)
2-element Vector{Real}:
 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.

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.

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 = SimpleVarInfo(x=1.0)
SimpleVarInfo((x = 1.0,), 0.0)

julia> # Uses the `inverse` of `MyBijector`, which we have defined as `identity`
       vi_linked = link!!(vi, model)
Transformed SimpleVarInfo((x = 1.0,), 0.0)

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> varinfo = VarInfo(model);

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

julia> for (i, vn) in enumerate(keys(varinfo))
           varinfo[vn] = i
       end

julia> varinfo[[@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`.
       varinfo_subset1 = subset(varinfo, [@varname(m),]);


julia> keys(varinfo_subset1)
1-element Vector{VarName{:m, typeof(identity)}}:
 m

julia> varinfo_subset1[@varname(m)]
2.0

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

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

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

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

julia> # Merge the two.
       varinfo_subset_merged = merge(varinfo_subset1, varinfo_subset2);

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

julia> varinfo_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.
       varinfo_merged = merge(varinfo, varinfo_subset_merged);

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

julia> varinfo_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 multiple threads 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'.

PrefixContext(vn::VarName[, context::AbstractContext])
PrefixContext(vn::Val{sym}[, context::AbstractContext]) where {sym}

Create a context that allows you to use the wrapped context when running the model and prefixes all parameters with the VarName vn.

PrefixContext(Val(:a), context) is equivalent to PrefixContext(@varname(a), context). If context is not provided, it defaults to DefaultContext().

This context is useful in nested models to ensure that the names of the parameters are unique.

See also: to_submodel

ConditionContext{Values<:Union{NamedTuple,AbstractDict},Ctx<:AbstractContext}

Model context that contains values that are to be conditioned on. The values can either be a NamedTuple mapping symbols to values, such as (a=1, b=2), or an AbstractDict mapping varnames to values (e.g. Dict(@varname(a) => 1, @varname(b) => 2)). The former is more performant, but the latter must be used when there are varnames that cannot be represented as symbols, e.g. @varname(x[1]).

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.

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::Union{AbstractDict{<:VarName},NamedTuple},
    fallback::Union{AbstractInitStrategy,Nothing}=InitFromPrior()
)

Obtain new values by extracting them from the given dictionary or NamedTuple.

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

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

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

determine_suitable_varinfo(model; only_dppl::Bool=true)

Return a suitable varinfo for the given model.

See also: DynamicPPL.Experimental.is_suitable_varinfo.

Arguments

• model: The model for which to determine the varinfo.

Keyword Arguments

• only_dppl: If true, only consider error reports within DynamicPPL.jl.

Examples

julia> using DynamicPPL.Experimental: determine_suitable_varinfo

julia> using JET: JET  # needs to be loaded for full functionality

julia> @model function model_with_random_support()
           x ~ Bernoulli()
           if x
               y ~ Normal()
           else
               z ~ Normal()
           end
       end
model_with_random_support (generic function with 2 methods)

julia> model = model_with_random_support();

julia> # Typed varinfo cannot handle this random support model properly
       # as using a single execution of the model will not see all random variables.
       # Hence, this this model requires untyped varinfo.
       vi = determine_suitable_varinfo(model);
┌ Warning: Model seems incompatible with typed varinfo. Falling back to untyped varinfo.
└ @ DynamicPPLJETExt ~/.julia/dev/DynamicPPL.jl/ext/DynamicPPLJETExt.jl:48

julia> vi isa typeof(DynamicPPL.untyped_varinfo(model))
true

julia> # In contrast, a simple model with no random support can be handled by typed varinfo.
       @model model_with_static_support() = x ~ Normal()
model_with_static_support (generic function with 2 methods)

julia> vi = determine_suitable_varinfo(model_with_static_support());

julia> vi isa typeof(DynamicPPL.typed_varinfo(model_with_static_support()))
true

is_suitable_varinfo(model::Model, varinfo::AbstractVarInfo; kwargs...)

Check if the model supports evaluation using the provided varinfo.

Arguments

• model: The model to verify the support for.
• varinfo: The varinfo to verify the support for.

Keyword Arguments

• only_dppl: If true, only consider error reports occuring in the tilde pipeline. Default: true.

Returns

• issuccess: true if the model supports the varinfo, otherwise false.
• report: The result of report_call from JET.jl.