Finding the approximate value of improper integral using the Monte-Carlo method

improper-integralsintegrationmonte carloprobability distributions

Find the approximate value of the improper integral
$$
\int_{-3}^{\infty} \left( \int_{0.5}^{\infty} \frac{2+\sin(x+y)}{e^{0.4x}+0.4y^2} \,dx \right) \, dy
$$
using the Monte-Carlo method.

I have managed to define the necessary functions in R in order to apply the MC method but the problem lies in the function under the integral itself. I have tried to use well-known distributions (exponential and normal) to rewrite the function and get an idea from which distributions should I generate $x$ and $y$ but due to lack of luck in doing that it seems like a wrong direction, perhaps there is an easier way.

So any help, hints, tips and tricks on how to solve the problem would be greatly appreciated.

Best Answer

Follows a python script which calculates an approximation to the integral. The script is written with the purpose of explain the Monte Carlo method.

import numpy as np

xmin = 0.5
xmax = 100
ymin = -3
ymax = 100
zmin = 0
zmax = 3
N = 1000000

def ok(x,y,z):
    f = (2.0 + np.sin(x + y))/(np.exp(0.4*x) + 0.4*y*y)
    if f > z:
        return True
    else:
        return False

def guess():
    x = np.random.uniform(xmin,xmax,1)
    y = np.random.uniform(ymin,ymax,1)
    z = np.random.uniform(zmin,zmax,1)
    return (x,y,z)


s = 0;
for i in range(N):
    (x,y,z) = guess()
    if ok(x,y,z):
        s = s + 1
    

vol = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(s/N)

print(vol)

NOTE

As

$$ \frac{2+\sin(x+y)}{e^{0.4x}+0.4y^2}\le \frac{3}{e^{0.4x}+0.4y^2} $$

we can estimate the amount left by cutting the integration range: thus as

$$ \int \left( \int \frac{3}{e^{a x}+ay^2} \,dx \right) \, dy = \frac{3 \log \left(a y^2 e^{-a x}+1\right)}{a^2 y}-\frac{6 e^{-\frac{a x}{2}} \tan ^{-1}\left(\sqrt{a} y e^{-\frac{a x}{2}}\right)}{a^{3/2}} $$

we can conclude that with $x_{max}\le x\lt \infty$ and $y_{max}\le y \lt \infty$ the integration "forgets" an amount to compute in the order of $10^{-8}$. Now the accuracy depends on $N$.

Related Question