Hey guys, in this video I want to go over how to do Rao-Blackwellization of a Gibbs sampler. And
so this going to be a continuation of the Gibbs sampling video that I already have posted and
I'll put a link to it in the description box below. And that video was also a continuation
of another video but you don't really have to watch those other videos to know what's going
on in this one. But I am going to do mostly a demo of how to do Rao-Blackwellization after
you've finished doing Gibbs sampling, or when you use Gibbs sampling to estimate a parameter. And
Rao-Blackwellization is basically a variance reduction technique. If you use it for Monte
Carlo sampling, you're probably going to use it in the context of reducing the variance of your
parameter estimate. In the example that I was doing before in the previous video, I had three
parameters that I was trying to estimate. This was from some blood typing example. I wanted
to estimate genetic frequencies of the gene A, gene B, and gene O. I labeled that PA, PB, PO,
and to do that, I said one way of doing that was using Gibbs sampling. The Gibbs sampling method
was basically a three-step Monte Carlo sampling method. I sampled for a latent variable that I
called Z1. This was the estimated count of people with AA alleles. Then I estimated another latent
variable Z2, which was basically an estimate of the people who have BB genetic alleles. And
then I used those two latent variables samplings to plug in to a posterior distribution and the
posterior distribution estimated PA, PB, and PO, using a Dirichlet probability distribution. This
whole process was explained in a previous video. And so basically when you have this three step or
however many step Gibbs sampler method figured out what the actual distributions are and what you're
sampling from, to plug that into R you basically have a sample and then you do Monte Carlo
sampling. You do a couple number of iterations until you reach convergence and then you remove
your burn-in and then you estimate basically from a sampling distribution your parameter estimates.
And then you can also get a standard error for your estimates as well. Rao-Blackwellization is
a method mainly of improving the standard error, or the variance, of your parameter
estimates. These are PA hat, PB hat, and PO hat that we estimate and then we probably also want
to give the standard error for each of those estimates as well. The standard error of PA,
PB, and PO. Let's go ahead and get started with explaining how you could do Rao-Blackwellization
on this particular example. In Gibbs sampling like I just said, we are estimating using a posterior
distribution, which is this distribution here, which was Dirichlet. With Rao-Blackwellization,
the main idea is that you're going to get the marginal distribution--the marginal distribution
for PA, the marginal distribution for PB, and the marginal distribution for PO--and you're going to
use the properties from that marginal distribution to estimate your parameter instead of estimating
directly from the posterior distribution. In this example, the posterior distribution for
these three parameters that were estimating is a three parameter Dirichlet distribution. The
marginal distribution from a three parameter Dirichlet is a beta distribution. In general,
if you have q is distributed Dirichlet with three parameters, let's say alpha1, alpha2, alpha3,
then the formula is that for each parameter qi, the marginal distribution is Beta. And the
beta distribution has two parameters. The first parameter would be alpha i and then the second
parameter is the sum of all the alphas that's not equal to i. The way you could write that would be,
the sum of i prime where i prime is not equal to i of alpha i prime. For this example, that's going
to look like the following three equations. I just pasted the three distributions here. For PA, the
marginal conditioning on all the other parameters and all the sample data, which was in Y and all
the latent variables, basically everything...PA, the marginal of that conditioned on everything,
that would be distributed beta with parameters m1 plus alpha1 and then the second parameter
would be everything else. That would be m2 plus m3 plus alpha2 plus alpha3. That comes from here.
Again if PA is the parameter m1 plus alpha1, that would be the first parameter in the beta
distribution and then everything else would be the second parameter. We have that here. For PB,
it would be the same thing. The first parameter is m2 plus alpha2 and then the second parameter
is the addition of everything else. Same thing for PO. You can calculate always the marginal
distributions for all the parameters that you want to estimate from a posterior distribution,
in this example. Instead of estimating from this distribution, you can use the marginal
distributions. And the reason you're going to want to do this is because once you have these marginal
distributions, you can use the expected value properties from these marginal distributions.
When you have a beta distribution that has two parameters, let's say alpha, beta. The
expected value by definition of a beta distributed random variable is going to be alpha over alpha
plus beta. To do Rao-Blackwellization in this example, we're going to calculate the marginals
first and then use the expected value formula from the marginal distributions--in this case,
its beta distribution--and we're going to estimate this and use this as our parameter estimates. When
we do this, it's going to basically give us the same PA hat, PB hat, and PO hat that we had before
but the variances of those parameter estimates will be reduced. If we use this expected value
formula for this example, we would get these three equations here and I'm just going to focus on the
first part of the equations for the PA, PB, and PO. The parameter estimates would be again alpha
over alpha plus beta. If we go back to the first marginal distribution, that would be m1 plus
alpha1 over m1 plus m2 plus m3 plus alpha1 plus alpha2 plus alpha3. From the blood-typing example
that we've been doing, we have been setting alpha1 is equal to alpha2 is equal to alpha3 as equal to
1. In this case, we have alpha1 plus alpha2 plus alpha3 is equal to 1 plus 1 plus 1, which would
give you 3. This was explained in a previous video and basically it's the hyperparameters
of the Dirichlet distribution. If you set those alpha parameters equal to 1, you're assuming no
information or a noninformative prior. That's how we get the formulas using Rao-Blackwellization.
Now, these m values are parameters that we estimate using the Monte Carlo sampling method.
We are still using Gibbs sampling and we're still using this three-step Gibbs sampling to estimate
these parameter estimates and we do that to plug into these m values. That's what the second half
of these equations is showing because we don't actually know the true values of m1, m2, and m3.
We estimate that given our sample using the Gibbs sampling method. If we did that, we basically have
m, which is a number of iterations in our Gibbs sampler, and then we just plug in--this is just
expanding out what m1 was equal to, what m2 was equal to, and what m3 was equal to--and getting
a sample estimate based on a certain number of iterations in the Monte Carlo sampling method.
This is kind of like a quick explanation of how you would do Rao-Blackwellization. I think will
be a lot more interesting to show the actual implementation of this in R. Let me go ahead
and move over to R. Okay so I'm moving on to the R demo and I will put the R code in the
description box below. The first half of the R code just runs the Gibbs sampler. I already did
a video on all of this code and I will put a link again to the video in the description box
below. But the first part of the code just does the Gibbs sampling. The Gibbs sampler was based
on 2,250 observations and I'm going to remove the first half as burn-in. I'm just going to use the
final 1,125 iterations to actually estimate my parameters. Then when I estimate PA, PB, and PO,
I get these values. This is again just using the Gibb sampling method. Those are the PA hat,
PB hat, PO hat. And then from that Gibbs sampling distribution that I created, if I just take the
standard deviation of that Gibbs sampler, then I get the standard errors using the traditional
Gibbs sampling method. The standard errors of my three parameter estimates are .008, .005,
and .009. It's really not that bad but I'm just going to do Rao-Blackwellization in this example
anyway and we're going to see how that's going to reduce the standard errors of these parameter
estimates. For Rao-Blackwellization, I'm going to show the density plots of what the parameters look
like just doing the Gibbs sampling method and then I'm going to overlay what it looks like if you do
Rao-Blackwellization. For PA if you just do the Gibbs sampling method, the sampler will produce
a sampling distribution that's shown in the black line. And then the red dotted line is what
the sampling distribution looks like when you do Rao-Blackwellization. You can see that the Gibbs
sampler was not bad but it's not completely...it's kind of unimodal, the mean is probably somewhere
around here, as you can see. But it has kind of like a spread out peak whereas the peak for the
Rao-Blackwellization method is a lot pointier and obviously the variance or standard error is
going to be greatly decreased. That's a plot. And then if you actually plug in to get the parameter
estimates, the mean would be .2658. That's not that different from using the Gibb sampling
method. It's basically the same thing. But then if we plug in for the standard deviation of that
sampling distribution, we're going to get the standard error. And so for the Rao-Blackwellized
parameter estimate, the standard error is .004, whereas before the standard error was .008. It
pretty much cuts in half the standard error for PA, that was one of the parameter estimates. You
can do the same thing for PB. Just plotting the different sampling distributions for the Gibbs
sampler in black and the Rao-Blackwellized sample in red. Again you see that it's going to reduce
the variance substantially and also probably give you more confidence in the parameter
estimate itself because it's going to be much more peakier and unimodal. And the actual values are
as follows: PB hat would be .09 and the standard error is .003. We can compare .003 to the previous
standard error of .005 and we see that it's been reduced. Finally for the third parameter, PO. We
see the same thing. This sampling distribution was not super peaky. It was kind of spread out again
like before and Rao-Blackwellization makes the sampling distribution a lot less spread out. And
we can see that for the standard error, the new standard error using Rao-Blackwellization is .006,
whereas before it was .009. Hopefully this short demo showed that Rao-Blackwellization can be
applied to a Gibbs sampler and it's a really good idea to apply it to the end of Gibbs sampling
technique because it's going to improve your parameter estimates and especially it's going to
reduce the variance of your parameter estimates in particular. And that's it for this video.
*chiptune music*