[ANN] ZigZagBoomerang.jl

So what about sparsity?

The Zig-Zag allows (that is what our paper is about) spike and slab priors for variable selection. Let’s demo this with Soss:

Let’s use MeasureTheory.jl 's SpikeMixture (a very spiky spike and slab with a Dirac measure as spike) to indicate that we have ground to believe the coefficients are actually zero:

m = @model x begin
    α ~ SpikeMixture(Normal(), 0.2) # 0.2*Normal() + 0.8*Dirac(0)
    β ~ SpikeMixture(Normal(), 0.2)
    yhat = α .+ β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end


x = randn(20);
obs = -0.1 .+ 2x + 1randn(20); 
T = 100.0
posterior = m(x=x) | (y=obs,)

sampling with the sticky ZigZag


include("sparsezigzag.jl")
trace, final, (num, acc) = @time sparse_zigzag(posterior, T, c=50)

ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1))
xs = xform(posterior).(xs) # make named tuples `NamedTuple{(:β, :α)`

now we can answer for example what is the posterior probability that the parameters are zero

julia> mean(getindex.(xs, :α) .== 0)
0.7211394302848576

julia> mean(getindex.(xs, :β) .== 0)
0.028985507246376812

The nice thing is that this scales very well…

4 Likes