Interesting, I posted a similar question yesterday:
What you specify, can be pretty easily modeled in PyMC3:
def beta_spike_slab(shape,spike):
inclusion_prop = 0.05
beta_spike = pm.Normal(‘beta_spike’, 0, spike, shape=shape)
beta_slab = pm.Normal(‘beta_slab’, 10, shape=shape)
gamma = pm.Bernoulli(‘gamma’, inclusion_prop, shape=shape)
beta_spike_slab = pm.Deterministic(‘beta_spike_slab’,(beta_spike * (1-gamma)) + ((beta_slab * gamma)))
return beta_spike_slab
but this leads to slow sampling, because it is needed to marginalize over the bernoulli distribution.
I have not found a way to do this without divergence/sampling problems, so if anyone has any ideas, I would love to hear them!