Integration of Incomplete Beta Functions and Hypergeometric Functions

urcing · #1 $Old$ November 6th, 2009, 02:30 PM

QUESTION

Let

I = double integral {x<y} [F(x)]^r * [F(y)]^s * [F(x)] * [1 - F(y)] dx dy

Verify/Show that for the 5 parameter Generalized Lambda Distribution (GLD),

I =

(a^2)(b^2)
--------------------------------------------
(r + b + 1)(r + s + 1 + 2b)(r + s + 2 + 2b)

abcd Gamma(d) Gamma(2 + r)
+ --------------------------------------------
(s + b) (s + 1 + b) Gamma(2 + r + d)

abcd Gamma(d) Gamma(r + s + 2 + b)
- --------------------------------------------
(s + b) Gamma(r + s + 2 + b + d)

abcd Gamma(d) Gamma(r + s + 3 + b)
+ --------------------------------------------
(s + b + 1) Gamma(r + s + 3 + b + d)

abcd Gamma(d + 1) Gamma(r + s + 2 + b)
+ --------------------------------------------
(r + b + 1) Gamma(r + s + 2 + b + d)

(c^2)(d^2) Gamma(1 + 2d) Gamma(r + 2) * F(-2,1+d,1+2d;2+d,r+3;1)
+ -------------------------------------------- -------------------------
(1 + d) Gamma(3 + r + 2d)

where b > -1/2 and d > -1/2

************************************************** *******
In the attachment I have provided the question with the all my working.

I am stuck at the integrations of the 3rd and 4th terms of my solution as I don't know how to integrate the incomplete beta function. Furthermore, their respective integrals I suspect should give 2 terms for a total of 4 for 'I' if you include the ones I have worked out. Yet in the question, 'I' is made up of 6 not 4 terms.

Also how does the hypergeometric function factor in the integration?

Finally is it possible there is an error in this question? I made sure that is was copied correctly.

Any hints, suitable resources or and help is greatly appreciated.

Thank you in advance

shawsend · #2 $Old$ November 6th, 2009, 06:00 PM

It would help if you wrote it clearly. You write:

$\mathop\int\int\limits_{\hspace{-10pt}x<y} F(x)^r F(y)^s F(x)(1-F(y))dxdy$

What are the actual limits on the variable x? Since it's symmetrical, it looks like:

$\int_{x_0}^{x_1}\int_{x_0}^y F(x)^r F(y)^s F(x)(1-F(y))dxdy$

Can you give the explicit values for $x_0,x_1$ or state what form the limits take?

I looked at your work and have some concerns about your new limits in the u-v coordinate system but for now, let's assume that part is ok. Is the following the integral you're trying to solve:

$\int_0^{1}\int_{0}^{v} u^{r+1}v^s(1-v) [ab u^{b-1}+cd(1-u)^{d-1}] [ab v^{b-1}+cd(1-v)^{d-1}] du dv$

If so, then what are the bounds on all the constants?

Hopefully then we can clearly see what integral you're trying to solve and can then proceed to work through it.

shawsend · #3 $Old$ November 6th, 2009, 07:05 PM

Here's what I believe are the integrals for others in case you're interested:

$I_1=\int_0^1 \frac{(ab)^2 v^{r+2b+s}(1-v)}{r+b+1}dv$

$I_2=\int_0^1 \frac{abcd v^{r+b+s+1}(1-v)^d}{r+b+1}dv$

$I_3=\int_0^{1} abcd v^{s+b-1}(1-v) \text{B}(v;r+2,d)dv$

$I_4=\int_0^{1} (cd)^2 v^s(1-v)^d \text{B}(v;r+2,d)dv$

where $\text{B}(v;s,t)$ is the incomplete Beta function

I went through the first one and apparently r,b s are integers since you're taking factorials and I get the same answer as you. You write further down in solving for the second one:

"However this is not the 5th term of I that I was asked to show. If you look closely, you will see that in the denominator in the Gamma function of my working there is a 3 instead of a 2. I just can't see where I'm going wrong."

I get a 3 also and so does Mathematica:

Code:

Integrate[((a*b*c*d)/(r + b + 1))*
   v^(r + b + s + 1)*(1 - v)^d, {v, 0, 1}]

