How do we compute an integral over a unit simplex?

Question

Let $ n \ge 2 $ and $ T > n $ be integers. The joint-distribution of eigenvalues in the Wishart ensemble subject to the underlying covariance matrix being equal to an identity matrix is given as follows:

\begin{eqnarray} &&f_{n,T}\left(\lambda_1,\cdots,\lambda_n\right) := \frac{(\sqrt{\pi})^{n-1} }{2^{(n T)/2} (\prod\limits_{j=0}^{n-1} \Gamma[\frac{T-j}{2}]) \cdot (\prod\limits_{j=0}^{n-3} \Gamma[\frac{n-j}{2}]) \cdot n! } \cdot \\ && \left(\prod\limits_{1 \le i < j \le n} \left| \lambda_i - \lambda_j \right| \right) \cdot \left(\prod\limits_{j=1}^n \lambda_j^{\frac{T-n-1}{2}} \right) \cdot e^{-\frac{1}{2} \sum\limits_{j=1}^n \lambda_j} \quad (i) \end{eqnarray}

The result above has been derived by using the Wishart distribution then transitioning from the matrix elements to the eigenvalues and the angles by using this link and then integrating over the angles.

Now, clearly the distribution in question has to be normalized to one. By integrating the above over ${\mathbb R}_+^n $ and decomposing the integral over a unit simplex $ \Delta _{n-1}:= \left\{ \lambda_i | 0 \le \lambda_1 \le \lambda_2 \le \cdots \le \lambda_{n-1} \le 1 \right\} $ we arrived at the following identity:

\begin{eqnarray} &&\int\limits_{\Delta_{n-1}} \prod\limits_{1 \le i < j \le n} \left| \lambda_i - \lambda_{i-1} - \lambda_j + \lambda_{j-1}\right| \cdot \prod\limits_{i=1}^n (\lambda_i -\lambda_{i-1})^{\frac{T-n-1}{2}} \cdot \prod\limits_{i=1}^{n-1} d\lambda_{i-1} = \\ && \frac{n! \prod\limits_{j=0}^{n-3} \Gamma[\frac{n-j}{2}] \prod\limits_{j=0}^{n-1} \Gamma[\frac{T-j}{2}] }{(\sqrt{\pi})^{n-1} \cdot \Gamma[\frac{n T}{2}]} \quad (ii) \end{eqnarray} subject to $\lambda_0 = 0 $ and $\lambda_n = 1$.

The Mathematica code below confirms the identity $(ii)$ numerically.

n = 3; T = 
 RandomInteger[{n + 1, 20}]; Clear[NN]; a1 =.; a2 =.; Clear[l];
NNangl[n_] := (Pi)^(n - 1) ( Sqrt[Pi])^Binomial[n - 1, 2]/
    Product[Gamma[(n - j)/2], {j, 0, n - 3}];
NN[T_, n_] := 
  1/n! NNangl[
    n] 1/(2^(n T/2) (Sqrt[Pi]^Binomial[n, 2] Product[
         Gamma[(T - j)/2], {j, 0, n - 1}]) );
NN[T, n] NIntegrate[
  Product[Abs[l[i] - l[j]], {i, 1, n}, {j, i + 1, n}] Product[
    l[j]^((T - n - 1)/2), {j, 1, n}]  Exp[-1/2 Sum[l[j], {j, 1, n}]], 
  Evaluate[Sequence @@ Table[{l[j], 0, Infinity}, {j, 1, n}]]]
