用scipy积分tplquad求多元高斯三重积分

Question

我试着做下面的三重积分：

f（v）是一个3变量的高斯概率密度函数。在

合成速度V的大小必须小于某个Vmax（因此Vx、Vy和Vz的范围都可以从0或-Vmax到+Vmax，但是如果Vx=Vmax，那么Vy和Vz必须为零。在

现在我们可以取sigma=1。我现在也忽略了积分中除以v的情况，所以我只是对f（v）积分。在

我一直想用整合tplquad函数，但是要得到一个数量级~1e-7的答案，如果Vmax足够大（我用Vmax=500），积分应该接近1，因为在整个（速度）空间的概率是1。在

这是我的代码：

from scipy.integrate import tplquad
from numpy import pi, exp, sqrt

def func(Vz,Vy,Vx):
    return sqrt( (1/(sigma*2*pi)**(3/2)) ) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) ) )

def Vymax(Vx):
    return sqrt(Vmax**2 - Vx**2)

def Vzmax(Vx,Vy):
    return sqrt(Vmax**2 - Vx**2 - Vy**2)

Vmax = 500
sigma = 1

integral = tplquad(func,0,Vmax,lambda a: 0, Vymax,lambda a, y: 0, Vzmax)
print(8*integral[0])

输出：

^{pr2}$

Vymax和Vzmax函数用于将Vy和Vz值保持在范围内，以便合成速度不超过Vmax。我在tplquad中为0使用lambda函数，因为tplquad需要一个函数作为第4个参数的输入。在

我看不出哪里出了问题，但我肯定我错过了一些简单的事情，或者完全是愚蠢的。在

如果tplquad不是解决这个问题的合适函数，还有其他选择吗？我甚至很乐意编写我自己的函数，但不确定什么是最好的方法-我尝试了蒙特卡罗方法，但没有完全弄明白。我不能仅仅把组件分开，因为最终我需要有一个偏移Vx+Voffset，这样它就不再以原点为中心了。我尽我所能在这里搜索了很多，发现了this，但是它没有恰当地描述如何做一个多元高斯函数，因为问题是（本来是）关于单变量高斯的。在

非常感谢你的帮助。在

谢谢。在

Answer 1

您的公式multivariate normal distribution不正确。你有sqrt调用，和指数系数幂的1/2因子。指数中还缺少系数1/sigma**2。在

这是一个修改版本，使用sigma的方式与问题中显示的公式相匹配：

def func(Vz,Vy,Vx):
    return (1.0/(sigma**2*2*pi)**(3./2)) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) )/sigma**2 )

添加print语句以打印sigma的值后，我得到：

^{pr2}$

编辑脚本并再次运行：

In [29]: run tplquad_question.py
sigma = 2.5
8*integral[0] = 1.00000000927

以下是完整的tplquad_question.py：

from __future__ import print_function, division

from scipy.integrate import tplquad
from numpy import pi, exp, sqrt

def func(Vz,Vy,Vx):
    return (1.0/(sigma**2*2*pi)**(3./2)) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) )/sigma**2 )

def Vymax(Vx):
    return sqrt(Vmax**2 - Vx**2)

def Vzmax(Vx,Vy):
    return sqrt(Vmax**2 - Vx**2 - Vy**2)

Vmax = 500
sigma = 1

integral = tplquad(func,0,Vmax,lambda a: 0, Vymax,lambda a, y: 0, Vzmax)

print("sigma =", sigma)
print("8*integral[0] =", 8*integral[0])