(1/(1 + b + r))*(a*b*c*d*
   If[Re[d] > -1 && Re[b + r + s] > -2, 
    (Gamma[1 + d]*Gamma[2 + b + r + s])/
     Gamma[3 + b + d + r + s], 
    Integrate[(1 - v)^d*
      v^(1 + b + r + s), {v, 0, 1}, 
     Assumptions -> Re[d] <= -1 || 
       Re[b + r + s] <= -2]])

Therefore, I believe your initial answer is incorrect and the 2 should be changed to a 3. And as far as integrating the beta function, what about setting it up as a double integral and switching the order of integration?

$I_3=\int_0^{1} abcd v^{s+b-1}(1-v) \text{B}(v;r+2,d)dv=$ $abcd\int_0^1 v^{s+b-1}(1-v)\int_0^v u^{r+1}(1-u)^{d-1}dudv$

urcing · #4 $Old$ November 11th, 2009, 06:57 AM

Thank you so much for the you quick response, moreover for the proving to me that I wasn't making some catastrophic mistake in my calculation of the the second term. Thanks again!

As for you suggestion of switching the order of integration according the theorem on the evaluation of double integrals, the limits on the outside integral should be fixed points and not variable; so I don't think that would be the best approach. Can you thing of another way to tackle this?

I am curious. Is this possible to integrate the last two integrals with the incomplete beta function in the mathematica program? If yes, does it come out the the remaining terms I am looking for? And could you post the results and any insightful explanations.

Once again any help you provide was, is and will be GREATLY appreciated!

Thank you again.

p.s
I apologize that the it wasn't written clearly. How do you enter text so that it comes up exactly as the math symbols?

shawsend · #5 $Old$ November 11th, 2009, 07:44 AM

Hello Urcing. Thought you moved on. That was a beautiful problem. I was able to obtain what I believe is the third one (using limits of 0 and 1) by switching the order of integration:

$I_3=\int_0^{1} abcd v^{s+b-1}(1-v) \text{B}(v;r+2,d)dv$ $=abcd\int_0^1 v^{s+b-1}(1-v)\int_0^v u^{r+1}(1-u)^{d-1}dudv$

$=abcd\int_0^1\int_u^1 u^{r+1}(1-u)^{d-1}(v^{b+s-1}-v^{b+s})dvdu$

$=\frac{abcd}{(b+s)(b+s+1)}\frac{\Gamma(r+2)\Gamma(d)}{\Gamma(r+2+d)}$

$+\frac{abcd}{(b+s+1)}\frac{\Gamma(b+s+r+3)\Gamma(d)}{\Gamma(b+s+r+d+3)}$
$-\frac{abcd}{(b+s)}\frac{\Gamma(b+s+r+2)\Gamma(d)}{\Gamma(b+s+r+2+d)}$

And that is what you reported based on the other answers. The fourth one was quite a challenge via switching the order of integration and term-wise integration of the hypergeometric series (from the first integration) although I obtained a slightly different format for the final hypergeometric expression in the solution. I didn't save my work but as I recall, I was able to confirm (or get close to) Mathematica's solution:

$\int_0^1 (cd)^2 v^s(1-v)^d\text{B}(v;r+2,d)dv$ $=\frac{(cd)^2\Gamma(1+d)\Gamma(2+r)\Gamma(3+r+s)}{\Gamma(3+r)\Gamma(4+d+r+s)}\text{F}[\{1-d,2+r,3+r+s\},\{3+r,4+d+r+s\},1]$
with F being the generalized hypergeometric series.

(. . . said it was a challenge

)

To enter math text, you need to write it as LaTex code. Click on the math formulas and a window appears with the latex code or better yet, just do a quote to see the full latex code.

urcing · #6 $Old$ November 11th, 2009, 11:14 PM

I totally misunderstood your hint! After seeing you reply I did the calculation and got the next 3 terms of the solution. Thank you soooooooooo much!

I'm working on the final term now but I'm finding it difficult. So I that I don't make the same mistake twice, could you explain what you mean by "term-wise integration of the hypergeometric series" especially the part "from the first integration". I don't understand how the hypergeometric function comes in.

