Definite integral over a simplex Let $T^d:=\{(x_1,..,x_d):x_i \geq 0, \sum_{i=1}^{d}x_i \leq 1\}$ be the standard simplex in $\mathbb{R}^d$. Compute the integral $$\int_{T^d} x_1^{\nu_1-1}x_2^{\nu_2-1}...x_d^{\nu_d-1}(1-x_1-...-x_d)^{\nu_0-1}$$ where $\nu_i>0$.
Remark: I know the answer is $$\frac{\prod_{i=0}^{d}\Gamma(\nu_i)}{\Gamma(\sum_{i=0}^{d}\nu_i)}.$$ I evaluated for the case $d=2$ by using the transformation $(p-1)\iiint\limits_{T^{3}} x^{m-1}y^{n-1}z^{p-2} \mathrm{d}z\mathrm{d}y\mathrm{d}x= 
\iint\limits_{T^{2}} x^{m-1}y^{n-1}(1-x-y)^{p-1}\mathrm{d}y\mathrm{d}x$
and the substitutions $ \left\{\begin{matrix}x=u^2& &\\y=v^2& &\\z=w^2& &\end{matrix}\right.$and    $ \left\{\begin{matrix}u=r\sin\varphi\cos\theta& &\\v=r\sin\varphi\sin\theta& &\\w=r\cos\varphi& &\end{matrix}\right.,$ but this method is complex for computing the general case.
 A: Another, elegant way to perform your integration is explained by Jaynes in §18.10 of "Probability Theory: The Logic of Science". It uses a Laplace transform:
Slightly modifying your notation, your integral can be rewritten as
$$I(r) := \int_{[0,\infty)^n} x_1^{a_1} \dotsm x_n^{a_n} \,\mathrm{\delta}(r-x_1-\dotsb-x_n)\, \mathrm{d}x_1\dotsm\mathrm{d}x_n\qquad\text{with $r=1$}.$$
Note that the delta allows us to extend the integration to $\infty$, and renders the integration bounds independent from one another.
Now take the Laplace transform of $I(r)$ and eliminate the delta from it:
$$\int_0^{\infty} \exp(-zr) I(r)\,\mathrm{d}r =
\int_{[0,\infty)^n} x_1^{a_1} \dotsm x_n^{a_n} \,\exp[-z(x_1+\dotsb+x_n)]\, \mathrm{d}x_1\dotsm\mathrm{d}x_n 
\\= \prod_{i=1}^n\int_0^\infty x_i^{a_i}\,\exp(-z x_i)\,\mathrm{d}x_i
=\prod_{i=1}^n\frac{a_i!}{z^{a_i+1}}.$$
Now take the inverse Laplace transform of the last expression to obtain back $I(r)$:
$$I(r) = \frac{1}{2\pi\mathrm{i}}\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}
\exp(rz)\,\prod_{i=1}^n\frac{a_i!}{z^{a_i+1}}\,\mathrm{d}z
\\=\frac{\prod_{i=1}^n a_i!}{(n+\sum_ia_i-1)!} r^{n+\sum_ia_i-1}$$
and evaluate it at $r=1$. The path of integration of the last integral "passes
to the right of the origin, and is closed by an infinite semicircle over the left half-plane" (Jaynes).
Hope this helps.
A: If one wishes to use a change of variables, an appropriate transformation would be 
$(x_1,x_2,...,x_d)\to(y_1,y_2,...,y_d)$ such that
$\qquad y_1=x_1+x_2+\cdots+x_{d-1}+x_d$
$\qquad y_1y_2=x_1+x_2+\cdots+x_{d-2}+x_{d-1}$
$\qquad y_1y_2y_3=x_1+x_2+\cdots+x_{d-3}+x_{d-2}$
$\qquad \vdots$
$\qquad \prod_{i=1}^{d-1} y_i=x_1+x_2$
$\qquad \prod_{i=1}^{d} y_i=x_1$
$\implies x_1=y_1y_2y_3\,...\,y_{d-1}y_d$
$\qquad x_2=y_1y_2y_3\,...\,y_{d-1}(1-y_d)$
$\qquad x_3=y_1y_2y_3\,...\,y_{d-2}(1-y_{d-1})$
$\qquad \vdots$
$\qquad x_{d-1}=y_1y_2(1-y_3)$
$\qquad x_d=y_1(1-y_2)$ $\qquad\qquad$, so that $0\le y_i\le 1$ for all $i=1,2,...,d$.
The Jacobian of the transformation can be calculated as $\det(J)=-\prod_{i=1}^{d-1} y_i^{d-i}$
The integral thus reduces to
$$\displaystyle I=\int_0^1 y_1^{({\sum_{i=1}^{d}\nu_i-1})}(1-y_1)^{\nu_0-1}\mathrm{d}y_1\int_0^1 y_2^{({\sum_{i=1}^{d-1}\nu_i-1})}(1-y_2)^{\nu_d-1}\mathrm{d}y_2\cdots\int_0^1 y_{d-1}^{(\nu_1+\nu_2-1)}(1-y_{d-1})^{\nu_3-1}\mathrm{d}y_{d-1}\int_0^1 y_{d}^{(\nu_1-1)}(1-y_{d})^{\nu_2-1}\mathrm{d}y_{d}$$
$$=B\left(\sum_{i=1}^d\nu_i,\nu_0\right)B\left(\sum_{i=1}^{d-1}\nu_i,\nu_d\right)\cdots B\left(\nu_1+\nu_2,\nu_3\right)B\left(\nu_1,\nu_2\right)$$
$\qquad\qquad\qquad\quad=\displaystyle\frac{\prod_{i=0}^d\Gamma(\nu_i)}{\Gamma\left(\sum_{i=0}^d\nu_i\right)}$
A: Put $(n_0,n_1,\ldots, n_d)=:n$, $\sum_{k=0}^d n_k=:|n|$ and for $\lambda\geq0$ define
$$Q(n,\lambda):=\int\nolimits_{\lambda T^d}\prod_{1\leq k\leq d} x_k^{n_k-1} \ (\lambda -x_1-\ldots -x_d)^{n_0-1}\ {\rm d}(x)\ .$$
The substitution $x:=\lambda y\ (y\in T^d)$ gives
$$Q(n,\lambda)=\lambda^{d+|n|-(d+1)}\int\nolimits_{T^d}\prod_{1\leq k\leq d} y_k^{n_k-1} \ (1 -y_1-\ldots -y_d)^{n_0-1}\ {\rm d}(y)=\lambda^{|n|-1} Q(n,1)\ .$$
Letting the "outer integration" be  with respect to the last variable we now obtain
$$\eqalign{Q(n,1)&=\int_0^1 x^{n_d-1}\int\nolimits_{(1-x)T^{d-1}}\prod_{1\leq k\leq d-1} x_k^{n_k-1} (1-x_1-\ldots-x_{d-1}-x)^{n_0-1}\ {\rm d}(x')\ dx\cr
&= \int_0^1 x^{n_d-1} Q(n',1-x)\ dx=Q(n',1)\int_0^1 x^{n_d-1}(1-x)^{|n'|-1}\ dx\ .  \cr}$$
Here the last integral evaluates to $B(|n'|,n_d)$, so that we end up with the recursion
$$Q(n,1)=Q(n',1) B(|n'|,n_d)\ .$$
From here it should not be too difficult to arrive at the desired result.
A: Another way is the following: we can rewrite the integral as
$$
I_n(a)=\int_{\mathbb R_+^n} x_1^{\nu_1-1}\cdots x_n^{\nu_n-1} \delta(x_1+\cdots+x_n-a) dx_1 \cdots dx_n
$$
with $n=d+1$ and $\nu_0\equiv \nu_n$. 
By rescaling each $x_i$ by $a$, we see that
$$
I_n(a)=a^{\nu_1+\cdots + \nu_n-1} I_n(1).
$$
Then, 
$$\begin{aligned}
I_n(1)
&=\int_0^1 dx\, x^{\nu_n-1} I_{n-1}(1-x) \\
&= I_{n-1}(1)\int_0^1 x^{\nu_n-1} (1-x)^{\nu_1+\cdots+\nu_{n-1}-1}dx\\
&=I_{n-1}(1) B(\nu_n,\nu_1+\cdots+\nu_{n-1}).
\end{aligned}$$
Proceeding by induction,
$$
I_n(1)=I_1(1) B(\nu_2, \nu_1) B(\nu_3,\nu_1+\nu_2)\cdots B(\nu_{n-1},\nu_1+\cdots + \nu_{n-2}) B(\nu_n, \nu_1+\cdots+\nu_{n-1})
$$
and after a few cancellations between numerator and denominator,
$$
I_n(1)=\frac{\Gamma(\nu_1)\cdots \Gamma(\nu_n)}{\Gamma(\nu_1+\cdots+\nu_n)}\,.
$$
A: For any $a \in \mathbb{R}-\{0\}$  and $m,n \in (0,+\infty)$ one has
\begin{align*}
\frac{1}{a^{m+n-1}}\int_{0}^{a}y^{m-1}(a-y)^{n-1}\mathrm{d}y&=\frac{1}{a^{m+n-1}}\int_{0}^{a}a^{m+n-2}\left(\frac{y}{a}\right)^{m-1}\left(1-\frac{y}{a}\right)^{n-1}\mathrm{d}y\\&=\frac{1}{a}\int_{0}^{a}\left(\frac{y}{a}\right)^{m-1}\left(1-\frac{y}{a}\right)^{n-1}\mathrm{d}y\\&=\int_{0}^{1}x^{m-1}(1-x)^{n-1}\mathrm{d}x
\end{align*}
Thus $\int_{0}^{a}y^{m-1}(a-y)^{n-1}\mathrm{d}y=a^{m+n-1}B(m,n)$.
With above observation we can integrate $x_1^{\nu_1-1}\cdots x_d^{\nu_d-1}(1-x_1-...-x_d)^{\nu_{0}-1}$ over $T^{d}$ by integrating out variables one at each step
\begin{align*}
&\mathrel{\phantom{=}} \int_{T^{d}}x_1^{\nu_1-1}...x_d^{\nu_d-1}(1-x_1-...-x_d)^{\nu_{0}-1}\mathrm{d}\boldsymbol x\\
&=\int_{0}^{1}\int_{0}^{1-x_1}...\int_{0}^{1-x_1-...-x_{d-1}}x_1^{\nu_{1}-1}x_2^{\nu_{2}-1}...x_d^{\nu_{d}-1}(1-x_1-...-x_d)^{\nu_{0}-1}\mathrm{d}x_d\cdots\mathrm{d}x_2\mathrm{d}x_1\\
&=\int_{0}^{1}x_1^{\nu_{1}-1}\int_{0}^{1-x_1}x_2^{\nu_{2}-1}...\int_{0}^{1-x_1-...-x_{d-1}}x_d^{\nu_{d}-1}(1-x_1-...-x_d)^{\nu_{0}-1}\mathrm{d}x_d\cdots\mathrm{d}x_2\mathrm{d}x_1\\&=B(\nu_0,\nu_d)\int_{0}^{1}x_1^{\nu_{1}-1}\int_{0}^{1-x_1}x_2^{\nu_{2}-1}...\int_{0}^{1-x_1-...-x_{d-2}}x_{d-1}^{\nu_{d-1}-1}(1-x_1-...-x_{d-1})^{\nu_{0}+\nu_{d}-1}\mathrm{d}x_{d-1}\cdots\mathrm{d}x_2\mathrm{d}x_1\\&=B(\nu_{d-1},\nu_{0}+\nu_{d})B(\nu_0,\nu_d)\int_{0}^{1}x_1^{\nu_{1}-1}...\int_{0}^{1-x_1-...-x_{d-2}}x_{d-2}^{\nu_{d-2}-1}(1-x_1-...-x_{d-2})^{\nu_{0}+\nu_{d-1}+\nu_{d}-1}\,\mathrm{d}x_{d-2}\cdots\mathrm{d}x_1\\&=\cdots\\&=B(\nu_1,\nu_0+\nu_d+\nu_{d-1}+\cdots+\nu_{2})B(\nu_2,\nu_0+\nu_d+\nu_{d-1}+\cdots+\nu_{3})\cdots B(\nu_{d-1},\nu_0+\nu_d)B(\nu_{0},\nu_d)\\&=\frac{\Gamma(\nu_{0})\Gamma(\nu_{1})\cdots\Gamma(\nu_{d})}{\Gamma(\nu_{0}+\nu_{1}+\cdots+\nu_{d})}\\&=\frac{\Gamma(\boldsymbol {\nu})}{\Gamma(|\boldsymbol {\nu}|)}
\end{align*}
A: Here is the solution form my professor at TU Darmstadt: We start by realizing
\begin{align*}
\int_{[0,\infty]^{n}}\prod_{i=1}^{n}\exp(-x_{i})x_{i}^{\alpha_{i}-1}\mathrm{d}\mathbf{x} & =\prod_{i=1}^{n}\int_{0}^{\infty}\exp(-x_{i})x_{i}^{\alpha_{i}-1}\mathrm{d}x_{i}=\prod_{i=1}^{n}\Gamma(\alpha_{i}).
\end{align*}
We simplify
\begin{align*}
 & \exp\underbrace{\left(-\sum_{i=1}^{n}x_{i}\right)}_{=-t}\prod_{i=1}^{n}x_{i}^{\alpha_{i}-1}\\
 & =\exp\left(-t\right)\prod_{i=1}^{n-1}\left(tp_{i}\right)^{\alpha_{i}-1}\left(t\left(1-\sum_{j=1}^{n-1}p_{j}\right)\right)^{\alpha_{n}-1}\\
 & =\left(\exp\left(-t\right)t^{\left(\sum_{i=1}^{n}\alpha_{i}\right)-n}\right)\left(\prod_{i=1}^{n-1}p_{i}^{\alpha_{i}-1}\left(1-\sum_{j=1}^{n-1}p_{j}\right)^{\alpha_{n}-1}\right)
\end{align*}
using the transformations
\begin{align*}
t & =\sum_{j=1}^{n}x_{j},\\
x_{i} & =tp_{i},\\
x_{n} & =t\left(1-\sum_{i=1}^{n-1}p_{i}\right),
\end{align*}
which yields a Jacobian
\begin{align*}
\mathbf{J} & =\begin{bmatrix}t\mathbf{I} & \mathbf{p}\\
-t\mathbf{1}_{n-1}^{\intercal} & \left(1-\mathbf{1}_{n-1}^{\intercal}\mathbf{p}\right)
\end{bmatrix}.
\end{align*}
The blockdiagonal structure of the upper $(n-1)\times(n-1)$ matrix
makes the determinant 
\begin{align*}
\det\mathbf{J}=\left(\left(1-\mathbf{1}_{n-1}^{\intercal}\mathbf{p}\right)-\left(-t\mathbf{1}_{n-1}^{\intercal}\right)\left(t\mathbf{I}\right)^{-1}\mathbf{p}\right)\det\left(t\mathbf{I}\right)=t^{n-1}.
\end{align*}
Now, we use the substitution formula
\begin{align*}
 & \int_{[0,\infty]^{n}}\prod_{i=1}^{n}\exp(-x_{i})x_{i}^{\alpha_{i}-1}\mathrm{d}\mathbf{x}\\
 & =\int_{0}^{\infty}\left(\exp\left(-t\right)t^{\left(\sum_{i=1}^{n}\alpha_{i}\right)-n}\right)\int_{[0,1]^{n-1}}\prod_{i=1}^{n-1}p_{i}^{\alpha_{i}-1}\left(1-\sum_{j=1}^{n-1}p_{j}\right)^{\alpha_{n}-1}t^{n-1}\mathrm{d}\mathbf{p}~\mathrm{d}t\\
 & =\int_{0}^{\infty}\left(\exp\left(-t\right)t^{\left(\sum_{i=1}^{n}\alpha_{i}\right)-1}\right)\mathrm{d}t\int_{[0,1]^{n-1}}\prod_{i=1}^{n-1}p_{i}^{\alpha_{i}-1}\left(1-\sum_{j=1}^{n-1}p_{j}\right)^{\alpha_{n}-1}\mathrm{d}\mathbf{p}\\
 & =\Gamma\left(\sum_{i=1}^{n}\alpha_{i}\right)\int_{[0,1]^{n-1}}\prod_{i=1}^{n-1}p_{i}^{\alpha_{i}-1}\left(1-\sum_{j=1}^{n-1}p_{j}\right)^{\alpha_{n}-1}\mathrm{d}\mathbf{p}
\end{align*}
and, thus, we have
\begin{align*}
\int_{[0,1]^{n-1}}\prod_{i=1}^{n-1}p_{i}^{\alpha_{i}-1}\left(1-\sum_{j=1}^{n-1}p_{j}\right)^{\alpha_{n}-1}\mathrm{d}\mathbf{p}=\frac{\prod_{i=1}^{n}\Gamma(\alpha_{i})}{\Gamma\left(\sum_{i=1}^{n}\alpha_{i}\right)}.
\end{align*}
