Understanding Monte Carlo and Importance Sampling

Question

I'm having some trouble understanding Monte Carlo method with and without Importance Sampling. Please alert me of any mistakes! We desire $I=\int_{a}^{b}f(x)dx$.

Simple Monte Carlo, hit'n'miss: By taking the quantity of randomly and uniformly generated points under the function and dividing this by the total number of generations, one approaches the ratio of the integrand to the domain area, A, in which the random number were generated.$$\frac{I}{A}=\lim_{n\to\infty}\frac{n_{hit}}{n}$$ One could code this process along the following lines ($y_{max}$ and $y_{min}$ are the bounds that help define A):

hits = 0
do i = 0, n
    x = ran()*(b-a)+a
    y = ran()*(ymax-ymin)+ymin
    if (y <= f(x)) hits = hits + 1
end do
I = (b-a)*(ymax-ymin)*(hits/n)

Importance Sampling: Accuracy can be improved by changing the integrand (to say $g(x)$) such that the variance is decreased. To preserve the integrals' value however, we must offset this change by altering the distribution of random numbers. Say $p(x)$ is this normalized distribution. $$I=\int_{a}^{b}f(x)dx=\int_{a}^{b}\bigg(\frac{f(x)}{p(x)}\bigg)p(x)dx=\int_{a}^{b}g(x)p(x)dx$$ In our Simple Monte Carlo discussion, $p(x)= (b-a)^{-1}$. Loosely speaking, we'd like to pick $p(x)$ such that $g(x)$ is as constant as possible.

My Questions:

• Is my understanding of the two forms of Monte Carlo Integration correct?
• In a code utilizing the Importance Sampling, is there anything else to gain aside from improved accuracy?
• What might the pseudo-code of the Importance Sampling method look like?

This last question is my greatest concern. Is the implementation as simple as writing $g(x)$ in place of $f(x)$, and $x=rand()\cdot(b-a)+a$ for whatever the new distribution gives by inverting, $$rand()=\int_{a}^{x}p(x')dx'$$ Is $y$ still generated uniformly random?

Answer 1

So this "hits" approach to Monte Carlo is taught quite frequently and is actually rather bad for a number of reasons. The better way to do it is to just draw uniform values of $x$ and then average the corresponding values of $f(x)$. This gives you an estimate of $\frac{1}{b-a} \int_a^b f(x) dx$ which you can multiply by $b-a$ to get an estimate of the integral itself. This method converges faster and is more similar to most other methods than this "hits" method. The "hits" method is usually only good when you can't evaluate $f$ at all, for example if you are trying to evaluate $\pi$ as the area of the unit circle without using any operations like fractional exponents etc. (In this example, $\pi$ is $4$ times the limit of the fraction of times that $x^2+y^2 \leq 1$ when $x,y$ are independent and uniform on $(0,1)$.)

As for importance sampling, it is also going to draw random values of $f(x)$ and average them, but not uniformly. You instead want to draw more numbers where $|f|$ is larger and fewer where it is smaller. In the process you lower the variance of your generator, and since the error estimate is roughly $\frac{\sigma}{\sqrt{n}}$, that improves the error, sometimes quite dramatically. This is also essential when doing Monte Carlo integration on infinite intervals, because such domains have no uniform distribution.

It is not so trivial to write general pseudocode for importance sampling because it requires you to be able to sample from this distribution $p$ which at the moment is arbitrary. The most general approach would be to define the quantile function $Q$. The quantile function $Q$ is the inverse of the CDF $F$ if that inverse exists, and otherwise it is "the best generalization". (You can look it up or ask for clarification in a comment.) In this approach you take u=rand() and average f(Q(u))/p(Q(u)). But this is not the only way to sample (for example, this is not how we sample from the normal distribution in practice).

There is no analogue of the "y" in your "hits" method in importance sampling. You just generate $x$ and then average up the values of $f(x)$ you get (weighted according to $p(x)$, if you're doing importance sampling). Thus the pseudocode looks something like:

I=0
do i=0,n
  x=samplefromp()
  I=I+f(x)/p(x)
end do
I=I/n

(assuming that "do i=0,n" runs the loop n times). Don't get too locked into a particular implementation of "samplefromp()"; at this level you should think of that as a black box to be implemented later.