# Question on how to model spike and slab priors

**URL:** <https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277>\
**Category:** Questions\
**Created:** [June 15, 2020, 6:29am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277 "2020-06-15T06:29:34Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![xing\_wu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/xing_wu/32/2816_2.png) [@xing\_wu](https://discourse.pymc.io/u/xing_wu)\
**Post date:** [June 15, 2020, 6:29am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/1 "2020-06-15T06:29:34Z")

</div>

Hi,  
I am new to the pymc3 community, and I have a question on how to apply spike and slab prior on variables.

Here is my model:

Y = X\beta + \epsilon

\beta\_i \sim (1-\gamma\_i)N(0,0.001^2) + \gamma\_i N(0,\sigma\_1^2)

\gamma\_i \sim \text{Bernoulli }(\pi)

\pi \sim \text{Beta}(1,1)

\epsilon \sim N(0,\sigma\_e^2)

\sigma\_1 \sim \text{HalfCauchy}(5)

\sigma\_e \sim \text{HalfCauchy}(5)

\beta and \gamma are vectors with length = n, in other words, length(beta) == length(gamma) ==n .

I cannot find examples showing mixture on the variables, and I am not sure what I wrote is correct or not.

here is my pymc3 model:

```auto
sigma_1 = pm3.HalfCauchy('sigma_1',5)
sigma_e= pm3.HalfCauchy('sigma_e',5)
pi = pm3.Beta('pi',1,1)
gamma = pm3.Bernoulli('gamma',p = pi,shape = n)
mixture_sd = pm3.math.switch(gamma > 0.5, sigma_1, 0.001)
beta = pm3.Normal('beta',mu = 0,sigma = mixture_sd, shape = n)
mu = tt.dot(X,beta)
likelihood = pm3.Normal('y',mu = mu, sigma = sigma_e,observed = y)

```

The mean(posterior) for \pi and \gamma are too small to be true. The program almost always reports that all of the \beta\_i comes from N(0,0.001^2).

I would really appreciate it if someone can help me on this.

Thanks

Xing

---

<div class="post-metadata">

**Author:** ![ld92](https://avatars.discourse-cdn.com/v4/letter/l/d9b06d/32.png) [@ld92](https://discourse.pymc.io/u/ld92)\
**Post date:** [June 15, 2020, 6:41am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/2 "2020-06-15T06:41:17Z")

</div>

Interesting, I posted a similar question yesterday:

> [@Spike and slab logsitic regression - marginalizing over bernoulli distribution](https://discourse.pymc.io/t/spike-and-slab-logsitic-regression-marginalizing-over-bernoulli-distribution/5273):
>
> Hello all, Im trying to implement a spike and slab prior for the beta coefficients in logistic regression. This works, but is very slow: with pm.Model() as uniform\_model: alfa = pm.Normal(‘alfa’, mu=0, sd=3) beta = beta\_spike\_slab(X.shape[1], 0.001) #spike/slab X\_shared = pm.Data(‘x\_shared’, X) mu = alfa + pm.math.dot(X\_shared,beta) θ = pm.Deterministic('θ', 1 / (1 + pm.math.exp(-mu))) y\_pred = pm.Bernoulli('y\_pred', p=θ, observed=y) trace = pm.sample(1000,tune=1000, cores=1, ta…

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!

---

<div class="post-metadata">

**Author:** ![xing\_wu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/xing_wu/32/2816_2.png) [@xing\_wu](https://discourse.pymc.io/u/xing_wu)\
**Post date:** [June 15, 2020, 6:56am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/3 "2020-06-15T06:56:15Z")

</div>

Thank you very much. In your code, did you mean

```auto
beta_slab = pm.Normal(‘beta_slab’, 0, 10, shape=shape)

```

In addition, did you spot anything wrong with my implementation?

---

<div class="post-metadata">

**Author:** ![ld92](https://avatars.discourse-cdn.com/v4/letter/l/d9b06d/32.png) [@ld92](https://discourse.pymc.io/u/ld92)\
**Post date:** [June 15, 2020, 9:36am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/4 "2020-06-15T09:36:51Z")

</div>

Yes, that is what I meant.  
I see youve updated your post with your model.

What are your sampling options? Because you cannot use NUTS for discrete variables, you need to (dramatically) increase your tuning steps. When I choose

`
trace = pm.sample(1000, tune=12000, nuts={"target_accept": 0.99})
`

with (this is adapted with your code)

`
```
sigma_1 = pm.HalfNormal('sigma_1',10, shape=shape)
pi = pm.Beta('pi',2,7)
gamma = pm.Bernoulli('gamma',p = pi,shape = shape)
mixture_sd = pm.Deterministic('mixture_sd', pm.math.switch(gamma > 0.5, sigma_1, 0.001))
beta = pm.Normal('beta',mu = 0,sigma = mixture_sd, shape = shape)

```
`

I do get the right beta coefficients with my toy dataset.

And I think the next step would be to not use match.switch, but replace it with a continous function like tt.net.sigmoid, to improve sampling, see [Frequently Asked Questions](https://discourse.pymc.io/t/frequently-asked-questions/74/3)

---

<div class="post-metadata">

**Author:** ![xing\_wu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/xing_wu/32/2816_2.png) [@xing\_wu](https://discourse.pymc.io/u/xing_wu)\
**Post date:** [June 15, 2020, 9:45pm UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/5 "2020-06-15T21:45:46Z")

</div>

I used NUTS for the continues variables and metropolis for gamma. I found it to be very slow.

---

<div class="post-metadata">

**Author:** ![Dushyant\_Sahoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/dushyant_sahoo/32/3839_2.png) [@Dushyant\_Sahoo](https://discourse.pymc.io/u/Dushyant_Sahoo)\
**Post date:** [March 26, 2021, 5:47pm UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/6 "2021-03-26T17:47:36Z")

</div>

A continuous version is much faster than the discrete one. Here is the code-

```auto
def conti_spike_slab(shape):
    mu_hat = 0
    sigma_hat = 10
    tau = 5
    lambda_hat = pm.Normal('lambda_hat', mu = mu_hat, sigma = sigma_hat, shape = shape)
    spike_raw = pm.Normal('spike_raw', mu = 0, sigma = 1, shape = shape)
    spike = pm.Deterministic('spike',tau*spike_raw*pm.invlogit(lambda_hat))
    return spike

```

Let me know what you guys think.

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [March 26, 2021, 6:27pm UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/7 "2021-03-26T18:27:52Z")

</div>

Perhaps the main point of the spike and slab prior is that it exactly zeros out many coefficients in the posterior samples. If you move to a continuous prior, you give that up. There are better continuous priors favoring many sparse entries. The horseshoe prior is one of them.

---

<div class="post-metadata">

**Author:** ![chartl](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chartl/32/1515_2.png) [@chartl](https://discourse.pymc.io/u/chartl)\
**Post date:** [March 26, 2021, 7:24pm UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/8 "2021-03-26T19:24:52Z")

</div>

Unless you’re drawing a single sample, having exact 0s is not particularly useful. And if you want the MAP why are you bothering with sampling at all? Otherwise the general process is “take the trace and ask how often this coefficient was 0” – which can be effectively read off from `lambda_hat` anyway.

---

<div class="post-metadata">

**Author:** ![Dushyant\_Sahoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/dushyant_sahoo/32/3839_2.png) [@Dushyant\_Sahoo](https://discourse.pymc.io/u/Dushyant_Sahoo)\
**Post date:** [March 26, 2021, 9:18pm UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/9 "2021-03-26T21:18:09Z")

</div>

I recently read a paper- [https://royalsocietypublishing.org/doi/pdf/10.1098/rsif.2018.0572](https://royalsocietypublishing.org/doi/pdf/10.1098/rsif.2018.0572) which claims that their continuous version of spike and slab outperforms horseshoe prior, it is based on iverse logit which I implemented above.

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [March 27, 2021, 5:03am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/10 "2021-03-27T05:03:26Z")

</div>

Interesting. I stand corrected with regard to my statement about the horseshoe - I look forward to taking a look at that paper.

---

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [March 27, 2021, 7:43am UTC](https://discourse.pymc.io/t/question-on-how-to-model-spike-and-slab-priors/5277/11 "2021-03-27T07:43:08Z")

</div>

That’s a really interesting paper, will have to try it later, thanks!