Dear all!

We would like to announce our package together with our new paper “Sticky PDMP samplers for sparse and local inference problems” (Arxiv link)

# ZigZagBoomerang.jl

*Sleek implementations of the ZigZag, Boomerang and other assorted piecewise deterministic Markov processes for Markov Chain Monte Carlo including Sticky PDMPs for variable selection*

## 1. What is this about

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.

## 2. What can it do

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))
```

## 4. Using ZigZag and friends as backends in your probabilistic programming environment.

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)

Let’s plot

```
using Plots
p = plot(ts, getindex.(xs, 1))
plot!(ts, getindex.(xs, 2), color=:red)
p2 = plot(first.(xs), last.(xs))
p
```