(\[Pi]^(1/2 (-1 + n)) Gamma[(n T)/2])/(n! (\!\(
\*UnderoverscriptBox[\(\[Product]\), \(j = 0\), \(\(-3\) + n\)]\(Gamma[
\*FractionBox[\(1\), \(2\)]\ \((\(-j\) + n)\)]\)\)) \!\(
\*UnderoverscriptBox[\(\[Product]\), \(j = 0\), \(\(-1\) + n\)]\(Gamma[
\*FractionBox[\(1\), \(2\)]\ \((\(-j\) + T)\)]\)\))
  NIntegrate[
  Product[Abs[
     If[i == n, 1, l[i]] - If[i == 1, 0, l[i - 1]] - 
      If[j == n, 1, l[j]] + If[j == 1, 0, l[j - 1]]], {i, 1, n}, {j, 
     i + 1, 
     n}] Product[(If[j == n, 1, l[j]] - 
       If[j == 1, 0, l[j - 1]])^((T - n - 1)/2), {j, 1, n}], 
  Evaluate[Sequence @@ 
    Table[{l[j], If[j == 1, 0, l[j - 1]], 1}, {j, 1, n - 1}]]]

Now, my question is straightforward. How would you be proving identity $(ii)$ otherwise?

Answer 1

As a matter of fact this identity is straightforward to verify by doing the consecutive integrations over $\lambda_{j} \in (\lambda_{j-1},1) $ for $j=n-1,n-2,\cdots, 1$ and splitting the integration range accordingly so that all the expressions under the absolute values are always positive. It is just that calculations are tedious. However the key to the verification procedure is the following identity for the incomplete beta function:

\begin{eqnarray} B_x \left(a+n,b+m\right) = \frac{a^{(n)} b^{(m)}}{(a+b)^{(n+m)}} \cdot \left[ \right. \\ \left. B_x\left(a,b\right) + \right. \\ \left. \frac{x^a (1-x)^{b-1}}{a+b-1} \sum\limits_{k=1}^m \frac{(a+b-1)^{(k)}}{b^{(k)}} \cdot (1-x)^k + \right. \\ \left. - \frac{x^{a-1} (1-x)^{b+m}}{a+b+m-1} \cdot \frac{(a+b)^{(m)}}{b^{(m)}} \sum\limits_{k=1}^n \frac{(a+b+m-1)^{(k)}}{a^{(k)}} \cdot x^k \right. \\ \left. \right] \quad (1) \end{eqnarray} where $a\ge 0$, $b \ge 0$, $x \in {\mathbb R}$ and $n \in {\mathbb N}$,$m \in {\mathbb N}$. This formula is a generalization of the identity from here..

Making use of $(1)$ we will be doing symbolic integrations using the following code:

l1 =.; a =.; b =.; x =.;
subst1 := 
  Beta[z_, x + n_., x + m_.] :> 
   Pochhammer[x, n]/
     Pochhammer[2 x + m, 
      n] (Pochhammer[x, m]/
        Pochhammer[2 x, m] (Beta[z, x, x] + 
         z^(x) (1 - z)^(x - 1)/(2 x - 1) Sum[
           Pochhammer[2 x - 1, k]/Pochhammer[x, k] (1 - z)^k, {k, 1, 
            m}]) - z^(x - 1) (1 - z)^(x + m)/(2 x + m - 1) Sum[ 
        Pochhammer[2 x + m - 1, k]/Pochhammer[x, k] z^k, {k, 1, n}]);
Simplify[(Simplify[(Expand[(l1 - a) (l1 - b) (b - 
          a) l1^(x - 1) a^(x - 1) b^(x - 1)] /. 
      a^n1_ b^n2_ :> (1 - l1)^(n1 + n2 + 1) (Beta[1/2, n1 + 1, 
           n2 + 1] - Beta[l1/(1 - l1), n1 + 1, n2 + 1]))] /. subst1), 
 Assumptions -> ... && x > 0]

Let us set $n=3$.Here the integrand reads:

\begin{equation} {\mathfrak I}(\lambda_1,\lambda_2) = |-2 \lambda_1+\lambda_2| |\lambda_1+\lambda_2-1| |2 \lambda_2-\lambda_1-1| \cdot \prod\limits_{j=1}^3 (\lambda_j -\lambda_{j-1})^{x-1} \end{equation}

Before we start we set $\lambda:= \lambda_1$ and $x = (T-n+1)/2$.

Now we have to consider three cases.

