# Monte Carlo to evaluate infinite integral on R

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.

• I think you made an error in your change of variable. In any case, the better way to evaluate this integral with Monte Carlo would be to take the average of $f(X)$ where X has PDF $e^{x/5}/5$ and $f$ is chosen appropriately. In terms of a change of variable this amounts to taking $x/5=-\ln(u)$. – Ian Sep 17 '18 at 1:48
• Assuming your transformed integral is in fact correct (and I see that you are at least on the right track actually because the exponential saves the day at $u=0$), you still have a serious problem in that your integrand is a wildly oscillatory function, so the random variable whose expectation you are trying to calculate has an extremely high variance, causing very slow convergence of the Monte Carlo integration. By contrast, with $x/5=-\ln(u)$, you will find that the variance is far smaller. – Ian Sep 17 '18 at 2:25
• The reason for this is that although you are giving more weight through your substitution to smaller values of $x$ than larger ones, you are still giving way more weight to the larger ones than the actual integral does. Making a substitution like $x/5=-\ln(u)$ respects the fact that the integrand decays exponentially fast at infinity, which dampens out the wild oscillation. (There are still infinitely many sign changes near $u=0$, but there is no way to avoid that.) – Ian Sep 17 '18 at 2:27

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.

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