Thank you again

shawsend · #7 $Old$ November 12th, 2009, 06:31 AM

Quote:

Originally Posted by urcing $View Post$

I'm working on the final term now but I'm finding it difficult. So I that I don't make the same mistake twice, could you explain what you mean by "term-wise integration of the hypergeometric series" especially the part "from the first integration". I don't understand how the hypergeometric function comes in.

Thank you again

I've been re-doing it and have some concerns with the results. There is of course the issue of justifying switching order of integration, summation, and term-wise integration which may not be valid but for now, let's assume that it is just to see what we get and let's just drop the $(cd)^2$ term for now:

$\begin{aligned} I_4&=\int_0^{1} v^s(1-v)^d \text{B}(v;r+2,d)dv \\ &=\int_0^1 u^{r+1}(1-u)^{d-1}\int_u^1 v^s(1-v)^d dv du \end{aligned}$

Now use the Binomial theorem for the $(1-v)^d$ term:

$(1-v)^d=\sum_{k=0}^{\infty} \binom{d}{k} (-v)^k=\sum_{k=0}^{\infty}\frac{(-d)_k}{k!} v^k$

That last sum is a hypergeometric function right? So then:

$\begin{aligned} I_4&=\int_0^1 u^{r+1}(1-u)^{d-1} \int_u^1v^s \sum_{k=0}^{\infty}\frac{(-d)_k}{k!} v^k dv du \\ &=\int_0^1 u^{r+1}(1-u)^{d-1} \left(\sum_{k=0}^{\infty} \frac{(-d)_k}{k!} \frac{v^{s+k+1}}{s+k+1}\biggr|_u^1\right) du \\ &=\int_0^1 u^{r+1}(1-u)^{d-1}\left\{\sum_{k=0}^{\infty}\frac{(-d)_k}{k!} \frac{1}{s+k+1}\left(1-u^{s+k+1}\right)\right\}du \end{aligned}$

And at that point, we can then express the integral as an infinite sum of beta functions right?

$I_4=\left[\sum_{k=0}^{\infty} \frac{(-d)_k}{k!} \frac{\beta(r+2,d)}{s+k+1}-\sum_{k=0}^{\infty} \frac{(-d)_k}{k!} \frac{\beta(r+s+k+3,d)}{s+k+1}\right]$

I can't go further on my own at this point, but Mathematica evaluates that sum to a hypergeometric expression however, Mathematica cannot resolve the equality of that expression to the expression obtained by direct integration of the original beta expression:

Via direct integration by Mathematica:
$I_a=\int_0^{1} v^s(1-v)^d\beta(v;r+2,d)dv,\quad d>-1/2, s>-1, r+s>-3, d\notin\mathbb{Z}$
$=\frac{\Gamma(1+d)\Gamma(2+r)\Gamma(3+r+s)}{\Gamma(3+r)\Gamma(4+d+r+s)} \,F\left[\{1+d,2+r,3+r+s\},\{3+r,4+d+r+s\},1\right]$
and from the above analysis using Mathematica to evaluate the last infinite sum expression:
$I_b=\int_0^{1} v^s(1-v)^d\beta(v;r+2,d)dv$
$=\left[\sum_{k=0}^{\infty} \frac{(-d)_k}{k!} \frac{\beta(r+2,d)}{s+k+1}-\sum_{k=0}^{\infty} \frac{(-d)_k}{k!} \frac{\beta(r+s+k+3,d)}{s+k+1}\right]$
$=\frac{\beta(2+r,d)\Gamma(1+d)\Gamma(2+s)}{(1+s)\Gamma(2+d+s)}-\frac{\text{F}\left[\{-d,1+s,3+r+s\},\{2+s,3+d+r+s\},1\right]}{(1+s)\Gamma(2+d+s)}$

Now, let d=1/2 and take 2500 points with $0<s<5$ and $0<r<5$ and plot the difference $I_a-I_b$ which is shown below. Note the difference is of the order of $10^{-16}$ . That of course doesn't mean everything is ok but does give me some confidence that the method above is probably valid at least in a select range of values.

