Julia adaptation Mochihashi, Daichi (2020) Unbounded Slice Sampling

SliceSampling.jl

# Mochihashi, Daichi (2020) Unbounded Slice Sampling
# https://arxiv.org/abs/2010.01760
using Distributions
using Random
using Plots# Mochihashi, Daichi (2020) Unbounded Slice Sampling
# https://arxiv.org/abs/2010.01760
module Slice
# Unbounded slice sampling
function slice1(x, likfun; A=100.0, maxiter=1000)
    # shrink/expand transformation
    shrink(x) = 1 / (1 + exp(-x / A))
    expand(p) = -A * log(1 / p - 1)

    st, ed = 0.0, 1.0
    r = shrink(x)
    slice = likfun(x) - log(A * r * (1-r)) + log(rand())
    for iter in 1:maxiter
        rnew = st + (ed - st) * rand()
        xnew = expand(rnew)
        newlik = likfun(xnew) - log(A * rnew * (1-rnew))
        if newlik > slice
            return xnew
        elseif rnew > r
            ed = rnew
        elseif rnew < r
            st = rnew
        else
            return xnew
        end
    end
    @warn "slice_sample: max iteration $maxiter reached"
    return x
end


# Slice sampling restricted to positive reals
function slice1_positive(x, likfun; maxiter=1000)
    # shrink/expand transformation for x>0
    shrinkp(x) = x / (1 + x)
    expandp(p) = p / (1 - p)

    st, ed = 0.0, 1.0
    r = shrinkp(x)
    slice = likfun(x) - 2 * log(1 - r) + log(rand())

    for iter in 1:maxiter
        rnew = st + (ed - st) * rand()
        xnew = expandp(rnew)
        newlik = likfun(xnew) - 2 * log(1 - rnew)
        if newlik > slice
            return xnew
        elseif rnew > r
            ed = rnew
        elseif rnew < r
            st = rnew
        else
            return xnew
        end
    end
    @warn "slice_sample_positive: max iteration $maxiter reached"
    return x
end

function slicesampler(Iter, likfun)
x = 0.0
rv = zeros(Iter)
for i in 1:Iter
    x = slice1(x, likfun)
    rv[i] = x
end
return rv    
end

function slicesampler_pos(Iter, likfun)
xp = 1.0
rv = zeros(Iter)
for i in 1:Iter
    xp = slice1_positive(xp, likfun)
    rv[i] = xp
end
return rv   
end

end

using Plots
using Distributions
using Random
using QuadGK

likfun(x) = ( − (x−3)*(x−1))
rv = Slice.slicesampler(10000, likfun)
#無理やり正規化
lb = minimum(rv)
ub = maximum(rv)
Z, _ = quadgk(x ->  exp(likfun(x)), lb, ub)
histogram(rv, normalize=:pdf, label=false, color="lightgrey")
plot!(x-> exp(likfun(x))/Z, label=false, color="royalblue", linewidth=2)

##正のみ
likfun_pos(x) = 4*log(x) − x 
rv = Slice.slicesampler_pos(10000, likfun_pos)

#無理やり正規化
lb = minimum(rv)
ub = maximum(rv)
Z, _ = quadgk(x ->  exp(likfun_pos(x)), lb, ub)

histogram(rv, normalize=:pdf, label=false, color="lightgrey")
plot!(x-> exp(likfun_pos(x))/Z, label=false, color="royalblue", linewidth=2)


# Unbounded slice sampling
function slice_sample(x, likfun; A=100.0, maxiter=1000)
    # shrink/expand transformation
    shrink(x) = 1 / (1 + exp(-x / A))
    expand(p) = -A * log(1 / p - 1)

    st, ed = 0.0, 1.0
    r = shrink(x)
    slice = likfun(x) - log(A * r * (1-r)) + log(rand())
    for iter in 1:maxiter
        rnew = st + (ed - st) * rand()
        xnew = expand(rnew)
        newlik = likfun(xnew) - log(A * rnew * (1-rnew))
        if newlik > slice
            return xnew
        elseif rnew > r
            ed = rnew
        elseif rnew < r
            st = rnew
        else
            return xnew
        end
    end
    @warn "slice_sample: max iteration $maxiter reached"
    return x
end


# Slice sampling restricted to positive reals
function slice_sample_positive(x, likfun; maxiter=1000, args...)
    # shrink/expand transformation for x>0
    shrinkp(x) = x / (1 + x)
    expandp(p) = p / (1 - p)

    st, ed = 0.0, 1.0
    r = shrinkp(x)
    slice = likfun(x, args...) - 2 * log(1 - r) + log(rand())

    for iter in 1:maxiter
        rnew = st + (ed - st) * rand()
        xnew = expandp(rnew)
        newlik = likfun(xnew, args...) - 2 * log(1 - rnew)
        if newlik > slice
            return xnew
        elseif rnew > r
            ed = rnew
        elseif rnew < r
            st = rnew
        else
            return xnew
        end
    end
    @warn "slice_sample_positive: max iteration $maxiter reached"
    return x
end


# 例: 標準正規分布からサンプリング
likfun(x) = -0.5 * x^2  # log-likelihood（定数項は省略）
x = 0.0
Iter = 10000
rv = zeros(Iter)
for i in 1:Iter
    x = slice_sample(x, likfun)
    rv[i] = x
end

histogram(rv, normalize=:pdf, label=false, color="lightgrey")
plot!(x->pdf(Normal(),x), label=false, color="royalblue", linewidth=2)

# 例: 正のみに制約
likfun_pos(x) = -x  # e^{-x}, x>0 の密度の log
xp = 1.0
Iter = 10000
rv = zeros(Iter)
for i in 1:Iter
    xp = slice_sample_positive(xp, likfun_pos)
    rv[i] = xp
end

histogram(rv, normalize=:pdf, label=false, color="lightgrey")
plot!(x->pdf(Exponential(),x), label=false, color="royalblue", linewidth=2)