I'm trying to approximate $\pi$ using Monte Carlo integration; I am approximating the integral
$$\int\limits_0^1\!\frac{4}{1+x^2}\;\mathrm{d}x=\pi$$
This is working fine, and so is estimating the error (variance), $\sigma$. However, when I then try to use importance sampling with a Cauchy(0,1) distribution, things start to go wrong:
$$\frac{1}{n}\sum\limits_{i=0}^n\frac{f(x_i)}{p(x_i)}=\frac{1}{n}\sum\limits_{i=0}^n\frac{\frac{4}{1+x^2}}{\frac{1}{\pi(1+x^2)}}=\frac{1}{n}\sum\limits_{i=0}^n\frac{4\pi(1+x^2)}{1+x^2}=\frac{1}{n}\sum\limits_{i=0}^n4\pi=4\pi$$
Obviously something's wrong, since the mean is computed independently of the random variables I generate. Where is this going wrong? Is the distribution too close to $f$?
Best Answer
This is quite a common error when doing Monte Carlo integration. The support of the random variable you choose must be equal to the range of integration. Though the Cauchy distribution has support on $\mathbb{R}$ we can adapt it slightly to make it work here.
Note that: $\int_0^1 \frac{1}{\pi(1+x^2)} dx = \frac{1}{\pi}\tan^{-1}(1)=\frac{1}{4}$
Thus $g(x) = \frac{4}{\pi(1+x^2)}$ for $x\in(0,1)$ is a probability density with support on $(0,1)$.
using this we have
$\frac{1}{n}\sum_{i=0}^n \frac{f(x_i)}{g(x_i)} = \pi$
This is not a problem! Since $g(x) = \frac{1}{\pi}f(x)$ you have found exactly the right probability distribution to use to evaluate $\int_0^1f(x)dx$! The error is zero for any sample, regardless of the size.
Note however, I needed to know how to integrate $\int_0^1\frac{1}{\pi(1+x^2)}$ in the first place to form the probability distribution. Thus for arbitrary functions it is impossible to get this situation.
Also note the Monte Carlo is not a very good way of approximating integrals in general. Far better deterministic methods are quadrature rules.