Scipy tplquad 语法

Question

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

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

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

您的解决方案是正确的。积分的解析公式为（在 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

Answer 2

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

您可以根据 Beta 函数进行分析积分。

我将幂称为 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)