urcing · #8 $Old$ November 13th, 2009, 04:21 PM

When I computed it by hand I got this

$I_4=\frac{\Gamma(1+d)\Gamma(3+r+s)}{(2+r)\Gamma(4+d+r+s)} \,F\left[\{1-d,2+r,3+r+s\},\{3+r,4+d+r+s\},1\right]$

Differs from Mathematica by the 1+d agrument compared to my 1-d but everything else is the same. I'll post my working tomorrow night for you to see and look through.

Also is
$\int_u^{1} v^s(1-v)^ddv=\int_0^{1-u} v^d(1-v)^s\$

a legitimate manipulation? If yes, why?

shawsend · #9 $Old$ November 13th, 2009, 06:09 PM

Hello Urcing. Yeah, I just got it. Took a while for me. I was making it too complicated above. First I noted that:

$\beta(a,b)=\int_0^1 u^{a-1}(1-u)^{b-1} du=\int_0^1 u^{a-1}\sum_{k=0}^{\infty} \frac{(1-b)_k}{k!} u^k du$
$=\sum_{k=0}^{\infty} \frac{(1-b)_k}{k!(a+k)}=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$

via the generalized Binomial Theorem previously mentioned above. Then:

$\begin{aligned} I_4&=\int_0^{1} v^s(1-v)^d\beta(v;r+2,d)dv\\ &=\int_{0}^1 v^s (1-v)^d\sum_{k=0}^{\infty}\frac{(1-d)_k}{k!}\frac{v^{r+k+2}}{r+k+2}dv\\ &=\int_0^1 \sum_{k=0}^{\infty}\frac{(1-d)_k}{k!(r+k+2)}v^{r+s+k+2}(1-v)^d dv\\ &=\sum_{k=0}^{\infty} \frac{(1-d)_k}{k!(r+k+2)} \beta(r+s+k+3,d+1) \end{aligned}$
$\hspace{25pt}=\sum_{k=0}^{\infty} \frac{(1-d)_k}{k!(r+k+2)}\frac{\Gamma(r+s+k+3)\Gamma(d+1)}{\Gamma(r+s+k+d+4)}$

Note the first three:

$k=0:\frac{1}{r+2}\frac{\Gamma(r+s+3)\Gamma(d+1)}{\Gamma(r+s+d+4)}$
$k=1:\frac{(1-d)\Gamma(r+s+4)\Gamma(d+1)}{(r+3)\Gamma(r+s+d+5)}=\frac{(1-d)(3+r+s)\Gamma(3+r+s)\Gamma(d+1)}{(r+3)(4+r+s+d)\Gamma(4+r+s+d)}$
$k=2:\frac{(1-d)(1-d+1)\Gamma(r+s+5)\Gamma(d+1)}{(r+4)\Gamma(r+s+d+6)}$
$=\frac{(1-d)(1-d+1)(4+r+s)(3+r+s)\Gamma(3+r+s)\Gamma(d+1)}{(r+4)(5+r+s+d)(4+r+s+d)\Gamma(4+r+s+d)}$

Where we can begin to see the rising factorial terms emerging via the relationship $\Gamma(n+x)=\Gamma(1+(n-1)+x)=(n-1+x)\Gamma(n-1+x)$ . Mathematica is factoring out the term:

$\frac{1}{r+2}\frac{\Gamma(r+s+3)\Gamma(1+d)}{\Gamma(r+s+d+4)}$ and then adding the expression $\frac{(2+r)_k}{(3+r)_k}$ in the hypergeometric expression to account for the $(r+k)$ terms. Nice trick I think. I made a typo above. It should be $1-d$ :

$=\frac{\Gamma(1+d)\Gamma(2+r)\Gamma(3+r+s)}{\Gamma(3+r)\Gamma(4+d+r+s)} \,F\left[\{1-d,2+r,3+r+s\},\{3+r,4+d+r+s\},1\right]$

And the other integral is just a $w=1-v$ substitution right?

urcing · #10 $Old$ November 15th, 2009, 05:48 AM

I guess I won't need to put on the working anymore! lol.
The upper and lower limits seem ok. I thought so too but then wouldn't v = w+1? After reading

Beta function - Wikipedia, the free encyclopedia

