How does Rao-Blackwellization of the Metropolis-Hastings algorithm work?

Question

I've read the paper A vanilla Rao-Blackwellization of Metropolis-Hastings algorithms, but I don't get what their actually suggested estimator is.

To give some detail, we are considering the following setup: Let

• $(E,\mathcal E,\lambda)$ be a measure space;
• $p:E\to[0,\infty)$ be $\mathcal E$-measurable with $$c:=\int p\:{\rm d}\lambda\in(0,\infty)$$ and $\mu$ denote the measure with density $\frac pc$ with respect to $\lambda$;
• $q:E^2\to[0,\infty)$ be $\mathcal E{\otimes2}$-measurable with $$c_x:=\int q(x,\;\cdot\;)\:{\rm d}\lambda\in(0,\infty)\;\;\;\text{for all }x\in E$$ and $Q(x,\;\cdot\;)$ denote the measure with density $\frac{q(x,\;\cdot\;)}{c_x}$ with respect to $\lambda$;
• $$\alpha(x,y):=\left.\begin{cases}\displaystyle\min\left(1,\frac{p(y)q(y,x)}{p(x)q(x,y)}\right)&\text{, if }p(x)q(x,y)\ne0\\1&\text{, otherwise}\end{cases}\right\}\;\;\;\text{for }x,y\in E;$$
• $$r(x):=1-\int Q(x,{\rm d}y)\alpha(x,y)\;\;\;\text{for }x\in E;$$
• $\delta_x$ denote the Dirac measure on $(E,\mathcal E)$ at $x\in E$ and $$\kappa(x,B):=\int_BQ(x,{\rm d}y)\alpha(x,y)+r(x)\delta_x(B)\;\;\;\text{for }(x,B)\in E\times\mathcal E;$$
• $(\Omega,\mathcal A,\operatorname P)$ be a probability space;
• $(X_n)_{n\in\mathbb N_0}$ denote the Markov chain generated by the Metropolis-Hastings algorithm with proposal kernel $Q$ and target distribution $\mu$, $(Y_n)_{n\in\mathbb N}$ denote the corresponding proposal sequence and $$Z_n:=(X_{n-1},Y_n)\;\;\;\text{for }n\in\mathbb N.$$

By construction, there is a $[0,1)$-valued independent identically $\mathcal U_{[0,\:1)}$-distributed process $(U_n)_{n\in\mathbb N}$ on $(\Omega,\mathcal A,\operatorname P)$ with

1. $(U_1,\ldots,U_n)$ and $(Z_1,\ldots,Z_n)$ are independent;
2. $$X_n=\begin{cases}Y_n&\text{if }U_n\le\alpha(Z_n);\\X_{n-1}&\text{otherwise}\end{cases}$$

for all $n\in\mathbb N$. Now, there is the following result:

$(Z_n)_{n\in\mathbb N}$ is a time-homogeneous Markov chain with transition kernel $$\kappa_{\text{aug}}((x,y),A\times B):=(1-\alpha(x,y))\delta_x(A)Q(x,B)+\delta_x(A)\alpha(x,y)Q(y,B)\;\;\;\text{for }x,y\in E\text{ and }A,B\in\mathcal E.$$ Moreover, $$\nu:=\mu\otimes Q$$ is invariant with respect to $\kappa_{\text{aug}}$.

This result is the justification why $$W_nf:=\sum_{i=1}^n\operatorname E\left[f(X_i)\mid Z_i\right]\;\;\;\text{for }n\in\mathbb N$$ is a consistent estimator of $\mu f$ for $f\in\mathcal L^1(\mu)$. Using this estimator was what I actually understood by "Rao-Blackwellization of the Metropolis-Hastings algorithm".

However, in the paper they do something else. They define $\tau_0:=0$ and $$\tau_k:=\inf\left\{n>\tau_{k-1}:U_n\le\alpha(Z_n)\right\}\;\;\;\text{for }k\in\mathbb N.$$ Moreover, $$\xi_k:=\begin{cases}\tau_k-\tau_{k-1}&\text{, if }\tau_k<\infty;\\0&\text{, otherwise}\end{cases}$$ and $$A_k:=\begin{cases}X_{\tau_k}&\text{, if }\tau_k<\infty;\\X_{\tau_{k-1}}&\text{, otherwise}\end{cases}$$ for $k\in\mathbb N_0$. Furthermore, let $$\varrho(x):=\int Q(x,{\rm d}y)\alpha(x,y)\;\;\;\text{for }x\in E.$$

For simplicitly, let's assume that $\tau_k<\infty$ for all $k\in\mathbb N_0$ (don't know how to generalize the following otherwise, but I also don't know whether this would be interesting at all).

We can now show that $(A_k)_{k\in\mathbb N_0}$ is a Markov chain with stationary distribution $\frac{p\varrho}{\tilde c}\lambda$, where $\tilde c:=\lambda(p\varrho)$. Moreover, $$\operatorname E\left[\xi_k\mid A_k\right]=\frac1\varrho(A_k)=\tilde c\underbrace{\frac p{\frac{p\varrho}{\tilde c}}(A_k)}_{=:\:w(A_k)}\tag1.$$

