Monte-Carlo integration in spherical coordinates

Question

I tried to implement a simple Monte-Carlo integration in in 6 dimensions. The integral only covers two 3D spheres, in the following the coordinates for the spheres are labeled r1 and r2. When using cartesian coordinates and ignoring anything outside of the spheres the integration works fine. Using spherical coordinates fails, when the integrand depends on angles between r1 and r2.

It seems that the r1 and r2 are likely to point in the same direction, whereas I would expect there alignment to be completely random.

The transformation used to generate spherical coordinates can be found here: http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html

The following code illustrates this with two different integrands and Monte-Carlo integration based on cartesian and spherical coordinates:

import numpy as np
from math import *

N=1e5
R=2
INVERT = False # if this is set to True, the results are correct
INTEGRAND = None


def spherical2cartesian(r, cos_theta, cos_phi):
    sin_theta = sqrt(1-cos_theta**2)
    sin_phi = sqrt(1-cos_phi**2)
    return np.array([r*sin_theta*cos_phi, r*sin_theta*sin_phi, r*cos_theta])

# Integrands (cartesian)

def integrand_sum_sqr(r1,r2):
    r1=np.linalg.norm(r1)
    r2=np.linalg.norm(r2)
    if r1>R or r2>R:
        return 0.0;
    return r1**2 + r2**2

def integrand_diff_sqr(r1,r2):
    delta = r1-r2
    r1=np.linalg.norm(r1)
    r2=np.linalg.norm(r2)
    if r1>R or r2>R:
        return 0.0;
    return np.linalg.norm(delta)**2

# integrand for spherical integration (calls cartesian version)

def integrand_spherical(r1,r2):
    # see the following link for the transformation used: http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html
    r1 = spherical2cartesian((3*r1[0])**(1.0/3.0), r1[1], cos(r1[2]))
    r2 = spherical2cartesian((3*r2[0])**(1.0/3.0), r2[1], cos(r2[2]))
    if INVERT:
        s = (-1)**np.random.choice(2,6)
        r1=s[0:3]*r1
        r2=s[3:6]*r2
    return INTEGRAND(r1,r2)

# monte carlo integration routines

def monte_carlo(name,func,samples,V):
    results=np.array([func(x[0:3],x[3:6]) for x in samples])
    avg = results.mean()*V
    std = results.std(ddof=1)*V/sqrt(len(samples))
    print name,": ",avg," +- ",std

def monte_carlo_cartesian():
    V=(2*R)**6
    samples = np.random.rand(N,6)
    samples = R*(2*samples-1)
    monte_carlo('cartesian',INTEGRAND,samples,V)

def monte_carlo_spherical():
    V=(4.0/3.0*R**3*pi)**2
    samples = np.random.rand(6,N)
    samples = np.array([R**3*samples[0]/3.0, 2*samples[1]-1, 2*pi*samples[2], R**3*samples[3]/3.0, 2*samples[4]-1, 2*pi*samples[5]])
    samples = samples.T
    monte_carlo('spherical',integrand_spherical,samples,V)

# run all functions with all different monte carlo versions
print "Integrating sum_sqr, expected:",32.0/15.0*R**8*pi**2
INTEGRAND=integrand_sum_sqr
monte_carlo_cartesian()
monte_carlo_spherical()

print "Integrating diff_sqr, expected:",32.0/15.0*R**8*pi**2
INTEGRAND=integrand_diff_sqr
monte_carlo_cartesian()
monte_carlo_spherical()

Typical output:

Integrating sum_sqr, expected: 5390.11995025
cartesian :  5406.6226422  +-  29.5405030567
spherical :  5392.72811794  +-  5.23871574928
Integrating diff_sqr, expected: 5390.11995025
cartesian :  5360.20055643  +-  34.8851044924
spherical :  4141.88319573  +-  9.69351527901

The last integral is obviously wrong. Why are r1 and r2 correlated? How can I fix that?

Monte-Carlo integration in spherical coordinates

Answers (1)

Related Questions