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
- $(U_1,\ldots,U_n)$ and $(Z_1,\ldots,Z_n)$ are independent;
- $$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$?