I wonder if this equation exploits the fact that beta function is symmetric? What do you think?

Apparently this transformation was used to get the the expression for the 6th term in the solution but I'm not entriely comvinced it is a legitimate step. Finally, is the difference between our solution and the given solution zero.

shawsend · #11 $Old$ November 15th, 2009, 09:26 AM

Hi. I think:

$\int_u^{1} v^s(1-v)^ddv=\int_0^{1-u} v^d(1-v)^s dv$ is correct. Note if $w=1-v$ then $v=1-w$ .

Glad you asked about the difference between out solution and the solution you gave for the sixth term. They are not the same. Note if I just numerically calculate $\int_0^1 v^s(1-v^d)\beta(v,r+2,d)dv$ for $d=1/2, r=2, s=5$ , the numerical results agrees with our solution but not the one you initially posted if I'm interpreting your notation correctly. That is:

$F[-2,1+d,1+2d;2+d,r+3;1]= _3\hspace{-4pt}F_2[\{-2,1+d,1+2d\},\{2+d,r+3\},1]$

Then this is the results from Mathematica:

Code:

(* Numerical Results *)
In[38]:=
NIntegrate[v^s*(1 - v)^d*Beta[v, 2 + r, d] /. 
   {d -> 1/2, r -> 2, s -> 5}, {v, 0, 1}]
(* Our results *)
N[((Gamma[1 + d]*Gamma[3 + r + s])/
     ((2 + r)*Gamma[4 + d + r + s]))*HypergeometricPFQ[
     {1 - d, 2 + r, 3 + r + s}, {3 + r, 4 + d + r + s}, 1] /. 
   {d -> 1/2, r -> 2, s -> 5}]
(* Posted results *)
N[((Gamma[1 + 2*d]*Gamma[2 + r])/((1 + d)*Gamma[3 + r + 2*d]))*
    HypergeometricPFQ[{-2, 1 + d, 1 + 2*d}, {2 + d, r + 3}, 
     1]] /. {d -> 1/2, r -> 2, s -> 5}

Out[38]=
0.014092256949399784

Out[39]=
0.014092256949399806

Out[40]=
0.02019047619047619

I believe what we've calculated is correct and what you initially posted is not for the following reason: We can do a very good job of numerically integrating the original double integral:
$I=\mathop\int\int\limits_{\hspace{-15pt}x<y} F(x)^r F(y)^s F(x)\left[1-F(y)\right]dxdy$

Let's do so for the following list of parameter values:

$a=1/2, b=1/3, c=1/4, d=3/4, m=1, r=2,s=1/4$ given $x=F^{-1}(u)=m+au^b-c(1-u)^d$ :

We do so by first generating a list of values $(u,F^{-1}(u))$ in the range $0\leq u\leq 1$ , then calculating an interpolation function for the values. The Mathematica command Interpolation[values], does that below and returns the interpolation function myF:

Code:

g[a_, b_, c_, d_, m_, f_] := m + a f^b - c (1 - f)^d;
x0 = N[g[1/2, 1/3, 1/4, 3/4, 1, 0]]
x1 = N[g[1/2, 1/3, 1/4, 3/4, 1, 1]]

p1 = Plot[g[1/2, 1/3, 1/4, 3/4, 1, f], {f, 0, 1}]
lns = Flatten[
   Cases[Normal[First[p1]], Line[pts_] -> pts, {0, Infinity}], 1];
