Gibbs sampling from a complex full conditional 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:


*

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

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

*If the MH step rejects the sample, do I then keep the new $\theta_i$ the same or keep going?
Thank you!
 A: 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.
