Hey guys, in this video I wanted to go over important sampling with an example
and some R code. And importance sampling is a type of Monte Carlo method.
You usually use it to estimate some moment of a distribution. So you have a
random variable X and it has some distribution and usually you want to
estimate some property of the distribution. So usually that's going to be
you want to estimate maybe the mean or you want to know the variance of X,
so some property from a target distribution and this is one Monte Carlo method of
getting those estimators. Important sampling, which is sometimes abbreviated
IS, is a simple extension of the idea of naive Monte Carlo.
The idea with naive Monte Carlo is X your random variable is distributed pi X,
where pi X is some target distribution, and if you wanted to estimate for
example the mean of either X or a function of X--so let's say the expected
value of h of X. H of X is just a function of the random variable X--if you
wanted to find the expected value, which is the mean, with naive Monte Carlo we
just know that this would be equal to the integral of h of X times pi of X dx.
If this is true, you can just estimate mu by getting the sample mean. You can
do resampling using Monte Carlo methods and then get the sample mean of a sample
that you created and that would give you an estimate of your mean and you can
also find the standard error or an estimate of the variance as well.
That's naive Monte Carlo and that's one way of doing Monte Carlo methods to get
information from a target distribution. Importance sampling is a simple
extension of this equation. The main idea is instead of sampling directly
from our target distribution pi of X, we're going to have a proposal
distribution g of X. When we sample X the random variable X is going to be
sampled from a proposal distribution g of X instead of the target. And then if
we want to know the mean still from our target distribution we can say that mu
is going to be equal to the integral of h of X times pi of X dx and we can say that
that's equal to the integral of h of X times pi of X times g of X divided by g
of X dx. Basically we're just multiplying and dividing g of X and when we do that
this is going to basically give us the expected value with respect to the
proposal distribution g of h of X pi of X over g of X. The reason we do this
extra step in important sampling of multiplying by g of X and dividing by g
of X is because now we can use a sample from a different distribution to
estimate the mean from the target distribution. And this is basically
important sampling in a nutshell. So this formula, if you know this then you know
enough to be able to write an important sampling algorithm to estimate certain
properties of a target distribution. And the reason you would do important
sampling, the reason you would do this step of creating a proposal distribution
and sampling from that instead of your target distribution, could be a couple of
reasons. One reason could be that the target distribution is hard to sample
from directly and another common reason of doing important sampling Monte Carlo
methods is variance reduction. You're going to see in the example that I'm
going to do in a minute that when you do important sampling instead of naive
Monte Carlo, you're going to get an estimator that has smaller standard
errors and also sometimes it will be more accurate as well as the naive Monte
Carlo method. So it's often used as a variance reduction technique as well.
In this video, I'm going to do this really simple example of important
sampling as a variance reduction technique. This is an example where you
could use naive Monte Carlo but I'm going to compare in R the naive Monte
Carlo algorithm and the importance sampling algorithm and show that the
importance sampling algorithm reduces the variance, so it's basically a better
way of estimating the, in this case, expected value of the random variable.
This example says suppose X is exponentially distributed with rate
lambda is equal to one and you want to estimate
the expected value of a function of X. The function is e to the negative X plus
cosine of X. We're going to do that two ways in the R code. We're going to
use naive Monte Carlo and we're going to also use importance sampling. I have
plotted here in the black line we have h of X times pi of X. Our target
distribution is exponential one and the h of X is e to the negative X plus
cosine X. From naive Monte Carlo, we know that this expected value is equal to the
integral of e to the negative X plus cosine X times the target distribution
PDF. That's just e to the negative X because that's the PDF of exponential
one. Basically we want to estimate the integral of this, which you could
simplify that to be e to the negative 2 X plus cosine X. So we want to estimate
this integral and I plotted h of x times pi of X here in the black line. If we
just did naive Monte Carlo, we would be sampling from a repeated sample from the
exponential one distribution, which is our target distribution. I plotted
here exponential one in pink. This is exponential one from zero to two point
five. You could estimate the expected value of e to the negative X plus cosine
X using just your target distribution and that would be naive Monte Carlo.
In naive Monte Carlo let's say you sample 10,000 samples from your target
distribution exponential one, then you would just plug in those values of the
sample into e to the negative X plus cosine X. You plug in your sample
where the Xs are. Then you could just take the sample mean of that--and I will
show that in the R code--to get an estimate of your expected value. That
would be the naive Monte Carlo method. However, an alternative method would be
to use importance sampling. So in important sampling, here I have plotted
exponential two. Exponential two is the orange line exponential one is the
pink line and if we compare this which I'm going to say is g of X, our proposal
distribution. The proposal distribution, if we sample from that and
do importance sampling, it's going to give us a better estimate of the
expected value because if you look at the support of h of X pi of X, which is
what we want to estimate from, the orange line more closely follows the black line
compared to the pink line. Because the orange line more closely follows the
support of h of X pi of X, it's going to give us a better estimator when we use
Monte Carlo methods. If you have a question of how do I find a good
proposal distribution, you want to plot either in a graphing calculator or in a
computer h of X pi of X over the support of your target distribution and you want
to find a proposal distribution that closely follows the plotted line. You
want it to closely follow h of X pi of X. You can see that in the beginning
before the exponential decline, the orange line more closely follows the
black line and still here when the orange line is above the black line it's
still closer than the pink line. In both cases when it's below the target
distribution and when it's above the target distribution, the orange line is
still better than the pink line and another thing to note with importance
sampling is that the orange line--so your proposal distribution can go above and
it can go below your black line. That's okay either way. It's okay if it's
not covering completely your distribution and it's okay if it's above
it as well. The main thing is you want to make sure that your proposal
distribution has the same support as your target distribution. Exponential
one goes from 0 to infinity and the proposal distribution exponential two will
also go from 0 to infinity, so it has the same support. It's a valid proposal
distribution. I'm going to do this example above in R but to write out the
algorithm, the algorithm using importance sampling is going to be: we want to draw
iid samples X1 to some value Xm. That could be a hundred samples; it could
be 10,000 samples; just some number m from the proposal distribution g of X.
Then we're going to calculate what are called important weights. We want to
calculate important weights and the importance weights are defined as w(i)--I'm
going to call them w--is equal to pi X of i divided by g X of i. For each random
sample that we get in step one, we're going to plug that into the target
distribution and divide by the proposal distribution for each sample and that's
going to give us weights for every sample. And then finally we just estimate
the properties that we want. If we want to estimate mu, we're going to
estimate mu using the following formula: mu hat is equal to--and then you're
basically for each sample just going to multiply the weight by h of X. So we have
the first weight multiplied by h of X1 and then after you multiply we're going
to add those all together. The weights times the h function for each sample and
then you divide by the number of samples and that would give your estimated mean
of your target distribution. Then let's go ahead and do this example in
R. Okay so in this window, I have R opened up and I'll try to put the sample
code in the description box below. The first lines of code just re-plot the
graph that I showed. We have the h of X times pi of X in the black line. We
have our proposal distribution exponential two in orange and then the
target distribution in pink. So that's just redoing the plots and then the
first algorithm that I have is the importance sampling algorithm. In this
example I'm going to start off by creating a thousand iid random samples.
So m is going to be 1,000 in this case and then the first step is we want 1,000
iid random samples from our proposal distribution. That's going to be a
thousand draws from exponential two. I'm going to do that in this line of code.
After I run that I have 1,000 samples from exponential two and then in this for
loop I'm going to create the important weights for each X(i) that I have. I'm
going to plug that into the target distribution, which is exponential one,
and then divide that from the proposal distribution, which is exponential two.
After I run this line of code in the for loop I'm going to have one thousand
weights and then in the same line of code I say to multiply each weight by h
of X. In this example h of X was e to the negative x plus cosine X, so I'm just
creating the numerator for the estimate of mu that I want to have. And also
before I run that for loop, I create an empty vector for each of the weights and
then I'm also creating an empty vector for the mu. Once that's created I can run
my for loop and then all I have to do to estimate the mu is sum up all the values
of my vector mu and divide them by m, the number of samples. I can just do that
by running mean of mu and that will give me the sample mean. When I run that,
it's going to give me that the sample mean is 1.157226. That's going to
be the estimate of our expected value of h of X. The expected value of e to the
negative X plus cosine X when your target distribution is exponential one.
That's our important sample estimator and also I got the standard error as
well. That's going to be the square root of the variance of mu divided by
the sample size and when I run that our standard error is .008.
As a reference the integral that we're trying to estimate our expected value
does not have a closed form solution but I looked up the calculation in
Wolfram Alpha. Approximately the correct answer should be one point one
six and so our estimator was one point one five seven so we're in the ballpark.
Using just 1,000 samples, that's not too bad. And then at the bottom for this line
of code, I just wanted to compare the important sample method to the naive
Monte Carlo method. To do the naive Monte Carlo it's a couple extra lines of
code that's very simple to run. For naive Monte Carlo we're going to sample
directly from the target distribution. Our X is directly sampled from
exponential one and then we just do the sample mean of h of X and that would
give us an estimator of 1.1437 and I also calculated
the standard error. The standard error would be .027. If the
true answer should be approximately 1.16 our important sampling estimator
was one point one five seven and the naive Monte Carlo method gave us an
estimator of one point one four three seven. You can see the importance
sampling gave us a better estimator using 1,000 samples than the naive Monte
Carlo method and it also gave us variance reduction like I said it would.
It gave us a smaller standard error by a lot compared to the naive Monte
Carlo method and you can compare this for different numbers of iterations as
well. If I change this to 10,000 then the estimator for the mean should get
better. You see it went from one point one five seven to one point one six zero
using importance sampling and the standard error gets smaller as well if
you increase the sample size. And the naive Monte Carlo method also improves
when you increase the number of iterations. So that's it for this
video and thank you guys for watching~~
*chiptune music*