Metropolis-Hastings as a Soss model

mh.jl

using Soss, MeasureTheory, Plots

# ℓ: log-posterior density we want to sample
# propose: the Markov kernel for the proposal
# x: current point
mhstep = @model ℓ,propose,x begin
    y ~ propose(x)
    u ~ Uniform()
    accept = log(u) < ℓ(y) - ℓ(x)
    return accept ? y : x
end


mh(ℓ,propose,d0) = Chain(d0) do x
    mhstep(ℓ, propose, x)
end

# Make sure we can sample from a standard normal
ℓ(x) = -0.5x^2

# We'll propose with a normal random walk
propose(x) = Normal(x,0.1)

# Build the Markov chain
mc = mh(ℓ, propose, Uniform())

# Draw a sample
r = rand(mc)

# collect into a vector
x = collect(Iterators.take(r, 100000))


# plot(x)

# julia> mean(x)
# 0.05759233759871822

# julia> std(x)
# 0.964331244694455