It’s more user-friendly and it’s much easier to work with programmatically.
For the former, it’s just nice to be able to tell the model to no longer consider a variable to be an observation (you can do this using decondition), e.g. if you want to look at the posterior predictive. The way of doing this using the arguments is to pass missing in its place, but this sooo often goes wrong and we’ve experienced users being generally quite confused about this.
When observations are passed as an argument, it’s not possible for me to query the model itself as to what is considered an observation without execution (we have to execute the model to determine which of the arguments end up on the left-hand side of a ~). In contrast, if you use the model | ... or condition(model, ...) syntax, the observations are effectively attached to the model in a way that we can then easily inspect using observations.
For example, we are working on simplifying the inner workings of our Gibbs sampler, where we now can simply use the condition(model, ...) syntax to sequentially fix different variables.
With all that being said, in the example with the Kalman filter you give above, it’s very understandable that you’d want to pass it as an argument. But in this case your ys never appear on the LHS of ~, and Turing won’t consider this argument as an observation.
So this has nothing to do with condition vs. pass-obs-as-arg, but it’s just the restrictions you get when using @addlogprob!. This, AFAIK, is the same in every other PPL, i.e. if you use the “manual accumulation of log-probability”, then you lose some of the convenience.
BUT, we can “fix” this! If you want to make the Turing consider your ys as an observations, for things like loo, we can implement a “distribution” that effectively does what @addlogprob! does but at the same time, since it’s a Distribution, signals to Turing.jl that “hey, this thing on the LHS of this ~ is an observation!”.
"""
DiracWithProb(value[, logp])
A Dirac delta distribution with a single point mass at `value` and `logpdf` equal to `logp`.
This is basically the same as `Distributions.Dirac` but without
being restricted to univariate values, in addition to having a non-trivial `logpdf`.
!!! warning
The non-trivial `logpdf` makes it an improper distribution, but it can be useful
in certain scenarios when used in a model.
"""
struct DiracWithProb{T1,T2,V} <: Distributions.DiscreteDistribution{V}
value::T1
logp::T2
end
variateform(x::Real) = Distributions.Univariate
variateform(x::AbstractArray{<:Real,N}) where {N} = Distributions.ArrayLikeVariate{N}
DiracWithProb(x) = DiracWithProb(x, zero(eltype(x)))
DiracWithProb(x, logp) = DiracWithProb{typeof(x),typeof(logp),variateform(x)}(x, logp)
const UnivariateDiracWithProb{T1,T2} = DiracWithProb{T1,T2,Distributions.Univariate}
const ArrayVariateDiracWithProb{T1,T2,N} = DiracWithProb{T1,T2,Distributions.ArrayLikeVariate{N}}
Base.eltype(::Type{DiracWithProb{T}}) where {T} = T
# Some required stuff from Distributions.jl.
Distributions.insupport(d::UnivariateDiracWithProb, x::Real) = x == d.value
Distributions.insupport(d::DiracWithProb, x::AbstractArray) = x == d.value
Base.minimum(d::DiracWithProb) = d.value
Base.maximum(d::DiracWithProb) = d.value
Distributions.support(d::DiracWithProb) = (d.value,)
# Evaluation
Distributions.pdf(d::UnivariateDiracWithProb, x::Real) = insupport(d, x) ? exp(d.logp) : zero(d.logp)
Distributions.logpdf(d::UnivariateDiracWithProb, x::Real) = insupport(d, x) ? d.logp : typeof(d.logp)(-Inf)
function Distributions.pdf(d::ArrayVariateDiracWithProb{<:Any,N}, x::AbstractArray{<:Real,N}) where {N}
return insupport(d, x) ? exp(d.logp) : zero(d.logp)
end
function Distributions.logpdf(d::ArrayVariateDiracWithProb{<:Any,N}, x::AbstractArray{<:Real,N}) where {N}
return insupport(d, x) ? d.logp : typeof(d.logp)(-Inf)
end
# Sampling
Base.rand(rng::Random.AbstractRNG, d::UnivariateDiracWithProb) = d.value
Base.rand(rng::Random.AbstractRNG, d::ArrayVariateDiracWithProb) = d.value
Equipped with this, we can do stuff like
julia> @model function demo(x, y)
s ~ InverseGamma(2, 3)
m ~ Normal(0, √s)
x ~ DiracWithProb(x, logpdf(Normal(m, √s), x))
y ~ DiracWithProb(y, logpdf(Normal(m, √s), y))
end
demo (generic function with 4 methods)
julia> model = demo(1.5, 2);
julia> chain = sample(model, NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 0.8
Sampling 100%|████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00
Chains MCMC chain (1000×14×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.08 seconds
Compute duration = 0.08 seconds
parameters = s, m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
s 1.9040 1.4629 0.0561 721.9248 664.9316 0.9992 9499.0100
m 1.1291 0.7297 0.0314 556.9236 473.0818 1.0012 7327.9420
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s 0.5666 1.0414 1.4280 2.2449 6.0441
m -0.4221 0.7014 1.1622 1.5815 2.5344
julia> pointwise_loglikelihoods(model, MCMCChains.get_sections(chain, :parameters))
Dict{String, Matrix{Float64}} with 2 entries:
"x" => [-1.11969; -0.859207; … ; -1.1811; -1.15306;;]
"y" => [-1.20779; -1.31884; … ; -0.877885; -0.944768;;]
which will allow it to work with loo, etc.