球坐标系下的蒙特卡罗积分

Question

我试图在6维中实现一个简单的蒙特卡罗积分。积分只覆盖两个三维球体，在下面的球体坐标标注为r1和r2。 当使用笛卡尔坐标而忽略球体之外的任何东西时，积分工作得很好。当被积函数依赖于r1和r2之间的角度时，使用球坐标是失败的。

看来r1和r2可能指向同一个方向，而我认为它们的排列完全是随机的。

用于生成球坐标的变换可以在这里找到：http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html

以下代码通过两个不同的被积函数以及基于笛卡尔坐标和球坐标的蒙特卡罗积分来说明这一点：

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()

典型输出：

^{pr2}$

最后一个积分显然是错误的。 为什么r1和r2相关？我该怎么解决呢？

Answer 1

上面代码的问题在下面的行中

sin_phi = sqrt(1-cos_phi**2)

这只会产生阳性结果，而sin(phi)则会产生phi > pi的阴性结果