Hey guys, in this video I'm going to go over an example of doing Gibbs sampling,
which is a type of Monte Carlo method and I'm going to do Gibbs sampling with
data augmentation. The example I'm going to do in the video is based on an
estimation problem that I'm going to try to solve using Monte Carlo methods.
I actually talked about the same problem in the last video I did because this is
like a project I was working on, so I will link in the description box below
the video if you want more background information on this problem but I think
you can just watch this video and understand the main idea of the problem
without watching the other video. Basically the example I'm applying Gibbs
sampling methods to is a blood typing problem, so we have sampled data from a
population of people's blood type and you have one out of four possible blood
types that you can be. You can have blood type A, B, AB or O and these are the
sample data below. And that has a probability as displayed in this column
here and your observed blood type depends on one of six genotypes.
You can have genotype AA, which would give you blood type A but AO would also
give you a blood type A as well. So that distribution of your genotype--underlying
genotype--also has a multinomial distribution with the following
underlying probabilities. We usually do not observe genotype, we usually just
observe your phenotype when we do blood testing. We observe this and we kind of
want to estimate the probability of having a certain genotype and what we
want to estimate in particular is the frequency--the genetic frequency--of the
different alleles. So the A, O, and B. The problem is we have this sample data and
we want to estimate pA, so let's say the proportion of A, the proportion of B, and
the proportion of O. That's the main idea with the problem and in this video
we're going to use Gibbs sampling with data augmentation to solve this. I'm
going to in this video estimate these parameters using Gibbs sampling to
estimate from a posterior distribution. So this idea comes from Bayesian
statistics, so this is also kind of like a Bayesian statistics video and the idea
for doing this type of analysis comes from this paper "The Calculation of the
Posterior Distributions by Data Augmentation" by
Tanner and Wong. I will also put a link in the description box below for
this paper and anyway let's go ahead and get started with the video. This is gonna
be kind of like a demonstration more than an in-depth tutorial. I'm kind of
assuming you know some background information on all of this, otherwise the
video would be like 10 hours long. So we want to estimate pA, and pB, and pO where
the sum of the alleles has to equal one. So pA plus pB plus pO has to equal one
and I said we're basically going to use Bayesian statistics to estimate those
three parameters. So we need to calculate or sample from the posterior
distribution. For now, let's assume that we could observe the genetic data, which
we don't observe, we observe the phenotypic data but we're gonna get to
that in a minute but let's say we know the underlining genotypes of our samples.
Then a posterior distribution could be defined as follows I'm going to say P of
pA, pB, pO given the genetic data is proportional to the probability of your
genetic data given the parameters pA, pB, pO times a prior distribution on the
parameters. In bayesian statistics, this is one way that you would define the
posterior distribution, so it's a distribution of the parameters given
your data and that's proportional to a likelihood, which is the probability of
your data given your parameters times a prior distribution on the parameters.
We would say the posterior is proportional to a likelihood times a
prior. The problem with this example is that we do not actually observe the
genetic information. Instead of observing the genetic data, we observe phenotypic
data. So instead of calling this X, I'm going to call this Y because we observe
one of four possibilities instead of one of six possibilities from the genetic
data. So we cannot actually sample from the posterior distribution yet.
Instead of observing X, we observe Y, which I'm going to define as y1 y2 y3 y4
and these are the observed sample data of the phenotypes. So I'm going to say y1
is 89, y2 is 642, y3 is 195, and y4 is 657.
If we considered only the phenotypic data that has four possible outcomes and
we tried to draw from a posterior distribution given Y instead of given X,
we would not have a complete likelihood. We'd have a problem here because we
would not be able to get a conjugate prior to the posterior distribution.
We would not be able to assign a conjugate prior and that would kind of
complicate our ability to estimate pA, pB, pO. Because we prefer having a complete
likelihood to be able to get a conjugate prior for the distribution, this is where
a data augmentation is going to help us. Data augmentation, as the name implies,
basically means we're going to try to add to our sample data that we have.
We're going to try to augment it to be able to form a complete likelihood. So if
we go back to the data that we do observe. We do observe AB, we do observe
O but for A we don't know if someone is AA or AO and for B we don't know if
someone is BB or BO. So in that case we have like two missing variables and so
with data augmentation what you do is you create latent variables for your
missing data--quote-unquote missing data. I'm going to define Z1 as my first
latent variable and I'm going to estimate that and that's going to be an
estimate of the number of people who have AA genotype and I'm going to do
the same thing for the observed phenotype B. I'm going to say Z2 is
going to be an estimate of the people who are BB. Using these latent variables
and estimating these values, I can then plug it in to my posterior distribution
and basically augment the posterior distribution, augment the likelihood, and
then I will be able to form a conjugate prior. So if we have the probability of
pA, pB, pO given Y and given Z then that will be proportional to the likelihood
of Y and Z given the parameters times a prior. Instead of trying to write all
that out, I just pasted it here and so I defined theta as a vector of my
parameters. So pB, pA, and pO are in the vector theta. So we have the posterior
distribution of our parameters given observed phenotypes plus the
latent variables is proportional to the likelihood, which is now a complete
likelihood, times some prior distribution on the parameters and then we can do
multiplication and simplification to basically get that we have that the
posterior distribution is proportional to pA to m1, pB to m2, pO to m3, where m1, m2, and
m3 are defined as follows. That's a mix of the sample data and the latent
variables which we have to estimate. And now because we have a complete
likelihood, we can form a conjugate prior to this distribution. This is a
complete likelihood for the multinomial distribution for the genetic data using
the latent variables and from that we can assume or we can place a conjugate
prior on P of theta. So we're assuming the prior is going to be Dirichlet
distributed because that's the conjugate prior for a multinomial distribution.
I'm going to place a Dirichlet prior on the parameters. That means that the
probability of pA, pB, pO is going to be proportional to pA to the alpha1 minus
1, pB to the alpha2 minus 1, and pO to the alpha3 minus 1, where alpha1, alpha2,
and alpha3 are Dirichlet parameters and in Bayesian statistics because this is a
prior distribution you would call these three parameters hyperparameters.
By substituting this prior distribution back into the equation for the posterior
distribution, you would then get that the posterior distribution--and this makes
sense because it's a conjugate prior--the posterior distribution would also be
distributed as a Dirichlet distribution. So now we're going to have that the
probability of theta given Y, Z is Dirichlet distributed with parameters
m1 plus alpha1, m2 plus alpha2, and m3 plus alpha3. In Bayesian statistics when
you don't know anything about alpha1, alpha2, and alpha3, which are the hyper-
parameters, you can use what's called a noninformative or flat prior. For the
Dirichlet distribution, the noninformative prior would just be setting
alpha1, alpha2, and alpha3 equal to 1. I'm going to use that in this problem.
I'm going to assume I don't know anything about the prior distribution, so
I'm going to set that to 1. Our posterior distribution is going to be
distributed Dirichlet with m1 plus 1, m2 plus 1, and m3 plus 1. Then now that we
have defined the posterior distribution we can actually sample from this using
Gibbs sampling. So this is actually the first part where Monte Carlo sampling is
going to take place. We have to use Gibbs sampling because we don't actually know
the true value of the latent variables so we have to estimate those. In Gibbs
sampling, you can estimate those by sampling from a conditional distribution.
Basically the way Gibbs sampling works is it's an iterative Monte Carlo
method and at each iteration we're going to sample from a conditional
distribution. We're going to condition on the parameters to estimate Z1.
Then using a different conditional distribution we're going to estimate Z2.
We're going to use those estimates to sample from the posterior distribution.
Then by sampling from the posterior distribution it's going to update the
values that you would use to sample from the conditional distributions for Z1 and
Z2. And basically you just iteratively update these values until they converge
to your final estimates. And then you don't really need the estimates for Z1
and Z2 anymore, you just use the final estimates of your posterior distribution
for whatever you want. In this case, we want to estimate the three parameter
values in the posterior distribution. So I have some R code that I'm going to
show from when I did that, but to show you what the distributions that I'm
sampling are: they're going to be Z1 conditional on Y and then the parameters,
pA, pB, pO. That conditional distribution is going to be binomial and it's going to
be binomial y2 and then pA squared over pA squared plus 2pA
pO. I'm not going to show how that was derived but you can derive that just
using math. And then the second conditional distribution is going to be
y given the parameters and that second conditional distribution is also going
to be distributed binomial with parameters y3 and then pB squared over
pB squared plus 2pBpO. So then when we actually run the Gibbs sampler, at each
iteration we're estimating a value of Z1, we're estimating a value of Z2, then
we're estimating the pA, pB, and pO values. So we sample from the conditional,
sample from a conditional, then sample from the posterior distribution. At the
next iteration, we plug back in the estimates of pA, pB, pO into the
parameters of the conditional and then the posterior distribution also uses
these estimates. So each iteration has to update all the values and then you do
that until the sampler converges. In the last part of the video, let me actually
go over the Gibbs sampling code that I have and I will put a link in the
description box below for this code but let's just go over it briefly in the
video as well. Because we need the Dirichlet distribution to sample our
posterior distribution, I defined a function that would do that sampling.
I think probably R does have packages or does have an ability to sample from the
Dirichlet distribution but I couldn't find a function that worked exactly how
I wanted it to. So since the Dirichlet is related to the gamma distribution you
can actually just define it manually and this worked for me. So this is a Dirichlet
with three parameters a, b, c, that you plug in. Then for the next step I had to
predefine some vectors and matrices for doing the Monte Carlo sampling. I kind
of just randomly picked that I was going to do 2,250 iterations of the Monte
Carlo sampling. Because we have to estimate the latent variables and then
the three parameters, I initialize an empty vector for Z1 for all the
iterations; for Z2 for all the iterations. So these are just two zero vectors.
And then I created a zero matrix with all the iterations and three columns for the
three parameter values I'm trying to estimate, pA, pB, and pO. I have my starter
values, which are the sample data. With Monte Carlo sampling you always have to
start with certain initial values, so for Z1 I just set the initial value to 500.
The same thing for Z2 and for my parameter estimates I set the initial values to
one-third for each value and the hyper parameters I'm using are 1. So alpha1 is
equal to 1, alpha2 is equal to 1, and alpha 3 is equal to 1 as well. So then
the next line of code actually runs the Gibbs sampler for the 2,250 iterations.
I conditionally estimate from the binomial distribution for Z1. I estimate
the conditional binomial distribution for Z2. Then I plug those values into the
parameters m1, m2, and m3 and then plug those values into the posterior
distribution at each iteration. So each iteration is going to update the values
of Z1, Z2, and pA, pB, and pO. In Bayesian statistics, once you have sampled from
the posterior distribution one thing you can do is get the posterior means and so
that could be one possible way of estimating the parameter values. You can
also look into estimating the mode. That's sometimes used in Bayesian
statistics but in my particular case I estimated the values using posterior
means. Then that just means I had to find the mean of the sample values from
the Monte Carlo sample. You usually do burn in with Monte Carlo methods so I
just threw away the first half of the samples and I took all the values of the
second half. So I took basically 1,125 values in the second half of the vectors
or matrices. So that gives me the estimates 0.266, that's one parameter
value; 0.094 as a second parameter estimate; and 0.64 for another
parameter estimate. And just to show that the values did converge at the end I
plotted the density plots, just to show that the density plots pretty much show
that there is one unimodal mean for the distributions. So for this problem, after
doing all that work, we have that pO hat is equal to 0.640; pA hat was equal
to .266; and pB hat was equal to .094.
If you watched my previous video, the true values are .638,
.268, and .095. So this gives me values that are
very close and tells me that the Gibbs sampling method worked and these are the
estimates for the underlying genetic frequencies from the population.
And that's it for this video and thank you guys for watching~~
*chiptune music*