title: Monte-Carlo integration in spherical coordinates
body: <p>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.</p>
  
  <p>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.</p>
  
  <p>The transformation used to generate spherical coordinates can be found here: <a href="http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html" rel="nofollow">http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html</a></p>
  
  <p>The following code illustrates this with two different integrands and Monte-Carlo integration based on cartesian and spherical coordinates:</p>
  
  <pre><code>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&gt;R or r2&gt;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&gt;R or r2&gt;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()
  </code></pre>
  
  <p>Typical output:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The last integral is obviously wrong.
  Why are r1 and r2 correlated? How can I fix that?</p>