Computing integrals in n dimensions with Monte Carlo methods

Let's try to find the volume of a unit sphere in n dimensions (the n-ball) numerically...

The n-ball of radius $R$ in $n$ dimensions has volume $V_{n}(R)$ given by:

$$ V_{n}(R)={\frac {\pi ^{\frac {n}{2}}}{\Gamma \left({\frac {n}{2}}+1\right)}}R^{n} $$

so for $R=1$ we have the exact answer given by (let's call it $U(n)$):

$$ U(n)={\frac {\pi ^{\frac {n}{2}}}{\Gamma \left({\frac {n}{2}}+1\right)}} $$

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import special
In [2]:
def vol_sphere(n):
    "Return the volume of the unit ball in n dimensions"
    return np.pi**(n/2)/ special.gamma(n/2+1)

Let's try to find this as an integral, since a volume is always the integral of a region.

In [3]:
def vol_sphere_mc(n, reps=10_000):
    points = np.random.uniform(-1, 1, size=(n, reps))
    radius = (points**2).sum(axis=0)
    
    return (2**n)*((radius<=1).sum()/reps)
In [8]:
reps = 10_000_000
for n in [2, 3, 4, 10, 20]:
    ex = vol_sphere(n)
    mc = vol_sphere_mc(n, reps)
    print(n, 2**n, ex, mc, 100*abs(ex-mc)/ex )
2 4 3.141592653589793 3.1410508 0.01724773544953504
3 8 4.188790204786391 4.1901416 0.032262184247470065
4 16 4.934802200544679 4.9330416 0.03567722622164131
10 1024 2.550164039877345 2.54208 0.31700077920218817
20 1048576 0.02580689139001405 0.0 100.0