Understanding the Delayed Rejection Metropolis algorithm (Mira, 2001a)

I'm having trouble understanding the algorithm as briefly described here, and I can't find the original paper by Mira since it seems to be from some obscure print journal (Metron Volume 59).

The first stage is typical, $$\alpha_1(x,y) = \min\left(1,\frac{N_1}{D_1}\right)$$ Where, $$N_1 = \pi(y)q_1(y,x)$$ $$D_1 = \pi(x)q_1(x,y)$$ $$x = \textrm{current value}$$ $$y = \textrm{proposed value}$$ $$\pi = \textrm{target distribution}$$ $$q_1(x,\cdot) = \textrm{distribution from which y is drawn}$$

The second stage is

$$\alpha_2(x,y,z) = \min\left(1,\frac{N_2}{D_2}\right)$$

Where $$N_2 = \pi(z)q_1(z,y)q_2(z,y,x)[1-\alpha_1(z,y)]$$ $$D_2 = \pi(x)q_1(x,y)q_2(x,y,z)[1-\alpha_1(x,y)]$$

Thanks to @Xi'an, I realized that I have to pay attention to the distributions $q_1$ and $q_2$ even if I assume the proposals are symmetric, because $q_1(z,y)$ is likely not equal to $q_1(x,y)$. The second thing that I realized is that "the second stage candidate is computed so that reversibility of the Markov Chain...is preserved", so I still don't get the intuition but I understand how this form came about.

Questions:

If somebody has some intuition I would appreciate it. I'm playing around with toy distributions now, things like what if $z = x$ and what if $z = x\pm \epsilon$.

Also is $q_2(z,y,x)$ also not necessarily equal to $q_2(x,y,z)$ if we assume $q_1 \neq q_2$? The paper uses a smaller covariance for the second stage proposal. Is it like $q_1(y|z)q_2(x|y)$ vs. $q_1(y|x)q_2(z|y)$?

• Well I thought I already understood MH, I implemented what I thought was standard MH as well as the Adaptive Metropolis by Haario et al. and they both seemed to work well at sampling the posterior. My understanding of the standard MH is that you propose a parameter value, evaluate the posterior at that value, and then your probability of moving is min(1,p2/p1) where p2 is the posterior at the proposed value, p1 is the posterior at the current value. Am I wrong? This didn't seem that complicated when I read it on wikipedia/etc. Jul 22, 2018 at 20:11
• Hm yes this whole time I was assuming a symmetric proposal such as uniform or Gaussian which both seem typical. So it's unnecessary for the first stage but in the second stage that may not be the case, I'll reevaluate. What else do you think is necessary to understand the second stages and beyond? Is there some intuition or do you think I just need to spend more time looking at it? Jul 22, 2018 at 20:26
• Well so for the next stage, q1(z,y) will not equal q1(x,y), but q2(x,y,z)=q2(z,y,x) doesn't it? So that simplifies a bit. Jul 22, 2018 at 20:29

The validation of the delayed rejection algorithm is that, when starting with a realisation of the variable $X_t\sim\pi(x)$, the outcome of the Markov move to $X_{t+1}$ still remains distributed as $X_{t+1}\sim\pi(x)$. If the first step leads to an acceptance, the validation is the same as with the classical Metropolis-Hastings algorithm, hence using the same acceptance probability, $\alpha_1(x_t,\cdot)$ and the accepted value is distributed as $$q_1(x_t,y)\alpha(x_t,y)$$ If the value $Y$ at the first stage is rejected, it is distributed as $$q_1(x_t,y)[1-\alpha(x_t,y)]$$up to a normalising constant. The second stage value is then generated from the proposal $q_2(x_t,y,z)$ which means that $Z$ is generated conditional on $(x_t,y)$ [and hence that a better notation would be $q_2(z|x_t,y)$]. In this second stage, $y$ becomes an auxiliary variable, with the joint distribution of $(Y,Z)$ given $x_t$ and a first rejection being$$q_1(x_t,y)[1-\alpha(x_t,y)]q_2(x_t,y,z)$$which explains why this block appears in $\alpha_2(x_t,y,z)$. Another way to explain the representation of $\alpha_2(\cdot,\cdot,\cdot)$ is that the marginal distribution of $Z$ is (proportional to)$$\int_\mathcal{X} q_1(x_t,y)[1-\alpha(x_t,y)]q_2(x_t,y,z)\,\text{d}y$$ and that an alternative acceptance probability would be \begin{align*}\alpha_2^\prime(x_t,z)=\min\bigg\{1,\,&\pi(z)\int_\mathcal{X} q_1(z,y)[1-\alpha(z,y)]q_2(z,y,x_t)\,\text{d}y\\&\left.\Big/\pi(x_t)\int_\mathcal{X} q_1(x_t,y)[1-\alpha(x_t,y)]q_2(x_t,y,z)\,\text{d}y\right\}\end{align*} But this probability cannot be computed in most cases [except when $q_2(x_t,y,z)=q_2(x_t,z)$ and hence is replaced by a valid ratio of unbiased estimators of the integrals as in Andrieu & Roberts (2009) pseudo-marginal technique.
This is for the "intuition" part, but the validation proceeds by detailed balance, namely that the distribution of $X_{t+1}$ given $X_t=x_t$ writes as $$q(x_t,x)\alpha_1(x_t,x)q_1(x_t,x)+\int_\mathcal{X} \alpha_2(x_t,y,x)\,q_1(x_t,y)[1-\alpha(x_t,y)]q_2(x_t,y,x)\,\text{d}y+\int_\mathcal{X^2} \alpha_2(x_t,y,x)\,q_1(x_t,y)[1-\alpha(x_t,y)]q_2(x_t,y,z)\,\text{d}y\,\text{d}z\,\delta_{x_t}(x)$$and that this density satisfies $$\pi(x_t)q(x_t,x)=\pi(x)q(x,x_t)$$