[Math] 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.