ipts = {#[[2]], #[[1]]} & /@ lns;
myF = Interpolation[ipts]

Now, the function $myF(x)$ in the range $(x_0,x_1)=(0.75,1.5)$ is a very good approximation to $F(x)$ . I've plotted below the function $f(x,y)=myF(x)^r myF(y)^s myF(x)\left[1-myF(y)\right]$ and note that the blue region above the line $x=y$ is the region of integration of the original integral. We can then calculate this numerically and obtain:

$I=\int_{x_0}^{x_1} \int_{x_0}^y myF(x)^r myF(y)^s myF(x)\left[1-myF(y)\right]dxdy\approx 0.00261072$

Note also, when we numerically integrate over the u-v domain we obtain good agreement with the original integral:

$I=\int_0^1 \int_0^v u^{r+1}v^{s}(1-v)\left[abu^{b-1}+cd(1-u)^{d-1}\right]\left[ab v^{b-1}+cd(1-v)^{d-1}\right]dudv$
$\approx 0.00261072$

And finally, when we assemble all the pieces of the above analysis, this numeric value agrees well with the results of that expression:

$I=\frac{(ab)^2}{(r+b+1)(r+s+2b+1)(r+s+2b+2)}$
$+\frac{abcd}{r+b+1}\frac{\Gamma(r+s+b+2)\Gamma(1+d)}{\Gamma(r+s+b+d+3)}$
$+\frac{abcd\Gamma(2+r)\Gamma(d)}{(b+s)(b+s+1)\Gamma(r+2+d)}$
$+\frac{abcd\Gamma(b+s+r+3)\Gamma(d)}{\Gamma(b+s+r+d+3)(b+s+2)}$
$-\frac{abcd\Gamma(b+s+r+2)\Gamma(d)}{(b+s)\Gamma(b+s+r+s+d)}$
$+\frac{(cd)^2\Gamma(1+d)\Gamma(2+r)\Gamma(3+r+s)}{\Gamma(3+r)\Gamma(4+d+r+s)} F[\{1-d,2+r,3+r+s\},\{3+r,4+d+r+s\},1]$
$\approx 0.00261072$

shawsend · #12 $Old$ November 15th, 2009, 01:22 PM

Urcing, this is the Mathematica code I used to numerically construct the function F(x) and compute the integral in F, the integral in u and v, and then the expression we derived. You may want to try and find a machine running Mathematica and run the code and change the parameters as you wish.

Code:

Clear[a, b, c, d, m, r, s, f, x0, x1]; 
a = 1/2; 
b = 3/2; 
c = 1/4; 
d = 3/4; 
m = 1; 
r = 3/2; 
s = 1/4; 

(* set up data to calculate F(x) *)

g[aval_, bval_, cval_, dval_, mval_, fval_] := 
   mval + aval*fval^bval - cval*(1 - fval)^dval; 
x0 = N[g[a, b, c, d, m, 0]]; 
x1 = N[g[a, b, c, d, m, 1]]; 
myPoints = Table[{g[a, b, c, d, m, f], f}, {f, 0, 1, 0.01}]; 
myInverse = Interpolation[myPoints]; 

(* now construct integrand with interpolated value of F(x) *)

fval[x_, y_] := myInverse[x]^r*myInverse[y]^s*myInverse[x]*
    (1 - myInverse[y]); 

(* integrate in x-y plane *)
NIntegrate[fval[x, y], {y, x0, x1}, {x, x0, y}]

(* integrate in u-v plane *)
NIntegrate[u^(r + 1)*v^s*(1 - v)*(a*b*u^(b - 1) + 
    c*d*(1 - u)^(d - 1))*(a*b*v^(b - 1) + c*d*(1 - v)^(d - 1)), 
  {v, 0, 1}, {u, 0, v}]

(* compute symbolic value *)
i1 = (a*b)^2/((r + b + 1)*(r + s + 2*b + 1)*(r + s + 2*b + 2)); 
i2 = ((a*b*c*d)/(r + b + 1))*(Gamma[r + s + b + 2]/
     Gamma[r + s + b + d + 3])*Gamma[1 + d]; 
i3 = a*b*c*d*Gamma[d]*(Gamma[2 + r]/((b + s)*(b + s + 1)*
       Gamma[r + 2 + d]) + Gamma[b + s + r + 3]/
      (Gamma[b + s + r + d + 3]*(b + s + 1)) - 
     Gamma[b + s + r + 2]/((b + s)*Gamma[b + s + r + 2 + d])); 
i4 = (c*d)^2*((Gamma[1 + d]*Gamma[2 + r]*Gamma[3 + r + s])/
     (Gamma[3 + r]*Gamma[4 + d + r + s]))*
    HypergeometricPFQ[{1 - d, 2 + r, 3 + r + s}, 
     {3 + r, 4 + d + r + s}, 1]; 
N[i1 + i2 + i3 + i4]