An estimator for $\mu f$ using the first $n$ distinct states of $X$ can now be written as $$\frac{\sum_{k=0}^{n-1}\xi_kf(A_k)}{\sum_{k=0}^{n-1}\xi_k}\tag2.$$ The variance reduction idea is now to replace $\xi_k$ by $$\operatorname E\left[\xi_kf(A_k)\mid A_k\right]=\tilde cf(A_k)w(A_k)\tag3.$$

The problem is that $w(A_k)$ (or $1/\varrho(A_k)$) cannot be evaluated. So, they suggest to estimate it instead.

Xi'an already answered and quoted the estimator which is used in the paper. However, I don't get the derivation of it.

Things should be as follows: At time $\tau_k$ the state $A_k$ was accepted. After $\xi_k$ further steps, $A_{k+1}$ got accepted. That means, we had $$U_{\tau_k+1}>\alpha(Z_{\tau_k+1}),\ldots,U_{\tau_k+\xi_k-1}>\alpha(Z_{\tau_k+\xi_k-1}),U_{\tau_k+\xi_k}\le\alpha(Z_{\tau_k+\xi_k})\tag4.$$ With this in mind, it is obvious that we can write $$\xi_k=1+\sum_{n=1}^{\xi_k-1}\prod_{i=1}^n\underbrace{1_{\left\{\:U_{\tau_k+i}\:>\:\alpha(Z_{\tau_k+i})\:\right\}}}_{=\:1}=1+\sum_{n\in\mathbb N}\prod_{i=1}^n1_{\left\{\:U_{\tau_k+i}\:>\:\alpha(Z_{\tau_k+i})\:\right\}}.\tag5$$ We can note further that $$Z_{\tau_k+i}=(A_k,Y_{\tau_k+i})\tag6$$ for all $i\in\{1,\ldots,\xi_k-1\}$ and hence write $$\xi_k=1+\sum_{n\in\mathbb N}\prod_{i=1}^n1_{\left\{\:U_{\tau_k+i}\:>\:\alpha(A_k,\:Y_{\tau_k+i})\:\right\}}\tag7.$$

In $(7)$ we can clearly replace $\left(U_{\tau_k+i}\right)_{i\in\mathbb N}$ by an arbitrary independent identically $\mathcal U_{[0,\:1)}$-distributed process. The resulting replacement of $\xi_k$ will have the same distribution as $\xi_k$.

However, what I don't get is that they claim that we can also replace $\left(Y_{\tau_k+i}\right)_{i\in\mathbb N}$ by an arbitrary independent process $(H_n)_{n\in\mathbb N}$ (they call it $(Y_n)_{n\in\mathbb N}$, but this symbol is already in use in my post) such that $$\operatorname P\left[H_n\in\;\cdot\;\mid A_k\right]=Q(A_k,\;\cdot\;)\tag8$$ for all $n\in\mathbb N$. Why is the resulting replacement of $\xi_k$ then still distributionally equivalent to $\xi_k$?

Answer 1

Let me try to summarise here the argument of our paper. (Note that I am mostly using different notations, cut&pasted from my set of slides.)

