Spike-and-Slab Implementation using Turing.jl?

New to Julia
,
1

So I was looking for an implementation of spike-and-slab prior using the Turing.jl package. Say I have a model as follows:

y_{ig} \sim N(\theta_g, \sigma^2) where g = 1, \cdots, G and i = 1, \cdots, n_g independently
And I would like the priors on \theta_g to be independently a spike-and-slab as \pi\delta_0 + (1-\pi)t_v, which is a mixture of a degenerate mass at 0 and a central t distribution.

Can somebody give me some guidance as to how I can do this?

2

Distributions.jl has a Dirac distribution but it didn’t seem to work here. One option could be something like this?

using Turing
@model function M₁(g; G = maximum(g))
	ν ~ Exponential(10)
	π ~ Beta(1, 1)
	θ ~ filldist(MixtureModel([Normal(0, 0.001), TDist(ν)], [π, 1-π]), G)
	σ² ~ Exponential(1)
	y ~ MvNormal(θ[g], σ²*I)
end
3

While you can probably define a spike-and-slab prior somehow, it will be hard to sample from – due to its non-differentiable density at the location of the point measure. There are some modern differentiable alternatives/approximation available, such as the horseshoe prior. Betancourt has a very detailed discussion of sparsity priors and their properties.

4

Sorry for reviving this old topic. I tried to write a spike-and-slab prior for a homework assignment and the model below returns results that agree with a horseshoe implementation of the same model. I thought I’d add it for posterity:

using Turing, Statistics, LinearAlgebra
using Downloads, CSV, DataFrames

standardize(x) = (x .- mean(x)) ./ std(x)

data = let
	path = "https://hastie.su.domains/ElemStatLearn/datasets/prostate.data"
	file = Downloads.download(path)
	tmp = CSV.read(file, DataFrame; delim='\t', header=1, drop=[1])
	covariates = names(tmp)[1:8]
	transform(tmp, covariates .=> standardize => identity)
end

@model function turing_spikeslab(y, X)
	n, p = size(X)
	
	ω ~ Beta(1, 1)
	β₀ ~ Normal(0, 100)
	β ~ filldist(Normal(0, 1), p) 
	z ~ filldist(Bernoulli(ω), p)
	σ² ~ InverseGamma(3, 2)
	
	μ = β₀ .+ X * (β .* z)
	y ~ MvNormal(μ, σ² * I)
end

model = turing_spikeslab(data.lpsa, Matrix(data[:, 1:8]))
sampler = Gibbs(:z => PG(20), (:ω, :β₀, :β, :σ²) => NUTS())
chain = sample(model, sampler, MCMCThreads(), 1500, 4; discard_initial=500)