CODEWITHSUNDEEP
pythonmathscipynumerical-integration
ProgSnob
ProgSnob

Reputation: 493

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?

Upvotes: 0

Views: 286

Answers (2)

Andrei
Andrei

Reputation: 2665

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

Upvotes: 1

dmuir
dmuir

Reputation: 4431

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)

Upvotes: 0

Related Questions