Rao-Blackwellization of Gibbs Sampler

Question

I am currently estimating a stochastic volatility model with Markov Chain Monte Carlo methods. Thereby, I am implementing Gibbs and Metropolis sampling methods.

Assuming I take the mean of the posterior distribution rather than a random sample from it, is this what is commonly referred to as Rao-Blackwellization?

Overall, this would result in taking the mean over the means of the posterior distributions as parameter estimate.

Answer 1

Assuming I take the mean of the posterior distribution rather than a random sample from it, is this what is commonly referred to as Rao-Blackwellization?

I am not very familiar with stochastic volatility models, but I do know that in most settings, the reason we choose Gibbs or M-H algorithms to draw from the posterior, is because we don't know the posterior. Often we want to estimate the posterior mean, and since we don't know the posterior mean, we draw samples from the posterior and estimate it using the sample mean. So, I am not sure how you will be able to take the mean from the posterior distribution.

Instead the Rao-Blackwellized estimator depends o the knowledge of the mean of the full conditional; but even then sampling is still required. I explain more below.

Suppose the posterior distribution is defined on two variables, $\theta = (\mu, \phi$), such that you want to estimate the posterior mean: $E[\theta \mid \text{data}]$. Now, if a Gibbs sampler was available you could run that or run a M-H algorithm to sample from the posterior.

If you can run a Gibbs sampler, then you know $f(\phi \mid \mu, data)$ in closed form and you know the mean of this distribution. Let that mean be $\phi^*$. Note that $\phi^*$ is a function of $\mu$ and the data.

This also means that you can integrate out $\phi$ from the posterior, so the marginal posterior of $\mu$ is $f(\mu \mid data)$ (this is not known completely, but known upto a constant). You now want to now run a Markov chain such that $f(\mu \mid data)$ is the invariant distribution, and you obtain samples from this marginal posterior. The question is

How can you now estimate the posterior mean of $\phi$ using only these samples from the marginal posterior of $\mu$?

This is done via Rao-Blackwellization.

\begin{align*} E[\phi \mid data]& = \int \phi \; f(\mu, \phi \mid data) d\mu \, d\phi\\ & = \int \phi \; f(\phi \mid \mu, data) f(\mu \mid data) d\mu \, d\phi\\ & = \int \phi^* f(\mu \mid data) d\mu. \end{align*}

Thus suppose we have obtained samples $X_1, X_2, \dots X_N$ from the marginal posterior of $\mu$. Then $$ \hat{\phi} = \dfrac{1}{N} \sum_{i=1}^{N} \phi^*(X_i), $$

is called the Rao-Blackwellized estimator for $\phi$. The same can be done by simulating from the joint marginals as well.

Example (Purely for demonstration).

Suppose you have a joint unknown posterior for $\theta = (\mu, \phi)$ from which you want to sample. Your data is some $y$, and you have the following full conditionals $$\mu \mid \phi, y \sim N(\phi^2 + 2y, y^2) $$ $$\phi \mid \mu, y \sim Gamma(2\mu + y, y + 1) $$

You run the Gibbs sampler using these conditionals, and obtained samples from the joint posterior $f(\mu, \phi \mid y)$. Let these samples be $(\mu_1, \phi_1), (\mu_2, \phi_2), \dots, (\mu_N, \phi_N)$. You can find the sample mean of the $\phi$s, and that would be the usual Monte Carlo estimator for the posterior mean for $\phi$..

Or, note that by the properties of the Gamma distribution $$E[\phi | \mu, y] = \dfrac{2 \mu + y}{y + 1} = \phi^*.$$

Here $y$ is the data given to you and is thus known. The Rao Blackwellized estimator would then be

$$\hat{\phi} = \dfrac{1}{N} \sum_{i=1}^{N} \dfrac{2 \mu_i + y}{y + 1}. $$

Notice how the estimator for the posterior mean of $\phi$ does not even use the $\phi$ samples, and only uses the $\mu$ samples. In any case, as you can see you are still using the samples you obtained from a Markov chain. This is not a deterministic process.

Answer 2

The Gibbs sampler can then be used to improve efficiency of (say) samples from a marginal posterior, call it $\pi_2(\theta_2|y)$. Note \begin{eqnarray*} \pi_2(\theta_2|y)&=&\int \pi(\theta_1,\theta_2|y)d\theta_1\\ &=&\int \pi_{2|1}(\theta_2|\theta_1,y)\pi_1(\theta_1|y)d\theta_1\\ &=&E(\pi_{2|1}(\theta_2|\theta_1,y)) \end{eqnarray*} Thus, the marginal density of $\theta_2$ at some value $\theta_2$ is the expected value of the conditional density of $\theta_2$ given $\theta_1$ at the point $\theta_2$.

This is interesting due to the Variance Decomposition Lemma $$ Var(X)=E[Var(X|Y)]+Var[E(X|Y)], $$ where the conditional variance $Var(X|Y)$ is $E\left\{(X-E(X|Y))^2|Y\right\}$. Also, $Var(E(X|Y))=E\left[(E(X|Y)-E(X))^2\right]$. In particular, $$ Var(X)\geq Var[E(X|Y)]. $$ A Gibbs sampler will give us realizations $(\theta_{1i},\theta_{2i})$. The upshot is that it is better to estimate $\pi_2(\theta_2|y)$ by $$ \hat{\pi}_2(\theta_2|y)=\frac{1}{M}\sum_{i=1}^M\pi_{2|1}(\theta_2|\theta_{1i},y) $$ than by some conventional kernel density estimate using the $\theta_{2i}$ for the point $\theta_2$ - provided we know the conditional distributions (which is of course why we use Gibbs sampling in the first place).

Example

Suppose $X$ and $Y$ are bivariate normal with means zero, variances 1 and correlation $\rho$. That is, $$ \pi(x,y)\propto\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2+y^2-2\rho x y)\right\} $$ Clearly, marginally, $Y\sim N(0,1)$, but let us pretend we do not know this. It is well-known that the conditional distribution of $Y$ given $X=x$ is $N(\rho x,1-\rho^2)$.

Given some $M$ realizations of $(X,Y)$ the "Rao-Blackwell" estimate of the density of $Y$ at $y$ then is $$ \hat\pi_Y(y)=\frac{1}{M}\sum_{i=1}^M\frac{1}{\sqrt{1-\rho^2}\sqrt{2\pi}}\exp\left\{-\frac{1}{2(1-\rho^2)}(y-\rho x_i)^2\right\} $$ As an illustration, let us compare a kernel density estimate to the RB-approach

library(mvtnorm)

rho <- 0.5
R <- 50
xy <- rmvnorm(n=R, mean=c(0,0), sigma= matrix(c(1,rho,rho,1), ncol=2))
x <- xy[,1]
y <- xy[,2]

kernel_density <- density(y, kernel = "gaussian")
plot(kernel_density,col = "blue",lty=2,main="Rao-Blackwell estimates from conditional normals",ylim=c(0,0.4))
legend(1.5,.37,c("Kernel","N(0,1)","Rao-Blackwell"),lty=c(2,1,3),col=c("blue","black","red"))
g <- seq(-3.5,3.5,length=100)
lines(g,dnorm(g),lty=1) # here's what we pretend not to know

density_RB <- rep(0,100)
for(i in 1:100) {density_RB[i] <- mean(dnorm(g[i], rho*x, sd = sqrt(1-rho^2)))}
lines(g,density_RB,col = "red",lty=3) 

We observe that the RB estimate does much better (as it exploits the conditional information):