Search code examples
pythonmathscipynumerical-integration

Scipy tplquad syntax

I wish to do the following integral in python using the tplquad function in scipy:

enter image description here

So far, I've used the syntax below:

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)

where func is the function we want to integrate: f(p_1, p_2, p_3), z is p_3, y is p_2, x is p_1

Can anyone tell me what the correct syntax for using tplquad here should be?

EDIT: I'm having trouble with the limits - if the limits, if the limits of the inner most integral were 0 to p1, what would I use as the limits in terms of lambda functions?

Solution

  • You solution is correct. The analytical formula for the integral is (computed in 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)

    It gives exactly the same result as tplquad. For example, if b1=b2=b3=b4=2 and x1=x2=x3=x4=2, I get in both cases 1.7421194876552e-11