We are trying to use Monte Carlo Simulation to find:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
As you have discovered, the infinite limits are problematic. We can do several things to overcome this
Approach 1: A naive approach is to plot the function given that it has a decaying exponential which tends to zero very quickly. If we plot this function, we see:
So, lets try to approximate it as:
$$\displaystyle \int_0^{50} \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
Using your approach (recall, now your code needs to change to generate random numbers $(a, b)$ and with the function of interest) with $n = 10000$, we get:
$$I = 387.256$$
If we calculate the exact integral we get:
$$I = \Gamma \left(\frac{20}{3}\right) \approx 389.035$$
Approach 2: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one.
A good example that works for functions that decrease towards infinity faster than $\frac{1}{x^2}$ is to use:
$$\displaystyle \int_a^b f(x)~dx = \int_{\frac{1}{a}}^{\frac{1}{b}} \dfrac{1}{t^2}f\left(\dfrac{1}{t}\right)~dt$$
With this altered approach, we can now do integrals when one of the limits is infinite
and the other limit is of the same sign. If you need to integrate say something that has $a = 0$, what you do is break the integral into two parts and use Monte Carlo on each one.
If we have infinite limits for both limits, we could split the integral into three or more parts. So, for this problem, we have to split it into two pieces as:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \dfrac 34 x^4e^{-x^{3/4}}\,dx + \int_0^1 \dfrac{1}{t^2} \left(\dfrac 34\right) \left(\dfrac 1t \right)^4 e^{{-(\frac 1t)}^{\frac 34}}~dt$$
Using the Monte Carlo for the first part, we get $I_1 = 0.0647226$ and for the second part, we get $I_2 = 387.551$, which yields:
$$ \displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = 387.616$$
Compare that to our previous and the exact result.
Approach 3: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one. This time, we will transform variables using $u = \frac{1}{1+x}$, yielding:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \frac{3 e^{-\left(\frac{1-u}{u}\right)^{3/4}} (1-u)^4}{4 u^6} ~ du$$
Running the Monte Carlo simulation with $n = 10000$ yields:
$$I = 390.18704821709946$$
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.
Best Answer
Hint: multiplying your answer by $\pi$ gives the actual integral's value.