The original Metropolis-Hastings estimate of $\mathbb E[h(X]$ writes as$$ {\delta_n} =\frac{1}{n}\sum_{t=1}^n h(x^{(t)})={\frac{1}{n}\,\sum_{i=1}^{M_n} \mathfrak n_i h(\mathfrak z_i)}\,, $$ after $n$ iterations where

• $(x^{(t)})_{1\le t\le n}$ is the original Markov chain produced by the Metropolis-Hastings kernel (targetting density $\pi(\cdot)$ with proposal density $q(\cdot|\cdot)$)
• $(\mathfrak z_i)_{1\le i\le M_n}$ is the sequence of distinct values in $(x^{(t)})_{1\le t\le n}$
• $(\mathfrak n_i)_{1\le i\le M_n}$ is the corresponding sequence of the number of occurrences of these distinct values in $(x^{(t)})_{1\le t\le n}$

So far, there is no difference in running the MCMC chain, simply an alternative representation of the basic estimate. It is rather straightforward to establish that

The sequence $(\mathfrak z_i,\mathfrak n_i)$ satisfies

1. $(\mathfrak z_i,\mathfrak n_i)$ is a Markov chain
2. $\mathfrak z_{i+1}$ and $\mathfrak n_i$ are independent given $\mathfrak z_i$;
3. $\mathfrak n_i$ [given $\mathfrak z_i$] is distributed as a geometric random variable with probability parameter \begin{equation} {p(\mathfrak z_i) := \int \alpha(\mathfrak z_i,y)\,q(y|\mathfrak z_i)\,\text dy}\,; \end{equation}
4. $(\mathfrak z_i)_i$ is a Markov chain with transition kernel $\tilde Q(\mathfrak z,\text dy)=\tilde q(y|\mathfrak z)\text dy$ and stationary distribution $\tilde \pi$ such that $$ {\tilde q(\cdot|\mathfrak z) \propto \alpha(\mathfrak z,\cdot)\,q(\cdot|\mathfrak z)\quad \mbox{and} \quad \tilde \pi(\cdot) \propto \pi(\cdot)p(\cdot)}\,. $$

Were the density $p(\mathfrak z_i)$ available in closed form, an improved estimator would be $$\delta^* = \frac{1}{n}\,\sum_{i=1}^{M_n} \dfrac{h(\mathfrak z_i) }{p(\mathfrak z_i)}$$ since $$\mathbb E[\mathfrak n_i|\mathfrak z_i]=\frac{1}{p(\mathfrak z_i)}$$ But since the density $p(\mathfrak z_i)$ is not available in closed form, a reduced variance alternative to $\mathfrak n_i$ is available by Rao-Blackwellisation. Using the representation $$ \mathfrak n_i =1+ \sum_{j=1}^\infty \prod_{\ell\le j} \mathbb{I}\left\{ u_\ell \ge \alpha(\mathfrak z_i,y_\ell) \right\}\,, $$ which, again, is not modifying the original Metropolis-Hastings output, we propose an improvement by integrating out the $u_\ell$ uniform (at the potential cost of having to possibly generate more proposed $y_j$'s).

If $(y_j)_j$ is an iid sequence (of proposals) with distribution $q(y|\mathfrak z_i)$, the quantity $$ {\hat\xi_i = 1+\sum_{j=1}^\infty \prod_{\ell\le j} \left\{ 1 - \alpha(\mathfrak z_i,y_\ell) \right\}} $$ is an unbiased estimator of $1/p(\mathfrak z_i)$ which variance, conditional on $\mathfrak z_i$, is lower than the conditional variance of $\mathfrak n_i$, $\{1-p(\mathfrak z_i)\}/p^2(\mathfrak z_i)$.

Therefore, the Rao-Blackwellised estimator $$\hat \delta_n = \frac{1}{n}\,\sum_{i=1}^{M_n} \hat\xi_i h(\mathfrak z_i)$$ is an improvement over $\delta_n$

Answer 2

In order for myself to understand this better, and hopefully for others as well, I have made the graphical representation below.

• bottom panel What the algorithm effectively does is creating a sample according to the distribution below that represents the transitions in the Markov Chain.

• In the middle panel we see the transition probabilities for a specific value between 0.78 and 0.8 and about 11% of the transitions are a repetition of the same value (the high peak in the histogram) and 89% of the transitions are to a new proposed value.

  So for some 'transitions' we actually 'remain in place' and we could alternatively view the process as

  1. Deciding how long to stay in place or how long to repeat the same value. Let's call this a waiting time $t_i$
  2. Deciding which new value to transition to.

• In the upper panel we see that this process generates a sample with a marginal distribution that resembles the density $f(x) \propto max(0,(1-x)^2)$

So given some current value $x_i$ the transitions in the Markov chain are according to a mixture distribution with a part $1-p$ that is a degenerate distribution centered at $x_i$ and a part $p$ that is distribution of the probability density of the next point.

The Metropolis Hastings MCMC will do both of these two processes (1. sample waiting time $t_i$ 2. sample new value) together by using uniform variables $u_i$ and proposal variables $y_i$.

By doing it in that way, by making comparisons with the uniform values, the algorithm will give only discrete values for the waiting time. This is nice when we want a sample (which requires a discrete number of values), but not neccesary for computing a derived value (where we could also use a weight with the average waiting time).

E.g. if a sample has 10% probability to stay in place, and 90% probability to transition to some other value, then we get as waiting times only discrete values of 1, 2, 3, etc. with respective frequencies of 0.9, 0.09, 0.009, etc. whereas we would like to know the average waiting time 0.9+0.092+0.0093+... ≈ 1.11 average.

The idea behind the Rao-Blackwellization is to use a more precise estimate of the waiting time instead of the discrete values.

---

Estimating the average stopping time.

When we consider the average waiting time for making a transition to a different value, then the observed waiting time as estimate might not be the most efficient estimate.

Let's consider a more abstract view of the part 1 of the process where the waiting time is sampled. That part 1, deciding the $n_i$ how long we 'stay in place' can also be considered as being (independently) decided with a follow process:

• We repeatedly draw a sample of transition/acceptance probability values $p_i \sim f(p)$ and a sample $u_i \sim U(0,1)$ of cutoff values.
• We stop the drawing at some time $t$ when $u_t < p_t$
• We choose the observed $t$, the stopping time of the process, as the estimate for the average stopping time.

So we have a sample $y_1, y_2, \dots, y_t$ and $u_1, u_2, \dots, u_t$ and with that sample wish to estimate the average stopping time.

This Metropolis Hastings method uses the observed stopping time, which is dependent of the values $u_i$, is that an efficient estimator of the average 'stay in place' time? Intuitively it seems not. Consider for example the following: we could have large values for $p_i$ (large probability to transition and stop), but by chance we also have very large values $u_i$ which makes us stop very late. Would it be good to estimate a large stopping time? Should we maybe ignore the values of $u_i$ and use the values of $p_i$ only?

... to be continued ... (I have to figure out the sufficient statistic and the distribution of the stopping time given that sufficient statistic)