Gibbs sampling with mixed prior using a Metropolis-Hastings step

Question

My questions are about a sampling procedure for fitting a Bayesian hierarchical model where one of the priors is a mixture distribution of discrete and continuous parts. The model is not my own but I am trying to implement it and fit it to some data. The prior for one of the variables θ is of the form:

$$ \theta \sim \rho I_{[\theta = 0]} + (1 - \rho) N(\mu_\theta, \sigma^2_\theta) $$ where $\rho$ is the (prior) probability that the variable is 0 and $I_{[\theta = 0]}$ is the indicator function. $\theta$ is a mixture of a discrete part (a "point-mass" at 0) and a normal distribution.

The joint posterior distribution can be written down and the complete conditional distributions may then be read off (many of them are conjugates) for use in a Gibbs sampler. What I'm interested in is the complete conditional for $\theta$. Like the prior it contains a "point-mass". It is proportional to:

$$f(\theta) = \frac{[\exp(\gamma + \theta)]^y}{(1 + \exp(\gamma + \theta))^{N_t}} \left[\rho I_{[\theta = 0]} + (1 - \rho)\frac{1}{\sqrt{2\pi\sigma^2_\theta}} \exp\left[-\frac{1}{2\sigma^2_\theta}(\theta - \mu_\theta)^2 \right]\right]$$

The $\gamma, \rho, \mu_\theta, \sigma^2_{\theta}$ parameters change value on each iteration of the gibbs-sampler. As I only know the distribution up to a constant I can't just use a straightforward sampling approach. At the moment I use a Metropolis-Hastings (MH) step as follows:

Proposal distibution: $q(\theta|x) = 0.5 I_{[\theta = 0]} + 0.5 N\left(x, \sigma^2_{MH}\right)$, where

$$q(\theta|x) = \left\{ \begin{array}{lr} 0.5 & \theta = 0 \\ 0.5 \; g(\theta | x, \sigma^2_{MH}) & \theta \ne 0 \end{array} \right .$$

where $g$ is the density of $N\left(x, \sigma^2_{MH}\right)$, $\sigma^2_{MH}$ a fixed value.

If $\theta_{curr}$ is the current value, the MH step is then:

Simulate a candidate value $\theta^*$ from $q(\theta|\theta_{curr})$. Form the ratio: $r = \frac{f(\theta^*)q(\theta_{curr}|\theta^*)}{f(\theta_{curr})q(\theta^*|\theta_{curr})}$

Simulate $u$ from a uniform and accept $\theta^*$ with probability $\min(1, r)$.

My questions are as follows:

1. Is the approach correct - i.e. is it theoretically valid?
2. Is there a better or more "standard" approach to this type of sampling situation that I should consider?
3. Is it possible to assess this type of sampling where there is a move from a discrete to a continuous distribution as part of the step in term of acceptance rates etc?

Apologies if this has been covered in another question but I haven't been able to locate an adequate anwser and this type of scenario seems to be rarely discussed in the general books about MCMC sampling that I've looked at. I realise that there are probably no simple answers here but any pointers would be appreciated.

Answer 1

1. Your approach is correct because your proposal is absolutely continuous wrt the dominating measure, which is the sum of the Dirac mass at zero and of the Lebesgue measure and because it is everywhere positive. Just to make things precise, when $\theta=0$, we have $$f(0) = \frac{[\exp(\gamma + \theta)]^y}{(1 + \exp(\gamma +\theta))^{N_t}}\rho I_{[\theta = 0]}$$and not$$f(0) = \frac{[\exp(\gamma + \theta)]^y}{(1 + \exp(\gamma + \theta))^{N_t}} \left[\rho I_{[\theta = 0]} + (1 - \rho)\frac{1}{\sqrt{2\pi\sigma^2_\theta}} \exp\left[-\frac{1}{2\sigma^2_\theta}\mu_\theta^2 \right]\right]$$since the Lebesgue measure gives no weight to $\theta=0$. The practical difficulty of the proposal is to achieve proper exploration of both sides of the target density, the Dirac mass and the continuous part, which means calibrating properly the weights of those sides in the proposal. This means running preliminary versions with a given weight, 0.5 for instance, and monitor the number of accepted values with $\theta=0$ and with $\theta\ne 0$. If there is a strong disproportion, this may be due to the choice of the weights or to the asymmetry between both sides in the posterior. If one of the two sides is never visited, more asymmetrical weights should be used.

2. An alternative to this global exploration is to use one version or another of reversible jump MCMC introduced by Green (1995). Without getting into the technical details, reversible jump MCMC considers two statistical models, one with $\theta=0$ and one with $\theta\ne 0$. And it proposes reversible switches between both models at each MCMC iteration. In your case, the calibration should be relatively easy. A slightly simpler version is Carlin and Chib (1995) saturation, where $\theta$ is simulated at each step, but either from the prior or from a pseudo-posterior.

3. I do not understand this part. The problem is standard once you realise the dominating measure is the sum of the Dirac mass at zero and of the Lebesgue measure.