Gibbs sampling from a complex full conditional

Question

I have a sampling question relating to Gibbs sampling of a complicated full conditional. Supposed I have a complicated full conditional that I want a single sample from $p(\theta_i$|$\theta_{-i}$, $data)$ that cannot be sampled from directly.
I have a three fold question:

1. At this point would one use a Metropolis-Hastings step?

2. Do we accept the first MH step that yields a new $\theta_i$ or is it typically better to use a "burn in" phase before accepting the new $\theta_i$. If so typically how many steps are required? How do you fine tune a good jumping distribution in this case?

3. If the MH step rejects the sample, do I then keep the new $\theta_i$ the same or keep going?

Thank you!

Answer 1

Here is an excerpt from our Monte Carlo Statistical Methods book:

10.3.3. Metropolizing the Gibbs Sampler

Hybrid MCMC algorithms are often useful at an elementary level of the simulation process; that is, when some components of the Gibbs sampler conditionals cannot be easily simulated. Rather than looking for a customized algorithm such as Accept--Reject in each of these cases or for alternatives to Gibbs sampling, there is a compromise suggested by Müller (1991), sometimes called ``Metropolis-within-Gibbs". In any step i of the Gibbs sampling algorithm with a difficult simulation from the full conditional $g_i(\theta_i|\theta_j,j\neq i)$, substitute a simulation from an instrumental distribution $q_i$. This means replacing step $i$ of the Gibbs sampler

i. Simulate $$ {\tilde \theta}_i \sim g_i(\theta_i^{(t)}|\theta_1^{(t+1)},\ldots,\theta_{i-1}^{(t+1)},\theta_{i+1}^{(t)},\ldots,\theta_p^{(t)}) $$

with

i.1. Simulate $$ {\tilde \theta}_i \sim q_i(\theta_i|\theta_1^{(t+1)},\ldots,\theta_i^{(t)},\theta_{i+1}^{(t)}, \ldots,\theta_p^{(t)}) $$ i.2. Take $$ \theta_i^{(t+1)} = \begin{cases} y_\theta^{(t)} & \hbox{ with probability }1-\rho, \cr {\tilde \theta}_i & \hbox{ with probability }\rho, \cr \end{cases} $$ where $$ \rho = 1 \wedge \left\{ \frac{\left( \displaystyle{{g_i({\tilde \theta}_i|θ_1^{(t+1)},\ldots,θ_{i-1}^{(t+1)},θ_{i+1}^{(t)}, \ldots,θ_p^{(t)}) \over q_i({\tilde θ}_i|θ_1^{(t+1)},\ldots,θ_{i-1}^{(t+1)},θ_i^{(t)},θ_{i+1}^{(t)}, \ldots,θ_p^{(t)})} }\right)}{\left( \displaystyle{ {g_i(θ_i^{(t)}|θ_1^{(t+1)},\ldots,θ_{i-1}^{(t+1)},θ_{i+1}^{(t)}, \ldots,θ_p^{(t)}) \over q_i(θ_i^{(t)}|θ_1^{(t+1)},\ldots,θ_{i-1}^{(t+1)},{\tilde θ}_i,θ_{i+1}^{(t)}, \ldots,θ_p^{(t)})}}\right)} \right\} \;. $$

An important point about this substitution is that the above Metropolis-Hastings step is only used once in an iteration from the Gibbs sampler. The modified step thus proposes a single simulation ${\tilde \theta}_i$ instead of trying to approximate $g_i(θ_i|θ_j,j\neq i)$ more accurately by producing $T$ simulations from $q_i$. The reasons for this choice are twofold:

First, the resulting hybrid algorithm is valid since $g$ is its stationary distribution.

Second, Gibbs sampling also leads to an approximation of $g$. To provide a more "precise" approximation of $g_i(y_i|y_j,j\neq i)$ does not necessarily lead to a better approximation of $g$ and the replacement of $g_i$ by $q_i$ may even be beneficial for the speed of excursion of the chain on the surface of $g$.

You may also take a look at this other X validated question.