# Gibbs sampling for spike and slab priors

In Spike and slab variable selection (equation 4) there is a model setup of the form

$\beta_k | \lambda_k, \tau_k \sim \text{Normal} (0, \lambda_k \tau_k^2)$

$\lambda_k | \nu_0, w \sim (1-w)\delta_{\nu_0}(\lambda_k) + w \delta_1(\lambda_k)$

where $\beta_k$ is the $k^{th}$ regression coefficient and $\delta_x$ is the dirac-delta function centred at $x$ (I have changed notation slightly).

I'm trying to derive a Gibbs sampler for a similar model. The Gibbs sampler for this algorithm is in the appendix of the above link (page 43). My confusion comes from the update for $\lambda_k$:

$p(\lambda_k | \cdot) \propto p(\beta_k | 0, \lambda_k \tau_k^2) p(\lambda_k | \nu_0, w)$

which if you follow through gives you an unnormalised density of the form

$\frac{1}{\sqrt{\lambda_k \tau_k^2}} \exp(-\frac{2}{\lambda_k \tau_k^2}\beta_k^2)[(1 - w) \delta_{\nu_0}(\lambda_k) + w \delta_1(\lambda_k)]$

Intuitively, I can see how multiplying the exponent factor with the first term gives a point mass at $\nu_0$ and with the second gives a point mass at $1$, which we then normalise to give the Gibbs update in the attached paper (ie all the $\lambda_k$s in the above equation get set to either $\nu_0$ or $1$ for the update). However, I feel some things don't entirely make sense:

1. Dirac-delta functions "pick out" the point mass values when integrating over the region around the point mass, but there is no such integration here.
2. How does one sample from such a conditional distribution anyway? Is it simply the weighted average of the two point masses, or one or the other point mass with probabilities given by the weights?
3. If it is the weighted average, isn't this similar to ARD rather than spike-and-slab since we're back at a continuous measure of sparsity?
• Having thought some more I think I can answer (2) myself - A distribution of point masses means it can only take one of the two values $\nu_0$ or $1$, ie you're sampling from a discrete distr. This makes qn 3 irrelevant. Oct 7 '16 at 14:13

The notation in the paper uses $$\mathcal J_k$$ instead of $$\lambda_k$$. I am going to use $$\lambda_k$$ as in the question. I am going to drop subscript $$k$$ for simplicity. The model is then \begin{align*} \beta \mid \lambda &\sim N(0, \lambda \tau^2) \\ \lambda &\sim (1-w) \delta_{\nu_0} + w \delta_1. \end{align*} The rest requires some measure theory knowledge and the Radon-Nikodym theorem. The trick to write a joint density for this model is to note that both of the point mass measures $$\delta_{\nu_0}$$ and $$\delta_1$$ have densities with w.r.t. $$\mu =$$ the counting measure on $$\{\nu_0, 1\}$$. With some abuse of notation let us write $$\delta_{\nu_0}(\lambda)$$ for the density of $$\delta_{\nu_0}$$ w.r.t. $$\mu$$ as well, and similarly $$\delta_1(\lambda)$$ denotes the density of $$\delta_1$$ w.r.t. $$\mu$$. It is easy to verify that $$\delta_{\nu_0}(\lambda) = \begin{cases} 1 & \lambda = \nu_0 \\ 0 & \lambda \neq \nu_0 \end{cases}$$ and similarly for $$\delta_1(\lambda)$$. Let $$\mathcal L$$ be the Lebesgue measure on the real line. Then the distribution of $$(\beta, \lambda)$$ is absolutely continuous w.r.t. $$\mathcal L + \mu$$ with density \begin{align*} p(\beta, \lambda) &\propto \underbrace{\frac{1}{\sqrt{\lambda \tau^2}} \exp \Bigl( - \frac{\beta^2}{2 \lambda \tau^2} \Bigr)}_{:= f(\lambda,\beta)}\cdot \bigl[(1-w) \delta_{\nu_0}(\lambda) + w \delta_1(\lambda)\bigr] \\ &= (1-w) f(\lambda, \beta) \delta_{\nu_0}(\lambda) + w f(\lambda, \beta) \delta_1(\lambda) \\ &= (1-w) f(\nu_0, \beta) \delta_{\nu_0}(\lambda) + (1-w) f(1, \beta) \delta_1(\lambda) \end{align*} where the last line follows since $$\delta_{\nu_0}(\lambda)$$ and $$\delta_1(\lambda)$$ are indicator functions.
Then, we have $$p(\lambda \mid \beta) \propto \underbrace{(1-w) f(\nu_0, \beta)}_{w_1} \delta_{\nu_0}(\lambda) + \underbrace{(1-w) f(1, \beta)}_{w_2} \delta_1(\lambda)$$ which is a density w.r.t. to $$\mathcal L + \mu$$. Since $$\beta$$ is a constant here, this is a discrete distribution taking values $$\nu_0$$ and $$1$$ with probabilities proportional to $$w_1$$ and $$w_2$$. This can be alternatively stated as $$\mathbb P(\lambda = \nu_0 \mid \beta) = \frac{w_1}{w_1 + w_2}, \quad \mathbb P(\lambda = 1 \mid \beta) = \frac{w_2}{w_1 + w_2}.$$