Experimental Gibbs seems to have higher variance than "stable" Gibbs (maybe)

[torfjelde]

Turing.Experimental.Gibbs often causes test-suites to fail due to high-variance estimates for models in DynamicPPL.TestUtils.DEMO_MODELS.

At the moment, it's somewhat unclear to me whether these are due to more extensive testing of Turing.Experimental.Gibbs (unceratin if Gibbs is actually being tested as rigourously on DEMO_MODELS), or if there's an actual statistical difference between the two.

I'm making this issue to keep track of it; will have a look at it soon.

[yebai]

We have tools for validating MCMC samplers: https://github.com/TuringLang/MCMCDebugging.jl

[torfjelde]

Is that not mainly about reversibilty? Still useful, just seems somwhat limited in what we want to test?

There's also @Red-Portal 's https://krkim.me/MCMCTesting.jl/dev/ which also looks quite nice:) (though I believe it provides similar functionality as MCMCDebugging.jl)

EDIT: I see there's also some code for MMD calculations in MCMCDebugging.jl 👍

[torfjelde]

A note for the future: we can replace the missing arguments used in MCMCDebugging.jl with the following

rand_x_given_θ(model::DynamicPPL.Model, θ::AbstractVector{<:Real}; kwargs...) = rand_x_given_θ(Random.default_rng(), model, θ; kwargs...)
function rand_x_given_θ(rng::Random.AbstractRNG, model::DynamicPPL.Model, θ::AbstractVector{<:Real}; varinfo=VarInfo(model))
    # Get the varinfo with the given values.
    varinfo = DynamicPPL.unflatten(varinfo, θ)
    # Fix the deconditioned model with the given values.
    θ_dict = DynamicPPL.values_as(varinfo, DynamicPPL.OrderedDict)
    model_fixed = DynamicPPL.fix(DynamicPPL.decondition(model), θ_dict)
    # Sample from the fixed model.
    return rand(rng, Vector, model_fixed)
end

rand_joint(model::DynamicPPL.Model; kwargs...) = rand_joint(Random.default_rng(), model; kwargs...)
function rand_joint(rng::Random.AbstractRNG, model::DynamicPPL.Model; varinfo=VarInfo(model))
    # Sample latent variables.
    θ = rand(Vector, model)
    # Sample observed variables.
    x = rand_x_given_θ(rng, model, θ; varinfo)
    return θ, x
end

Then it's just a matter of doing

@model function BetaBinomial()
    θ ~ Beta(2, 3)
    x ~ Binomial(3, θ)
    return θ, x
end

model = BetaBinomial() | (x = 3,)
rand_joint(model)
rand_x_given_θ(model, first(rand_joint(model))

which will allow us to perform this test on more models without too much hassle.