Firstly we have $\lambda \in (1/2,1) $ and then we have either $\lambda_2 \in (\lambda_1,(\lambda_1+1)/2)$ or $\lambda_2 \in ((\lambda_1+1)/2,1)$. It turns out that the integrals over $\lambda_2$ over both intervals are the same and their sum reads:

\begin{eqnarray} \int\limits_{\lambda_1}^1 {\mathfrak I}(\lambda_1,\lambda_2) d\lambda_2= \frac{4^{1-x} (1-\lambda )^{2 x} \lambda ^{x+1}}{x^2+x}-\frac{2^{1-2 x} (1-\lambda )^{2 x} \lambda ^x}{x^2+x}+\frac{2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x-1}}{x+1}+\frac{9\ 2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x+1}}{x+1}-\frac{3\ 4^{-x} (1-\lambda )^{2 x} \lambda ^x}{x+1} \quad (2a) \end{eqnarray}

Secondly we have $\lambda \in (1/3,1/2) $ and then we have either $\lambda_2 \in (\lambda_1,1-\lambda_1)$ or $\lambda_2 \in (1-\lambda_1,(1+\lambda_1)/2)$ or $\lambda_2 \in ((1+\lambda_1)/2,2 \lambda_1)$ or $\lambda_2 \in (2 \lambda_1,1)$. Again, it turns out that the first and the fourth integral are equal and the second and the third are also equal but different from the previous one. Adding all four of them gives the following result:

\begin{eqnarray} \int\limits_{\lambda_1}^1 {\mathfrak I}(\lambda_1,\lambda_2) d\lambda_2= \frac{4^{1-x} (1-\lambda )^{2 x} \lambda ^{x+1}}{x^2+x}-\frac{8 (1-2 \lambda )^x \lambda ^{2 x+1}}{x^2+x}-\frac{2^{1-2 x} (1-\lambda )^{2 x} \lambda ^x}{x^2+x}+\frac{4 \left((1-2 \lambda ) \lambda ^2\right)^x}{x^2+x}+\frac{2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x-1}}{x+1}+\frac{9\ 2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x+1}}{x+1}-\frac{3\ 4^{-x} (1-\lambda )^{2 x} \lambda ^x}{x+1} \quad (2b) \end{eqnarray}

Thirdly we have $\lambda \in (0,1/3) $ and then we have either $\lambda_2 \in (\lambda_1,2 \lambda_1)$ or $\lambda_2 \in (2\lambda_1,(1+\lambda_1)/2)$ or $\lambda_2 \in ((1+\lambda_1)/2,1- \lambda_1)$ or $\lambda_2 \in (1 - \lambda_1,1)$. Again, it turns out that the first and the fourth integral are equal and the second and the third are also equal but different from the previous one. Adding all four of them gives the following result:

\begin{eqnarray} \int\limits_{\lambda_1}^1 {\mathfrak I}(\lambda_1,\lambda_2) d\lambda_2= \frac{2 (1-2 \lambda )^x \lambda ^{2 x}}{x^2+x}+\frac{4^{1-x} (1-\lambda )^{2 x} \lambda ^{x+1}}{x^2+x}-\frac{4 (1-2 \lambda )^x \lambda ^{2 x+1}}{x^2+x}-\frac{2^{1-2 x} (1-\lambda )^{2 x} \lambda ^x}{x^2+x}-\frac{4 \lambda \left((1-2 \lambda ) \lambda ^2\right)^x}{x^2+x}+\frac{2 \left((1-2 \lambda ) \lambda ^2\right)^x}{x^2+x}+\frac{2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x-1}}{x+1}+\frac{9\ 2^{-2 x-1} (1-\lambda )^{2 x} \lambda ^{x+1}}{x+1}-\frac{3\ 4^{-x} (1-\lambda )^{2 x} \lambda ^x}{x+1} \quad (2c) \end{eqnarray}

