Scipy tplquad syntax

Question

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

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?

Possible duplicate of [How to use tplquad?](https://stackoverflow.com/questions/40431553/how-to-use-tplquad) — ev-br, Jun 21 '17 at 08:17
@ev-br I've seen the post, but the limits in that integral are constants, and I'm having trouble with the limits only, where my limits are functions. — ProgSnob, Jun 21 '17 at 08:20

score 1 · Accepted Answer · answered Jun 23 '17 at 15:27

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

score 0 · Answer 2 · answered Jun 21 '17 at 11:51

This is not an answer to your question, but it's too long for a comment, and I thought it might be helpful.

You can do the integral analytically, in terms of the Beta function.

I'll call the powers k (for b1 + x1 - 1), l, m, n. The inner most integral is

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

let r = 1-x-y. the substitution zeta = z/r transforms this to

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

where

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

a similar procedure gets the integral over y too, and we are left with

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

where

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

The x integral can be espressed using the incomplete beta function as

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

`B(k+1, l+2) - B(0.5; k+1, l+2)` this term is negative somehow. Is it correctly written? `B(0.5; k+1, l+2)` is the incomplete beta function from 0 to 0.5 for the given range, right? — ProgSnob, Jun 22 '17 at 10:04
@ProgSnob yes, it is the incomplete beta function (not the regularised one). How the difference could be negative is a mystery, as the dirrefence is the integral of a non-negative function over (0,5, 1) — dmuir, Jun 22 '17 at 11:48

Scipy tplquad syntax

2 Answers2