gdemo_test.jl

using Distributions

α₀ = 2.0
θ₀ = inv(3.0)

x = [1.5, 2.0]

function gdemo_statistics(x)
    # The conditionals and posterior can be formulated in terms of the following statistics:
    N = length(x) # number of samples
    x̄ = mean(x) # sample mean
    s² = var(x; mean=x̄, corrected=false) # sample variance
    return N, x̄, s²
end

function gdemo_cond_m(c)
    N, x̄, s² = gdemo_statistics(x)
    mₙ = N * x̄ / (N + 1)
    λₙ = c.λ * (N + 1)
    σₙ = √(1 / λₙ)
    return Normal(mₙ, σₙ)
end

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

@model gdemo(x) = begin
    λ ~ Gamma(α₀, θ₀)
    m ~ Normal(0, √(1 / λ))
    x .~ Normal(m, √(1 / λ))
end

m = gdemo(x)

# julia> sample(m, Gibbs(Turing.GibbsConditional(:λ, gdemo_cond_λ), Turing.GibbsConditional(:m, gdemo_cond_m)), 10000)
# Sampling 100%|███████████████████████████████████████████████████████████████| Time: 0:00:01
# Object of type Chains, with data of type 10000×3×1 Array{Float64,3}

# Iterations        = 1:10000
# Thinning interval = 1
# Chains            = 1
# Samples per chain = 10000
# internals         = lp
# parameters        = m, λ

# 2-element Array{ChainDataFrame,1}

# Summary Statistics
#   parameters    mean     std  naive_se    mcse        ess   r_hat
#   ──────────  ──────  ──────  ────────  ──────  ─────────  ──────
#            m  1.1682  0.9812    0.0098  0.0091  9641.3892  0.9999
#            λ  0.6345  0.4292    0.0043  0.0047  9630.6857  0.9999

# Quantiles
#   parameters     2.5%   25.0%   50.0%   75.0%   97.5%
#   ──────────  ───────  ──────  ──────  ──────  ──────
#            m  -0.8233  0.6486  1.1648  1.6960  3.1134
#            λ   0.0986  0.3187  0.5384  0.8478  1.7164


# julia> 7/6, 24/49  # expected values
# (1.1666666666666667, 0.4897959183673469)