Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
    
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? –  Steven Stadnicki Sep 23 '10 at 0:25
1  
"Monte Carlo is an extremely bad method; it should be used only when all alternative methods are worse." ~ Alan Sokal –  J. M. Sep 23 '10 at 14:14

1 Answer 1

up vote 6 down vote accepted

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.

share|improve this answer
1  
I've been told that Monte Carlo methods get more useful the higher the dimension. –  Qiaochu Yuan Sep 23 '10 at 7:53
    
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). –  alext87 Sep 23 '10 at 8:43
    
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.) –  alext87 Sep 23 '10 at 8:53
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. –  J. M. Sep 23 '10 at 9:57
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. –  Steven Stadnicki Sep 23 '10 at 22:17

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.