Scipy tplquad 语法

Question

我希望在 python 中使用 scipy 中的 tplquad 函数进行以下积分：

到目前为止，我使用了以下语法：

def func(z, y, x):
        return (x**(b1 + x1 - 1))*(y**(b2 + x2 - 1))*(z**(b3 + x3 - 1))*((1-x-y-z)**(x4+b4-1))

t1 = tplquad(func, 1/2, 1, lambda y: 0, lambda y: 1-y, lambda y, z: 0, lambda y, z: 1 - y - z)

其中func是我们要集成的函数：f(p_1,p_2,p_3) ，z 为 p_3，y 为 p_2，x 为 p_1

谁能告诉我在这里使用 tplquad 的正确语法应该是什么？

编辑： 我在限制方面遇到了问题 - 如果限制，如果最内层积分的限制是 0 到 p1，我将使用什么作为限制在 lambda 函数方面？

Answer 1

这不是您问题的答案，但对于评论来说太长了，我认为它可能会有帮助。

您可以根据 Beta function.

进行解析积分

我将调用 k（对于 b1 + x1 - 1）、l、m、n 的幂。最里面的积分是

x**k * y**l * Integral{ 0<=z<=1-x-y | z**m * ( 1-x-y-z)**n dz}

令 r = 1-x-y。替换 zeta = z/r 将其转换为

   x**k * y**l * r**(m+1) * J

哪里

 J = Integral{ 0<=zeta<=1 | zeta**m * ( 1-zeta)**n dzeta}
   = Beta( m+1, n+1)

一个类似的过程也得到了 y 的积分，我们剩下

K*Integral{ 0.5<=x<=1 | x ** k * (1-x)**(l+1) dx

哪里

K = Beta( m+1, n+1) * Beta( l+1, n+2)

x 积分可以使用不完整的 beta 函数表示为

B(k+1, l+2) - B(0.5; k+1, l+2)

Answer 2

您的解决方案是正确的。积分的解析公式为（在 Mathematica 中计算）

def ftest():
    a = b1 + x1 - 1
    b = b2 + x2 - 1
    c = b3 + x3 - 1
    d = b4 + x4 - 1
    return (gamma(1+a)*gamma(1+b)*gamma(1+d)*(-beta(1+c,3+a+b+d)*betainc(1+c,3+a+b+d,1/2)+(gamma(1+c)*gamma(3+a+b+d))/gamma(4+a+b+c+d)))/gamma(3+a+b+d)

它给出的结果与 tplquad 完全相同。例如，如果 b1=b2=b3=b4=2 和 x1=x2=x3=x4=2，我在这两种情况下都会得到 1.7421194876552e-11

Scipy tplquad 语法

Scipy tplquad syntax

python

math

integration

scipy