Hey guys, in this video I want to go over how to do rejection sampling which is
also known as the accept-reject method or accept-reject algorithm. This is
one type of Monte Carlo random sampling method and for this video I want to
first introduce the algorithm and the concept behind it and then do an example
in R and by hand afterwards. So let's go ahead and get started. The idea, or the
motivation behind rejection sampling, is you have some target distribution pi of
X and this is some distribution that you want to get random samples from but it's
going to be difficult to sample from pi of X directly. So maybe it's a complex
function and I'm going to just draw some random PDF and let's call that our
target distribution pi of X and we want to be able to sample from this target
distribution. The way the rejection sampling method works is instead of
sampling directly from pi of X, we're going to sample from some easy to sample
from distribution. So some well-defined or well-known distribution like the
normal distribution, the uniform distribution, or some distribution that's
easy to sample from. Then the rejection sampling method says as long as there's
some g of X distribution--proposal distribution--that's easy to sample from
that can be multiplied by a constant c, if c times g of X is always greater than
or equal to the target pi of X, then you can use a random sample from c times g
of X to get a random sample from pi of X and I'll go over how the exact algorithm
works in a second. But basically the idea behind this is going to be that you kind
of want to create an envelope over your target distribution pi of X. So to
describe how this rejection sampling method works visually, let's say this is
our PDF pi of X. You basically want your proposal distribution times a constant
to envelop your target distribution. By that I mean, let's say our proposal
distribution is the uniform distribution. So the uniform distribution, when you
graph it as a PDF is just a straight line. Then if you do this straight line,
whatever uniform distribution this is, times the constant c--
let's say we multiply this by a positive constant c, so this is c times g of X--
then for all points on the support of X of the pi of X distribution, the c times
g of X distribution envelops the target distribution and the reason you want
this envelope to be created for all the support of X is that then we're going to
sample and we're either going to accept or reject the new samples according to
where they lie in this envelope. So if they lie under the target distribution
and under c times the proposal distribution, then we're going to accept
the samples if they lie above the target distribution but below the proposal
distribution we're going to reject those samples so by doing that we have samples
that we can easily sample from and then we're just going to accept a subset of
those samples and then the subset of those samples will be distributed
according to our target distribution. Formally the algorithm is we're going to
draw samples from our easy to sample from distribution g of X. Then we're
going to draw random samples from the uniform [0, 1] distribution. Then we're
going to accept each random sample X of i according to the following formula:
we're going to accept X of i if U of i, which is each random sample from the
uniform distribution, is less than or equal to pi of X of i divided by c times
g of X of i and we're going to check this for all i, otherwise we reject.
From our sample of X, we're going to accept some of these X of i's and
reject some of them depending on whether or not they fall under the target
distribution and we can check that using this formula. Then our accepted random
sample of X of i's will follow the target distribution pi of X that we were
looking for in the first place. So that's the algorithm and the concept behind the
algorithm. Now let me go ahead and do an example. Okay so I'm going to go ahead
and go over this example. It says to design a rejection sampling algorithm
for generating random samples from this target distribution pi of X and when we
look at pi of X, it's a distribution that's difficult to sample from directly.
It's basically a third degree polynomial with support from zero to one and so
this is not a well-known distribution and it would be difficult to get a
random sample from here directly. So we're going to use the rejection
sampling method to sample from this distribution. I went ahead and I plotted
in R the PDF of our target distribution pi of X. This is pi of X
on the support from 0 to 1. This is the PDF that we want to randomly sample
from. So let's go over how we would do one way of doing the rejection sampling
algorithm for this target distribution. This distribution is continuous with
finite support from 0 to 1. This PDF only has to exist on the support from 0
to 1 and it's a continuous distribution. So this has a special case where the g
of X, our proposal distribution, can always come from the uniform
distribution. When you're doing rejection sampling and you want to know
what should I use as the proposal distribution, there's never one right
answer. It could be whatever comes to mind first or whatever works for you but
when you have a continuous distribution over finite support, the uniform
distribution can always be your proposal distribution. If we think of the
uniform distribution as we said before, if we plot that as a PDF it just looks
like a straight line. Since this support is from 0 to 1, a uniform 0 to 1
distribution would just look like this. This is uniform 0 to 1. The PDF would be
the value 1 at Y is equal to 1 from 0 to 1 and it would just be a straight line.
We have that and that works over our support from 0 to 1 but now we need to
make sure to multiply our proposal distribution, the uniform 0 to 1 random
variable, by a constant c to make sure that our proposal distribution envelops
the target distribution. In other words, our target distribution has to be
at least up to here so that all the values that we sample from we can make
sure are beneath our target distribution. If I look at what is the value of the
target distribution at the end of the support, which is at X is equal to 1 in
other words if I plug in 1 to the PDF, that will give me a
value of 3.125 and since this is basically a monotonically increasing
function, basically we see that the highest value or the largest value of pi
of X will only ever be 3.125. So if this is 1 and I want to make sure that I
get at least the highest value of the target distribution pi of X, then we can
say or we know that the constant that we should multiply our proposal
distribution by should be 3.125 or you can think of it as
it should be at least as big as 3.125. You can make it
bigger but that will just mean you end up rejecting more samples, so if you want
the tightest uniform distribution that will have you reject the fewest number
of samples it should be 3.125 and that's also the easiest
value to come up with. Now we have our proposal distribution or sometimes it's
called the instrumental distribution, which can be g of X is uniform [0,1]
and we know that if we multiply that by c, which is 3.125, then
we will always be under the target distribution and we can do rejection
sampling. If I wrote that out in words instead of trying to do it visually I
can say we're going to draw random samples X from the proposal or
instrumental distribution uniform [0,1]. Then we're going to draw uniformly
random samples from uniform [0,1] again. And then we're going to accept X
of i if the uniform random samples are less than or equal to pi X of i divided
by 3.125 times g of X of i and otherwise we'll reject.
Then we expect our accepted X random samples will be distributed according to the
target pi of X. Okay so that's our algorithm for this example, so now let's
go ahead and do this in R. So I switched over to R and I'm going to go over this
short R demo. I will try to put the code for everything in the description
box below but it's a pretty simple algorithm and it's pretty simple to code
it yourself.The first thing I did was just include the code for how I printed
out the PDF. Again by hand I usually do this with a graphing calculator but this
is how you would do it in R. You basically just passed
through the function of your target distribution pi of X to the function called
curve over the support 0 to 1, which is our support of X, and then I plotted the
ab line, which is just a straight line which would be our proposal distribution
g of X times c. So our proposal distribution would be--basically at y is
equal to 1 it would be a straight line from 0 to 1. That's the uniform [0, 1]
proposal distribution and then if we multiply that by the constant 3.125 it
will give us this line here and that covers all the support of the target
distribution pi of X. And the next part is the actual rejection sampling code or
the rejection sampling algorithm. I randomly sampled from our proposal
distribution g of X which is uniform [0, 1]. I'm calling those random samples X
and I took 4,500 random samples from the uniform distribution. We're going to
reject some of those samples so our accepted final random sample number
will be smaller than 4,500 but I just picked that number to make sure it was
more than a thousand. In this example, I'm trying to at least get a thousand
random samples that are distributed according to our target distribution.
The second step is we have to sample from the uniform [0, 1]. In my example just by
chance, our proposal distribution was also uniformly distributed from 0 to 1
but that doesn't always have to be the case. In this case, it just happens to be
the same. Then this line of code runs a function or saves a function as pi of X.
We're going to use this as part of our accept-reject formula, so I need to
have this function here. Then I set up a counter starting at 1 and accept is
going to be an empty vector that I fill with the accepted random samples from X.
So I'm going to run that and then this is how I wrote the accept-reject
algorithm. Basically I said that while we still have simple random samples, so
while our sample that we've been cycling through is less than 4,500 and while we
have fewer than 1,000 samples that are accepted
keep cycling through this loop. So basically I keep cycling until you run
out of samples or our sample size reaches a thousand which is what we want
and then within the while loop we basically want to make sure that we
accept X of i if each U of i is less than or equal to pi of X of i over c
times g of X of i. So our U of i is going to be test_u and that's just
going to be--we're going to pick one random sample from U and then our X of
i is going to be pi of X of i divided by 3.125, which is our c
times our g of X. Our g of X again is just from the uniform distribution. Then
I have an if loop, so if the U of i is less than or equal to the test
underscore x which is our fraction, we're going to add it to our accepted vector.
If it's not we're just going to increase the counter. So if we accept it, we accept
it and we add it to our vector of accepted random samples and then we
increase the counter, otherwise we just increase the counter meaning we reject
the sample. So I'm going to run that and the value of the counter will tell me
how many values it took for me to get 1,000 random samples. So from 4,500
random samples that I thought I would need I ended up only needing
4,188 and our accept vector now has a thousand
values and I'm going to plot the histogram of the accepted vector,
the vector of all the accepted random values, to show that it should look similar to
the PDF. If I run this last line of code this will give me the histogram of
the values in our accept vector, which we're saying is a random sample from the
target distribution and if I compare the histogram to the continuous PDF, we see
that empirically it does look like it is coming from this PDF. It looks like
the values are correctly being sampled from our target distribution. So it looks
like the rejection sampling method does work and that's it for this video and
thank you guys for watching!
*chiptune music*