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