# Implementation of Monte-Carlo Integration

After reading the Wikipedia page for Monte-Carlo integration, I have understood the basic idea but I am having trouble implementing it for a general case.

The integration that I am trying to do is $$\int_{-\infty}^{\infty} \frac{\cos x}{x^2}- \frac{\sin x}{x^3} dx$$

This integral is equal to $$-\frac{\pi}{2}$$, but I do not see how to do it with the Monte-Carlo method. What is the algorithmic procedure in this case?

From Wikipedia:

The naive Monte Carlo approach is to sample points uniformly on Ω[...]

There is an implicit assumption here that a uniform distribution on $$\Omega$$ exists. It is well-known that such a distribution does not exist on any unbounded subinterval of $$\mathbb{R}$$. Therefore, using basic Monte Carlo will not work for this problem. One option is to do a $$u$$-substitution so that $$u$$ is integrated over a finite interval then perform uniform Monte Carlo sampling on that interval. Since the integrands are highly oscillatory, this will likely take a while to converge.

• That works, but with the substitution u/(1-u), I was more than 33million samples in and I still only had four significant figures of pi. Commented Apr 1 at 19:20
• Monte Carlo error decreases at a rate $\propto N^{-1/2}$. An error of $10^{-4}$ corresponds to $N\approx10^8$, so your results are to be expected. Commented Apr 1 at 20:00
• @JamesK MC integration is not very efficient in low number of dimensions. It becomes more efficient than other methods in high number of dimensions where the number of points in a regular grid grows like $N^d$. Commented Apr 1 at 20:38
• Note that the integrand is an even function; that should let you halve the number of samples. Commented Apr 2 at 4:55

Note that the integrand is an even function, so we need to evaluate

$$I = 2 \int_0^{\infty} f(x)/x^2 dx,$$

where $$f(x) = \cos(x)-\sin(x)/x$$

Let's write it as

$$I/2 = \int_0^{1} f(x)/x^2 dx + \int_1^{\infty} f(x)/x^2 dx$$

The difficult part is the second half, and one way to deal with it is to draw the Monte-Carlo samples from the distribution $$P(\xi)=1/\xi^2$$ for $$\xi \in [1,\infty]$$. Then the improper integral is approximated by

$$\frac{1}{N} \sum_{k=1}^N f(x_k),$$

where $$N$$ is the number of samples.

This approach may improve the accuracy, i.e., reduce the cost of the calculation.