All what we need to do now is to integrate $(2c)$ over $\lambda_1\in(0,1/3)$ then integrate $(2b)$ over $\lambda_1\in(1/3,1/2)$ and then $(2a)$ over $\lambda_1\in(1/2,1)$ and add all the results together. But this is simple. All the terms in the right hand sides contain incomplete beta functions of the form $B_{1/3}(,)$, $B_{1/2}(,)$, $B_{1}(,)$ or $B_{2/3}(,)$, $B_{1}(,)$, $B_{2}(,)$. It is easy to write a Mathematica script that replaces the respective expressions $\lambda^{\cdots}(1-\lambda)^{\cdots}$ and $\lambda^{\cdots}(1-2\lambda)^{\cdots}$ by the respective incomplete beta functions. We have done this and miraculously all the incomplete beta functions disappeared leaving the complete beta functions only and then we obtained the following final result:

\begin{eqnarray} &&\int\limits_0^1 \int\limits_{\lambda_1}^1 {\mathfrak I}(\lambda_1,\lambda_2) d\lambda_2 d\lambda_1=\\ &&\frac{2^{-2 x-1} (x B_1(x,2 x+1)-2 (3 x+2) B_1(x+1,2 x+1)+9 x B_1(x+2,2 x+1)+8 B_1(x+2,2 x+1)+4 B_1(2 x+1,x+1)-4 B_1(2 x+2,x+1))}{x (x+1)}=\\ &&\frac{9\ 2^{1-2 x} \Gamma (2 x) \Gamma (x+2)}{\Gamma (3 x+4)}=\\ &&\frac{3 \Gamma \left(\frac{T-2}{2}\right) \Gamma \left(\frac{T-1}{2}\right) \Gamma \left(\frac{T}{2}\right)}{\sqrt{\pi } \Gamma \left(\frac{3 T}{2}\right)} \end{eqnarray}

This finishes the proof. As always we provide Mathematica code that verifies all the partial results:

(*Integrating over l2 \in (l1,1)*)
{x} = RandomReal[{1, 10}, 1];
l1 = RandomReal[{1/2, 1}];
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l2, l1, 1}]
(2^(-1 - 2 x) (1 - l1)^(2 x) l1^(-1 + x))/(1 + x) - (
 3 4^-x (1 - l1)^(2 x) l1^x)/(1 + x) + (
 9 2^(-1 - 2 x) (1 - l1)^(2 x) l1^(1 + x))/(1 + x) - (
 2^(1 - 2 x) (1 - l1)^(2 x) l1^x)/(x + x^2) + (
 4^(1 - x) (1 - l1)^(2 x) l1^(1 + x))/(x + x^2)
l1 = RandomReal[{1/3, 1/2}];
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l2, l1, 1}]
(2^(-1 - 2 x) (1 - l1)^(2 x) l1^(-1 + x))/(1 + x) - (
 3 4^-x (1 - l1)^(2 x) l1^x)/(1 + x) + (
 9 2^(-1 - 2 x) (1 - l1)^(2 x) l1^(1 + x))/(1 + x) - (
 2^(1 - 2 x) (1 - l1)^(2 x) l1^x)/(x + x^2) + (
 4^(1 - x) (1 - l1)^(2 x) l1^(1 + x))/(x + x^2) - (
 8 (1 - 2 l1)^x l1^(1 + 2 x))/(x + x^2) + (4 ((1 - 2 l1) l1^2)^x)/(
 x + x^2)
l1 = RandomReal[{0, 1/3}];
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l2, l1, 1}]
(2^(-1 - 2 x) (1 - l1)^(2 x) l1^(-1 + x))/(1 + x) - (
 3 4^-x (1 - l1)^(2 x) l1^x)/(1 + x) + (
 9 2^(-1 - 2 x) (1 - l1)^(2 x) l1^(1 + x))/(1 + x) - (
 2^(1 - 2 x) (1 - l1)^(2 x) l1^x)/(x + x^2) + (
 2 (1 - 2 l1)^x l1^(2 x))/(x + x^2) + (
 4^(1 - x) (1 - l1)^(2 x) l1^(1 + x))/(x + x^2) - (
 4 (1 - 2 l1)^x l1^(1 + 2 x))/(x + x^2) + (2 ((1 - 2 l1) l1^2)^x)/(
 x + x^2) - (4 l1 ((1 - 2 l1) l1^2)^x)/(x + x^2)
