Finding the approximate value of improper integral using the Monte-Carlo method 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.
 A: 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$.
