# Rao-Blackwellization of Gibbs Sampler

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.

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.

• So assuming the posterior distribution of the parameter is known (which to the best of my knowledge happens to be true when applying Gibbs sampling), taking the mean of the distribution rather than a random sample would be the Rao-Blackwellized estimator? I hope I understood your answer correctly. Thank you very much already! – mscnvrsy Jul 18 '16 at 21:43
• That is incorrect. In Gibbs sampling, you don't know the posterior distribution of the parameter, but know the full conditional posterior for each parameter. There is a big difference between the two. Above, the posterior is $f(\mu, \phi \mid data)$ which is unknown, and for the Gibbs sampler to work you need to know both $f(\mu \mid \phi, data)$ and $f(\phi \mid \mu, data)$. And you are also incorrect in your second understanding. You still need to take a sample from the marginal posterior of $\mu$, and then calculate the sample mean of $\phi^*$ using those samples to find the R-B estimator. – Greenparker Jul 18 '16 at 21:59
• @mscnvrsy I added an example to help – Greenparker Jul 18 '16 at 22:12
• Wow, thank you very much for clarifying this to me. So assuming that I know the full conditional distributions, I can work with the theoretical means of the conditional distributions and average over these theoretical means (such as E[phi|mu,y]) to obtain the RB-estimator? This would then minimize the variance of my parameter estimates? – mscnvrsy Jul 18 '16 at 22:16
• If you were obtaining independent samples, yes it would minimize the variance of estimators, however, since you are dealing with Markov chains, its generally known that RB does not necessarily reduce variance, and there are some instances where the variance even increases. This paper by Charlie Geyer gave some examples to this point. – Greenparker Jul 18 '16 at 22:21

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): 