5
$\begingroup$

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$?

  • 0
    I think there's some confusion over the range of your integration and thesupport of the Cauchy distribution; it's supported on (-\infty, \infty), not on (0,1). Your sum (presuming that x_i = i and the x's in the second formula are supposed to be i's) covers (0, \infty), so you're not computing the integral from 0..1 that you think you are. Can you be a bit more explicit about your sampling process?2010-09-23
  • 1
    "Monte Carlo is an *extremely* bad method; it should be used only when all alternative methods are worse." ~ Alan Sokal2010-09-23

1 Answers 1

7

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.

  • 1
    I've been told that Monte Carlo methods get more useful the higher the dimension.2010-09-23
  • 0
    Here are my thoughts on the multi-dimensional setting: Monte Carlo methods in practice are going to give a larger error than their quadrature counterparts. However, in high dimensions quadrature rules require more evaluations of the integrand and this may be unacceptable (if evaluations are expensive).2010-09-23
  • 0
    My take on it is: Monte Carlo methods only get more useful the higher the dimension because quadrature rules become less so. (As you can tell, I think Monte Carlo method's popularity are more down to their simplicity rather than their usefulness.)2010-09-23
  • 1
    (Quasi-)Monte Carlo methods are not as affected by the "curse of dimensionality" as deterministic quadrature rules. That's the only reason we're even interested in Monte Carlo in practical work, in the absence of any exploitable structure (e.g. symmetry) in the integrand.2010-09-23
  • 1
    The other place where Monte Carlo reigns supreme is in integration over complicated regions; imagine trying to properly formulate an integral over, e.g., the intersection of three tori, or even something as simple as the interior of a convex polyhedron. As long as you have some mechanism for deciding point-in-region, and some reasonable bounding sphere/box/'easy' shape for your region, though, you can just rejection-sample to limit your Monte Carlo samples to the interior of the region and then normalize the result by dividing out by the sum of sample weights.2010-09-23
  • 0
    This sort of thing, incidentally (integration over complicated regions) happens more rarely mathematically but actually shows up quite a bit in high-end computer graphics, where the notion of stratified sampling was also brought to the fore.2010-09-23
  • 0
    If you take a Monte Carlo integral with $n$ points, the error is $O(\sqrt{n})$, *regardless of dimension*. It does not require an absurdly large number of dimensions for this to surpass every technique which would be useful in low dimensions.2011-09-14