Conditionals überprüfen

[trappmartin]

Stimmt das: https://github.com/phipsgabler/LittleGibbs.jl/blob/4358b6a973a12b0546b31d511abb4bc0e1dae1c4/test/testmodels.jl#L48 ?

[trappmartin]

Aha, ich verstehe. Du machst das da Distributions.jl für eine gamma Verteilung einen scale statt einem rate parameter will.

[trappmartin]

Rechne bitte die conditionals einmal selber durch.

[phipsgabler]

Does the following make sense? It still fails the tests...

[trappmartin]

It looks in principle right. Could you please check if the sampling converges to a correct value if you keep one parameter fixed all time. I can do the derivation myself if you don’t find the problem.

[phipsgabler]

@trappmartin I corrected the value of alpha in

function gdemo_cond_λ(α₀, θ₀, x, m)
    N, x̄, s² = gdemo_statistics(x)
    αₙ = α₀ + (N - 1) / 2 + 1
    βₙ = (s² * N / 2 + m^2 / 2 + inv(θ₀))
    return Gamma(αₙ, inv(βₙ))
end

but the test for lambda still fails with

gdemo: Test Failed at /home/philipp/git/LittleGibbs/test/runtests.jl:55
  Expression: isapprox(λ̂, λ_true, atol=0.2, rtol=0.0)
   Evaluated: isapprox(0.9127713624702762, 0.4897959183673469; atol=0.2, rtol=0.0)

[trappmartin]

Can it be that something funny happens when you pass around the parameters of the prior? If I try the following on your master, then everything seems to work just fine.

using Test, LittleGibbs
include("test/testmodels.jl")

function test_gdemo(observations, α₀, θ₀, m_true, λ_true; m_init=0.0, λ_init=1.0)
    # Tests for expectations
    conditionals = gdemo(observations; α₀=α₀, θ₀=θ₀)

    model = BlockModel((m=conditionals[:m], λ=(m,) -> Normal(λ_true, 0)))
    sampler = Gibbs((m=m_init, λ=λ_init))
    chain = sample(model, sampler, 10_000)
    chain_m = map(t -> t.params[:m], chain[1_000:end])
    @test isapprox(mean(chain_m), m_true, atol=0.2, rtol=0.0)

    model = BlockModel((m=(λ,) -> Normal(m_true, 0), λ=conditionals[:λ]))
    sampler = Gibbs((m=m_init, λ=λ_init))
    chain = sample(model, sampler, 10_000)
    chain_λ = map(t -> t.params[:λ], chain[1_000:end])
    @test isapprox(mean(chain_λ), λ_true, atol=0.2, rtol=0.0)

    model = BlockModel(conditionals)
    sampler = Gibbs((m=m_init, λ=λ_init))

    chain = sample(model, sampler, 10_000)
    chain_m = map(t -> t.params[:m], chain[1_000:end])
    chain_λ = map(t -> t.params[:λ], chain[1_000:end])

    @test isapprox(mean(chain_m), m_true, atol=0.2, rtol=0.0)
    @test isapprox(mean(chain_λ), λ_true, atol=0.2, rtol=0.0)
end

N = 1_000
α₀′, θ₀′ = 1.5, 1.5
λ_true = rand(Gamma(α₀′, θ₀′))
σ_true = √(1 / λ_true)
m_true = rand(Normal(0, σ_true))
test_gdemo(rand(Normal(m_true, σ_true), N), α₀′, θ₀′, m_true, λ_true)

[trappmartin]

julia> test_gdemo(rand(Normal(m_true, σ_true), N), α₀′, θ₀′, m_true, λ_true)
Test Passed

[phipsgabler]

That happens most of the time, but not always. Sometimes even the test for the conditional fails.

[trappmartin]

But I’m still surprised it actually does work. And in my case it always did, but that might have been luck. I tried it several times. However, the one in runtest did not work for me.

Any idea why this one works and the one in runtest not?

[phipsgabler]

It works for me now in the GibbsConditionals PR. Hopefully just a numerics problem.