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

p5cysglq  于 7个月前  发布在  其他
关注(0)|答案(1)|浏览(43)

我正在尝试执行以下三重积分:
x1c 0d1x的数据
f(v)是一个3变量高斯概率密度函数。
合成速度V的大小必须小于某个Vmax(例如,Vx、Vy和Vz的范围可以从0或-Vmax到+Vmax,但如果Vx = Vmax,则Vy和Vz必须为零)。
现在我们可以取sigma = 1,我也会忽略积分中除以v的部分,所以我只对f(v)积分。
我一直在尝试使用scipy.integrate tplquad函数,但不断得到一个数量级为~ 1 e-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需要一个函数作为第4个参数的输入。
我看不出我错在哪里,但我肯定我一定是错过了一些简单的东西,或者是完全愚蠢的。
如果tplquad不是解决这个问题的正确函数,有替代方法吗?我甚至很乐意编写自己的函数,但不确定最好的方法是什么-我尝试了蒙特卡罗方法,但不能完全弄清楚。我不能只是分离出组件,因为最终我需要有一个偏移Vx +我已经尽可能多地搜索了这里,遇到了this,但它没有正确描述如何做多变量高斯,因为问题是(意味着)关于一个单一变量高斯。
任何帮助非常感谢。

nue99wik

nue99wik1#

您的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的值后,我得到:

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

相关问题