Monte Carlo Integration at Infinity

Question

I am trying to numerically approximate the following integral using a Monte Carlo simulation.

$I = \int_0^{\infty} dr r^2 \exp(-a r) J_0(q*r)$

I have arbitrarily chosen $a = 1, q = 10$. I am new to Monte Carlo but it seems that you run into trouble for uniform sampling when one of the endpoints goes to infinity. Namely the sampling will tend to only be carried over very large values. A way around this is to define another sampling distribution. My question is, what sampling distribution gives the fastest convergence and how can this be seen from the function?

I have tried this two way. One way with an exponential pdf and another way by splitting the integration from 0 to 1 and from 1 to infinity and making the substitution $s = 1/r$ (as user121049 had instructed). The code and results are below. The rate of convergence for both methods is approximately equal but is accurate only to about 10% after 10^6 samples. So I am hoping to choose a more efficient pdf to increase the rate of convergence. Feedback is greatly appreciated.

# Begin import
import numpy as np
from scipy.special import jv
import sympy as sp

# Define MC integration algorithm
def montecarlo(f, N, a):
  lst = np.full((1, N), 1.0)
  for i in range(N):
    x = -a*np.log(np.random.random())
    lst[0][i] = f(x)
  return np.mean(lst)
#Generate split method
def splitmontecarlo(f, N, a):
  lst = np.full((1, N), 1.0)
  for i in range(N):
    x = np.random.random()
    lst[0][i] = f(x)
  return np.mean(lst)

#Define exact
a, q, b = sp.symbols('a, q, b', real = True, positive = True)
exact = lambda a, q: sp.integrate(a*b**2*sp.exp(1-a*b)*sp.besselj(0,q*b),(b,0,sp.oo))
# Approximate rate of convergence for Exp pdf
a = 1.0
q = 10.0
N = 10000
trials = 500
f = lambda z: 1/a*a*z**2*np.exp(1)*jv(0, q*z)
errorlst = []
for i in range(trials):
  error = sp.N((exact(a, q))-montecarlo(f, N, a))/sp.N(exact(a, q))
  errorlst.append(error)
print('Exp pdf approximate rate of convergence is'), 1/np.mean(errorlst)/float(N)
# Approximate rate of convergence for split
f = lambda z: a*z**2*np.exp(1)*jv(0, q*z)*np.exp(-a*z)+a/(z**4)*np.exp(1)*jv(0, q/z)*np.exp(-a/z)
errorlst = []
for i in range(trials):
  error = sp.N((exact(a, q))-splitmontecarlo(f, N, a))/sp.N(exact(a, q))
  errorlst.append(error)
print('Split approximate rate of convergence is'), 1/np.mean(errorlst)/float(N)

Exp pdf approximate rate of convergence is -0.000929221448661505
Split approximate rate of convergence is -0.000857990058774096

Answer 1

You can write your integral as $I=\frac1{a}E\left(R^2J_0(qR)\right),$ where $R$ is a random variable with pdf $ae^{-ar}$ for $r\in[0,\infty)$, and estimate that expectation as the average over your sample. You can get such a sample by observing $R=a\ln(1/U),$ where $U$ is a random variable with uniform distribution in $[0,1].$