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:


*

*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.

*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?

*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?

 A: 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}. 
$$
