-1

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
T-Ray
  • 179
  • Split the integral into 0,1 and 1,$\infty$ make substitution $r=1/s$ in second part. – user121049 Jul 21 '17 at 07:27
  • Nice trick. The rate of convergence is slightly better but still pretty bad with a uniform sampling distribution. – T-Ray Jul 21 '17 at 07:52

1 Answers1

0

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

  • This was actually the first thing that I tried and the rate of convergence was pretty bad. There was a 10% error at 10^7 sampling position. So I was hoping to get a boost by choosing another pdf. The trick that user121049 discussed increased the rate of convergence to a 1% error at 10^7 sampling position but this is still much smaller than I would like. – T-Ray Jul 21 '17 at 17:58