One of the topics of interest in random matrix theory is the joint distribution of eigenvalues of Gaussian random matrices in the Wishart ensemble, see
Harnad, John (ed.), Random matrices, random processes and integrable systems, CRM Series in Mathematical Physics. Berlin: Springer (ISBN 978-1-4419-9513-1/hbk; 978-1-4419-9514-8/ebook). xviii, 524 p. (2011). ZBL1215.15002., for example.
Now, in order to obtain the density of eigenvalues there is a need to integrate over simplices with the integrand being a product of power functions and an exponential. This leads to a generalization of the Gamma function in $d $ dimensions and as such motivates the following definition below.
Let $d \ge 2$ be an integer. Let $\vec{n}:= (n_1,n_2,\cdots,n_{d-1}) \in {\mathbb N}_+^{d-1} $ where $n_j \ge 1 $ for $j=1,\cdots,d-1$. Finally let $\lambda_d \ge 0 $ be real.
Definition 1:
The multivariate lower, incomplete Gamma function is defined as follows:
\begin{equation} \gamma(\vec{n},\lambda_d) := \int\limits_{\Delta_{d-1}^{(-)}(\lambda_d)} \left( \prod\limits_{\xi=1}^{d-1} \lambda_\xi^{n_\xi-1} \right) \cdot e^{-\sum\limits_{\xi=1}^d \lambda_\xi} \cdot \prod\limits_{\xi=1}^{d-1} d\lambda_\xi \end{equation}
where $\Delta^{(-)}_{d-1} (\lambda) := \left\{ (\lambda_\xi)_{\xi=1}^{d-1} | 0 \le \lambda_1 \le \cdots \le \lambda_{d-1} \le \lambda\right\}$.
Result 1:
By using elementary integration and induction in the dimension $d$ and then differentiation with respect to parameters we have found the following formula below:
\begin{equation} \gamma(\vec{n},\lambda_d) = (-1)^{|\vec{n}|-d+1} \prod\limits_{\xi=1}^{d-1} \frac{d^{n_\xi-1}}{d a_\xi^{n_\xi-1}} \left.\left( \sum\limits_{j=1}^d (-1)^{d-j} \frac{e^{-(\sum\limits_{\xi=j}^d a_\xi) \cdot \lambda_d}}{\prod\limits_{\xi=1}^{j-1} (a_\xi+\cdots+a_{j-1}) \cdot \prod\limits_{\xi=j}^{d-1}(a_j+\cdots+ a_\xi)} \right)\right|_{(a_\xi)_{\xi=1}^d:>1} \tag{1a} \end{equation}
Definition 2:
The multivariate upper, incomplete Gamma function is defined as follows:
\begin{equation} \Gamma(\vec{n},\lambda_d) := \int\limits_{\Delta_{d-1}^{(+)}(\lambda_d)} \left( \prod\limits_{\xi=1}^{d-1} \lambda_\xi^{n_\xi-1} \right) \cdot e^{-\sum\limits_{\xi=1}^d \lambda_\xi} \cdot \prod\limits_{\xi=1}^{d-1} d\lambda_\xi \end{equation}
where $\Delta^{(+)}_{d-1} (\lambda) := \left\{ (\lambda_\xi)_{\xi=1}^{d-1} | \lambda \le \lambda_1 \le \cdots \le \lambda_{d-1} \le \infty\right\}$.
Result 2:
Using the same technique as above we have found the following formula below:
\begin{equation} \Gamma(\vec{n},\lambda_d) = (-1)^{|\vec{n}|-d+1} \prod\limits_{\xi=1}^{d-1} \frac{d^{n_\xi-1}}{d a_\xi^{n_\xi-1}} \left. \left( \frac{e^{-(\sum\limits_{\xi=1}^d a_\xi) \cdot \lambda_d}}{\prod\limits_{\xi=1}^{d-1} (a_\xi+\cdots+a_{d-1})} \right) \right|_{(a_\xi)_{\xi=1}^d:>1} \tag{2a} \end{equation}
The Mathematica code snippet below verifies the results above:
In[772]:= d = RandomInteger[{2, 6}];
L = RandomReal[{0, 5}];
n = RandomInteger[{1, 5}, d - 1];
(*Lower incomplete*)
NIntegrate[
Product[l[xi]^(
n[[xi]] - 1), {xi, 1, d - 1}] Exp[-Sum[l[xi], {xi, 1, d - 1}] - L],
Evaluate[
Sequence @@
Table[{l[j], If[j == 1, 0, l[j - 1]], L}, {j, 1, d - 1}]]]
D[Sum[(-1)^(Total[n] - j - 1)/1 Exp[-Sum[a[eta], {eta, j, d}] L]/
Product[Sum[a[eta], {eta, xi, j - 1}], {xi, 1, j - 1}] 1/
Product[Sum[a[eta], {eta, j, xi}], {xi, j + 0, d - 1}], {j, 1,
d}], Evaluate[
Sequence @@ Table[{a[xi], n[[xi]] - 1}, {xi, 1, d - 1}]]] /.
a[j_] :> 1
(*Upper incomplete*)
NIntegrate[
Product[l[xi]^(
n[[xi]] - 1), {xi, 1, d - 1}] Exp[-Sum[l[xi], {xi, 1, d - 1}] - L],
Evaluate[
Sequence @@
Table[{l[j], If[j == 1, L, l[j - 1]], Infinity}, {j, 1, d - 1}]]]
(-1)^(Total[n] - d + 1)
D[Exp[-Sum[a[eta], {eta, 1, d}] L]/1 1/
Product[Sum[a[eta], {eta, xi, d - 1}], {xi, 1 + 0, d - 1}],
Evaluate[
Sequence @@ Table[{a[xi], n[[xi]] - 1}, {xi, 1, d - 1}]]] /.
a[j_] :> 1
Out[775]= 0.0295689
Out[776]= 0.0295689
Out[777]= 0.00111378
Out[778]= 0.00111378
Now, having said all this the question would be twofold. Firstly, can we evaluate those derivatives in $(1a)$ and $(2a)$ in order to obtain some neat closed form expressions. Secondly, can we generalize those results to the case when the components of the vector $\vec{n}$ are real?
Here is an answer to the first question. Rather than computing those derivatives we used induction in the dimension to derive the following results below. Let us denote by $|\vec{n}| := n_1+\cdots+n_{d-1}$ the $L^1$-norm of our vector in question. Then we have:
\begin{equation} \Gamma(\vec{n}; \lambda) = e^{-d \cdot \lambda} \sum\limits_{j=0}^{|\vec{n}|-d+1} {\mathfrak A}^{(d-1)}_j(\vec{n}) \cdot \lambda^{|\vec{n}|-d+1-j} \tag{3} \end{equation}
where the coefficients in the sum above satisfy the following recurrence relations:
\begin{equation} {\mathfrak A}^{(d)}_j( \vec{n} ) = \sum\limits_{{\tilde j}=0}^j {\mathfrak A}^{(d-1)}_{\tilde j}(n_2,\cdots,n_d) \cdot (|\vec{n}| - d-{\tilde j})_{(j-{\tilde j})} \cdot \frac{1}{d^{1+j-{\tilde j}}} \tag{3a} \end{equation}
for $j=0,\cdots,|\vec{n}|-d$ subject to ${\mathfrak A}^{(1)}_j(n) := (n-1)_{(j)}$.
Let us the following notation $n_{p,q} := n_p+ n_{p+1} + \cdots + n_q$. Then we have:
\begin{equation} \gamma(\vec{n},\lambda) = \sum\limits_{p=0}^{d-1} e^{-p \cdot \lambda} \sum\limits_{j=0}^{n_{d-p,d-1}-p} {\mathcal A}^{(d-1)}_{p,j} (\vec{n}) \cdot \lambda^{n_{d-p,d-1}-p-j} \tag{4} \end{equation} where the coefficients in the sum satisfy the following recurrence relations:
\begin{eqnarray} &&{\mathcal A}^{(d)}_{p,j}(\vec{n}) = \\ &&\left\{ \begin{array}{lll} -\sum\limits_{{\tilde j}=0}^j {\mathcal A}^{(d-1)}_{p-1,{\tilde j}} (n_1,\cdots, n_{d-2}) \cdot (n_{d-p+1,d}-p-{\tilde j})_{(j-{\tilde j})} \cdot \frac{1}{p^{1+j-{\tilde j}}} & \quad \mbox{if $\quad p\ge 1$} \\ \sum\limits_{{\tilde p}=1}^d \sum\limits_{{\tilde j}=0}^{n_{d-{\tilde p}+1,d-1}-{\tilde p}+1} {\mathcal A}^{(d-1)}_{{\tilde p}-1,{\tilde j}}(n_1,\cdots,n_{d-2}) \cdot (n_{d-{\tilde p}+1,d} - {\tilde p}-{\tilde j})! \cdot \frac{1}{p^{n_{d-{\tilde p}+1,d} - {\tilde p}-{\tilde j}+1}} & \quad \mbox{if $\quad p=0$} \end{array} \tag{4a} \right. \end{eqnarray} for $p=0,\cdots,d$ and $j=0,\cdots,n_{d-p+1,d}-p$ subject to ${\mathcal A}^{(1)}_{p,j}(n):= (n-1)! \cdot 1_{p=0} - (n-1)_{(j)} \cdot 1_{p\ge 1}$.
The code below demonstrates how significant is the speed-up when computing the quantity in question from $(3)$ rather than from its definition.