Rrerku

I am following these lecture slides on Monte Carlo integration with importance sampling. I am just implementing a very simple example: $int_0^1 e^xdx$. For the importance sampling version, I rewrite $int_0^1 e^xdx = int_0^1 e^x/p(x)cdot p(x)dx$ where $p(x) = 2.5x^1.5$. Then

$$hatI = frac1Nsum_j=1^N fracf(x_j)p(x_j),$$

where $x_j$ are sampled from $p(x_j)$ (I use an inverse transform method here). For the variance, I have $sigma_I^2 = hatsigma_I^2/N$ and

$$hatsigma_I^2 = frac1N sum_j=1^N fracf(x_j)^2g(x_j)^2 - hatI^2.$$

I know I should expected the variance to decrease with importance sampling, but a plot of the variance with $N$ shows that not much happens. Can anyone explain to me what I'm doing incorrectly? I'm not sure how the they are able to achieve such a drastic decrease in variance in the lecture slides.

This is a good illustration of the dangers of importance sampling: while
$$int_0^1 frace^xp(x), p(x)textd x = int_0^1 e^x textd x = I$$
shows that $hatI_N$ is an unbiased estimator of $I$, this estimator does not have a finite variance since
$$int_0^1 left(frace^xp(x)right)^2, p(x)textd x = int_0^1 frace^2x2.5 x^1.5 textd x = infty$$
since the integral diverges in $x=0$. For instance,

> x=runif(1e7)^1/2.5
> range(exp(x)/x^1.5)
[1] 2.718282 83403.685972

shows that the weights can widely differ. I am not surprised at the figures reported in the above slides since

 > mean(exp(x)/x^1.5)/2.5
 [1] 1.717576
 > var(exp(x)/x^1.5)/(2.5)^2/1e7
 [1] 2.070953e-06

but the empirical variance is rarely able to spot infinite variance importance sampling. (The graph shows that both the standard Monte Carlo estimate and the importance sampling version see the empirical standard deviation is decreasing as $N^-1/2$.)