RS after MCMC

mcmc.jl

using Distributions

K = 10^5
μ1 = 1.0
μ2 = -1.0
σ1 = σ2 = 0.3
args = (μ1 = μ1, σ1 = σ1, μ2 = μ2, σ2 = σ2)

lp(x, args) = log(pdf(Normal(args.μ1, args.σ2), x) + pdf(Normal(args.μ2, args.σ2), x))


function mcmc!(xs, x, (lp, args), Q, s, S, verbose = false)
    while s < last(S)
        xᵒ = rand(Q(x))
        lα = lp(xᵒ, args) - lp(x, args) + logpdf(Q(xᵒ), x) - logpdf(Q(x), xᵒ)
        verbose && print("$lα $xᵒ")
        if rand() < exp(lα)
            x = xᵒ
        verbose &&    print("✓")
        end
        verbose && println()
        s += 1
        if s in S
            push!(xs, x)
        end
    end
    xs, x
end

xs, _ = mcmc!(Float64[], 0.0, (lp, args), x->Normal(x, 0.1), 0, 0:K)
std(xs)


function wmcmc!(xs, x, (lp, args), Q, s, S; weight=(x->1.0, 1.0), verbose = false)
    f, fmin = weight
    while s < last(S)
        xᵒ = rand(Q(x))
        lα = log(f(xᵒ)) - log(f(x)) + lp(xᵒ, args) - lp(x, args) + logpdf(Q(xᵒ), x) - logpdf(Q(x), xᵒ)
        verbose && print("$lα $xᵒ")
        if rand() <= exp(lα)
            x = xᵒ
        verbose &&    print("✓")
        end
        verbose && println()
        s += 1
        if s in S
            if rand() <= fmin/f(x)
                push!(xs, x)
            end
        end
    end
    xs, x
end
xs2, _ = wmcmc!(Float64[], 0.0, (lp, args), x->Normal(x, 0.1), 0, 0:K, weight=(x->0.1 + pdf(Normal(0.0, 0.5),x), 0.1))