Primitive Gibbs sampler, with two demos.

gibbs.jl

using Distributions


function mysampler(previous::Union{NamedTuple, Nothing})::NamedTuple
    µ = rand(Normal(3, 2))
    σ = rand(TruncatedNormal(1, 3, 0, Inf))
    (µ=µ, σ=σ)
end

function myloglikelihood(θ::NamedTuple)::Float64
    obs = [3.45, 7.37, 10.08, 4.33, 3.61, 4.23, 8.35, 5.74, 0.77, 4.28]
    sum(logpdf.(Normal(θ[:µ], θ[:σ]), obs))
end

function mygdemo_sampler(previous::Union{NamedTuple, Nothing})::NamedTuple
    s² = rand(InverseGamma(2, 3))
    m = rand(Normal(0, sqrt(s²)))

    (s²=s², m=m)
end

function mygdemo_loglikelihood(θ::NamedTuple)::Float64
    x = 1.5
    y = 2.
    sum(logpdf.(Normal(θ[:m], sqrt(θ[:s²])), [x,y]))
end

function sample(sampler::F, loglikelihood::G; n=1000)::Vector{NamedTuple} where {F<:Function, G<:Function}
    previous = sampler(nothing)
    U = Uniform(0, 1)
    samples = []
    i = 1
    while i <= n
        next = sampler(previous)
        ratio = min(1., exp((loglikelihood(next)-loglikelihood(previous))))
        # @show (loglikelihood(next), loglikelihood(previous), ratio)
        if rand(U) < ratio
            i += 1
            push!(samples, next)
            previous = next
        end
    end
    samples
end

# Extract samples using, e.g., µs = [i[:µ] for i in sample(...)]

turing_comparison.jl

using Turing

X = [3.45, 7.37, 10.08, 4.33, 3.61, 4.23, 8.35, 5.74, 0.77, 4.28]

@model function gibbs_analog(obs)
    mu ~ Normal(3, 2)
    sigma ~ TruncatedNormal(1, 3, 0, Inf)
    for x in eachindex(obs)
        obs[x] ~ Normal(mu, sigma)
    end
end

@model function gdemo(x, y)
  s² ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s²))
  x ~ Normal(m, sqrt(s²))
  y ~ Normal(m, sqrt(s²))
end