(Integrating over l1 \in (0,1) )
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l1, 0, 
  1/3}, {l2, l1, 1}]
(2^(-1 - 2 x) Beta[1/3, x, 1 + 2 x])/(1 + x) - (
 3 4^-x Beta[1/3, 1 + x, 1 + 2 x])/(1 + x) - (
 2^(1 - 2 x) Beta[1/3, 1 + x, 1 + 2 x])/(x + x^2) + (
 9 2^(-1 - 2 x) Beta[1/3, 2 + x, 1 + 2 x])/(1 + x) + (
 4^(1 - x) Beta[1/3, 2 + x, 1 + 2 x])/(x + x^2) + (
 2^(1 - 2 x) Beta[2/3, 1 + 2 x, 1 + x])/(x + x^2) - (
 2^(1 - 2 x) Beta[2/3, 2 + 2 x, 1 + x])/(x + x^2)
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l1, 1/3, 
  1/2}, {l2, l1, 1}]
1/(x (1 + x)) 2^(-1 - 
  2 x) (-x Beta[1/3, x, 1 + 2 x] + (4 + 6 x) Beta[1/3, 1 + x, 
     1 + 2 x] - 8 Beta[1/3, 2 + x, 1 + 2 x] - 
   9 x Beta[1/3, 2 + x, 1 + 2 x] + x Beta[1/2, x, 1 + 2 x] - 
   4 Beta[1/2, 1 + x, 1 + 2 x] - 6 x Beta[1/2, 1 + x, 1 + 2 x] + 
   8 Beta[1/2, 2 + x, 1 + 2 x] + 9 x Beta[1/2, 2 + x, 1 + 2 x] - 
   4 Beta[2/3, 1 + 2 x, 1 + x] + 4 Beta[2/3, 2 + 2 x, 1 + x] + 
   4 Beta[1, 1 + 2 x, 1 + x] - 4 Beta[1, 2 + 2 x, 1 + x])
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l1, 1/2, 
  1}, {l2, l1, 1}]
1/(x (1 + x)) 2^(-1 - 
  2 x) (-x Beta[1/2, x, 1 + 2 x] + (4 + 6 x) Beta[1/2, 1 + x, 
     1 + 2 x] - 8 Beta[1/2, 2 + x, 1 + 2 x] - 
   9 x Beta[1/2, 2 + x, 1 + 2 x] + x Beta[1, x, 1 + 2 x] - 
   4 Beta[1, 1 + x, 1 + 2 x] - 6 x Beta[1, 1 + x, 1 + 2 x] + 
   8 Beta[1, 2 + x, 1 + 2 x] + 9 x Beta[1, 2 + x, 1 + 2 x])
(The whole thing.)
NIntegrate[
 Abs[(2 l1 - l2) (l1 + l2 - 1)] Abs[-2 l2 + l1 + 
    1] l1^(x - 1) (l2 - l1)^(x - 1) (1 - l2)^(x - 1), {l1, 0, 1}, {l2,
   l1, 1}]
1/(x (1 + x)) 2^(-1 - 
  2 x) (x Beta[1, x, 1 + 2 x] - 2 (2 + 3 x) Beta[1, 1 + x, 1 + 2 x] + 
   8 Beta[1, 2 + x, 1 + 2 x] + 9 x Beta[1, 2 + x, 1 + 2 x] + 
   4 Beta[1, 1 + 2 x, 1 + x] - 4 Beta[1, 2 + 2 x, 1 + x])
(9 2^(1 - 2 x) Gamma[2 x] Gamma[2 + x])/Gamma[4 + 3 x]
(3 Gamma[1/2 (-2 + T)] Gamma[1/2 (-1 + T)] Gamma[T/2])/(
 Sqrt[[Pi]] Gamma[(3 T)/2]) /. T :> 2 x + 2