We would like to announce our package together with our new paper “Sticky PDMP samplers for sparse and local inference problems” (Arxiv link)
Sleek implementations of the ZigZag, Boomerang and other assorted piecewise deterministic Markov processes for Markov Chain Monte Carlo including Sticky PDMPs for variable selection
ZigZagBoomerang.jl contains a user-friendly and clean implementations of the most prominent piecewise deterministic Monte Carlo methods. These Monte Carlo methods are fast (non-reversible, with momentum) and allow for subsampling of data without creating bias and can be used for sparse and local inference problems.
Within a few lines of code, you can run challenging problems with large data size and in high dimensional space. The only ingredient the user has to add is (the gradient of ) the target log-density and perhaps a rough estimate of the posterior mean.
For example the gradient of a Gaussian density is linear and can be given by a sparse matrix
Γ and then it is
# Γ = ... # sparse precision matrix # Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i] ∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # more efficient that dot(Γ[:, i], x) # Random initial values t0 = 0.0 x0 = randn(n*n) θ0 = rand([-1.0,1.0], n*n) # Rejection bounds c = 1.0*ones(length(x0)) # Define ZigZag Z = ZigZag(Γ, x0*0) # Run sparse ZigZag for T time units and collect trajectory T = 20.0 @time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ; adapt=true) @time traj = collect(discretize(trace, 0.1))
We have written ZigZagBoomerang.jl to be used as sampling engine for probabilistic programming languages. Let’s demo this with Soss.jl where work on Zig-Zag integration is already under way:
using Soss m = @model x begin α ~ Normal() β ~ Normal() yhat = α .+ β .* x y ~ For(eachindex(x)) do j Normal(yhat[j], 2.0) end end # generate data from model x = randn(20); obs = -0.1 .+ 2x + 1randn(20); posterior = m(x=x) | (y=obs,)
using ZigZagBoomerang include("zigzag.jl") T = 100.0 trace, final, (num, acc) = @time zigzag(posterior, T) # sample Soss posterior # trace is a continous object, discretize to obtain samples ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1)) ts # time points xs # vector of points
(see also zigzag.jl)
using Plots p = plot(ts, getindex.(xs, 1)) plot!(ts, getindex.(xs, 2), color=:red) p2 = plot(first.(xs), last.(xs)) p