使用 scipyintegrattplquad 计算多元高斯的三重积分

Question

我正在尝试执行以下三重积分：

f(v) 是 3 变量高斯概率密度函数。

合成速度 V 的大小必须小于某个 Vmax（因此 Vx、Vy 和 Vz 的范围可以分别为 0 或 -Vmax 到 +Vmax，但如果 Vx = Vmax，则 Vy 和 Vz 必须为零，对于示例）。

现在我们可以取 sigma = 1。我现在也会忽略积分中除以 v 的部分，所以我只是对 f(v) 进行积分。

我一直在尝试使用 scipy.integrate 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])

输出：

>>> (executing cell "Triple integral a..." (line 15 of "Integral4.py"))
1.4644282619462532e-07
>>>

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

我看不出我哪里出了问题，但我确信我一定错过了一些简单的东西或者完全愚蠢。

如果 tplquad 不是解决此问题的正确函数，是否有替代方法？我什至很高兴编写自己的函数，但不确定最好的方法是什么 - 我尝试了蒙特卡洛方法，但不太明白。我不能只是分离出组件，因为最终我需要一个偏移量 Vx + Voffset，这样它就不再以原点为中心。我已经在这里尽可能多地搜索并遇到了this，但它没有正确描述如何进行多元高斯分布，因为问题是（应该是）关于单个变量高斯分布。

非常感谢任何帮助。

Answer 1

您的多元正态分布公式不正确。您同时拥有

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 )

添加打印语句打印

sigma

的值后，我得到：

In [28]: run tplquad_question.py
sigma = 1
8*integral[0] = 1.00000000196

编辑脚本并再次运行：

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