Monte Carlo to evaluate infinite integral on R

Question

I am using Monte Carlo method to evaluate the integral above: $$\int_0^\infty \frac{x^4sin(x)}{e^{x/5}} \ dx $$ I transformed variables using $u=\frac{1}{1+x}$ so I have the following finite integral: $$\int_0^1 \frac{(1-u)^4 sen\frac{1-u}u}{u^6e^{\frac{1-u}{5u}}} \ du $$ I wrote the following code on R:

set.seed (666)

n <- 1e6

f <- function(x) ( (1-x)^4 * sin ((1-x)/x) ) / ( exp((1-x)/(5*x) ) * x^6 )

x <- runif(n)

I <- sum (f(x))/n

But I get the wrong answer, but if I integrate f(x) using the R built-in function and not Monte Carlo I get the right answer.

Answer 1

Let us use for instance sage to get the exact value of the integral (and some numerical approximation of this exact value).

sage: var('x');
sage: integral( x^4*sin(x)*exp(-x/5), x, 0, oo )
4453125/371293
sage: _.n()
11.9935603418325

Now back to R. Indeed, Monte Carlo delivers a number far away from this one, but the R integral is a good approximation.

> set.seed(666)
> n = 1e6
> f = function(x) ( (1-x)^4 * sin((1-x)/x) * exp(-(1-x)/(5*x)) * x^-6 )
> x = runif(n)
> J = sum(f(x))/n
> J
[1] 101.9375
> integrate(f, 0, 1)
11.99356 with absolute error < 0.00071

The question asks for the reason. It is hard to analyze as a human the x sample above, but it is good to notice the big variance of the function after the substitution. For instance, on the following intervals...

> integrate(f, 0, 0.05)
2910.487 with absolute error < 0.34
> integrate(f, 0.05, 0.10)
-4003.36 with absolute error < 0.093
> integrate(f, 0.10, 0.15)
1425.772 with absolute error < 1.7e-11
> integrate(f, 0.15, 0.20)
-306.4369 with absolute error < 3.4e-12
> integrate(f, 0.20, 0.25)
-30.99198 with absolute error < 3.5e-13
> integrate(f, 0.25, 0.30)
8.218838 with absolute error < 9.1e-14

There are big numbers in the first lines of code, and small changes of the limits lead to "relatively big variations",

> integrate(f, 0.05, 0.10)
-4003.36 with absolute error < 0.093
> integrate(f, 0.0501, 0.1001)
-4013.99 with absolute error < 0.092

so even increasing n may not help. Here is a plot of the function under the integral

and note that we have 1e6 random points chosen on the whole interval, only a part of them on the tiny interval with variation in the same big range. This explains the relatively big difference between J (101.9375) and the true value (11.993560341832...) of the integral.

Answer 2

Your Monte Carlo method tries to compute the expectation of a random variable with an extremely high variance, as you can see by plotting the integrand.

A better Monte Carlo method will respect the decay rate of the integrand as $x \to \infty$, which is slightly slower than $e^{-x/5}$ (so far faster than $1/(1+x)$). One of these is to write the integral as

$$I=\int_0^\infty 5x^4 \sin(x) \left ( \frac{1}{5} e^{-x/5} \right ) dx$$

which is to say $E[5X^4 \sin(X)]$ where $X$ is exponentially distributed with mean $5$. This still has some considerable variance, but it is driven by just a logarithmic singularity at $x=0$ (caused by the fact that $x^4$ wasn't bounded in the original integral) so the convergence behavior should be much better.

An even better method is to write

$$I=5^5 4! \int_0^\infty \sin(x) \left ( \frac{x^4 e^{-x/5}}{5^5 4!} \right ) dx$$

Thus it is $E[\sin(X)]$ where $X$ is Gamma distributed with shape parameter $k=5$ and scale parameter $\theta=5$.