Dirichlet's integral $\int_{V}\ x^{p}\,y^{q}\,z^{r}\ \left(\, 1 - x - y - z\,\right)^{\,s}\,{\rm d}x\,{\rm d}y\,{\rm d}z$ I found such an exercise:
Calculate the Dirichlet's integral:
$$
\int_{V}\ x^{p}\,y^{q}\,z^{r}\
\left(\, 1 - x - y - z\,\right)^{\,s}\,{\rm d}x\,{\rm d}y\,{\rm d}z
\quad\mbox{where}\quad p, q, r, s >0
$$
and $V=\left\{\,\left(\, x,y,z\,\right) \in {\mathbb R}^{3}_{+}:
x + y + z\ \leq\ 1\right\}$

I thought that I could put $x + y + z = \alpha$. I got a clue, that it is a correct approach, but I should also put $y + z = \alpha\beta$ and $z=\alpha\beta\gamma$. So:
$z=\alpha\beta\gamma\,,\quad y=\alpha\beta\left(\, 1 - \gamma\,\right)\,,\quad x=\alpha\left(\, 1 - \beta\,\right)$
Should I change $x,y,z$ under the integral sign to $\alpha,\beta,\gamma$ now ?.
 A: For Type I Dirichlet integrals, one has the formula:
$$\int_{\Delta_n} f\left(\sum_{k=1}^n t_k\right) \prod_{k=1}^n t_k^{\alpha_k-1}\prod_{k=1}^n dt_k = \frac{\prod_{k=1}^n \Gamma(\alpha_k)}{\Gamma(\sum_{k=1}^n\alpha_k)}\int_0^1
f(\tau) \tau^{(\sum_{k=1}^n\alpha_k)-1} d\tau$$
where $$\Delta_n = \bigg\{ (x_1,\ldots,x_n) \in [0,\infty)^n :  \sum_{k=1}^n x_k \le 1 \bigg\}$$
is the standard $n$-simplex. For a proof of a very similar formula where $\Delta_n$ is replaced by $[0,\infty)^n$, see this answer.
It will show you how to carry out the computation in your original approach.
Apply it to your integral with 
$$f(w) = (1-w)^s\quad\text{ and }\quad
\begin{cases}
\alpha_1 = p + 1\\
\alpha_2 = q + 1\\
\alpha_3 = r + 1
\end{cases},
$$
one find
$$\begin{align}
  \int_{\Delta_3}(1-x-y-z)^s x^p y^q z^r dxdydz
= & \frac{\Gamma(p+1)\Gamma(q+1)\Gamma(r+1)}{\Gamma(p+q+r+3)}\int_0^1 (1-\tau)^s t^{p+q+r+2} d\tau\\
= &\frac{\Gamma(p+1)\Gamma(q+1)\Gamma(r+1)\Gamma(s+1)}{\Gamma(p+q+r+s+4)}
\end{align}
$$
A: Why not to map $(x,y,z)$ into $(u^2,v^2,w^2)$ and integrate over a spherical sector? 
With the first change of variables we have:
$$ I = 8\iiint_S u^{2p+1} v^{2q+1} w^{2r+1} (1-(u^2+v^2+w^2))^s\,d\mu $$
where $S=\{(u,v,w)\in(0,+\infty)^3: u^2+v^2+w^2\leq 1\}$. 
By setting $u=\rho\cos\theta\sin\phi, v=\rho\sin\theta\sin\phi, w=\rho\cos\phi$, we get:
$$ I = 8\int_{0}^{1}\iint_{(0,\pi/2)^2}\rho^{2p+2q+2r+5}(1-\rho^2)^s\cos^{2p+1}\theta\sin^{2q+1}\theta\cos^{2r+1}\phi\sin^{2p+2q+3}\phi\,d\mu\,d\rho,$$
but since, due to the properties of the Euler Beta function:
$$ \int_{0}^{1}\rho^{2p+2q+2r+5}(1-\rho^2)^s\,d\rho = \frac{\Gamma(3+p+q+r)\Gamma(1+s)}{2\Gamma(4+p+q+r+s)},$$
$$\int_{0}^{\pi/2}\cos^{2p+1}\theta\sin^{2q+1}\theta\,d\theta = \frac{\Gamma(p+1)\Gamma(q+1)}{2\Gamma(2+p+q)},$$
$$\int_{0}^{\pi/2}\cos^{2r+1}\phi\sin^{2p+2q+3}\phi\,d\phi=\frac{\Gamma(2+p+q)\Gamma(1+r)}{2\Gamma(3+p+q+r)},$$
it follows that:
$$ \color{red}{I=\frac{\Gamma(p+1)\Gamma(q+1)\Gamma(r+1)\Gamma(s+1)}{\Gamma(4+p+q+r+s)}}. $$
A: $\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\dsc}[1]{\displaystyle{\color{red}{#1}}}
 \newcommand{\expo}[1]{\,\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\Li}[1]{\,{\rm Li}_{#1}}
 \newcommand{\norm}[1]{\left\vert\left\vert\, #1\,\right\vert\right\vert}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
$\ds{\int_{V}\ x^{p}\,y^{q}\,z^{r}\
\pars{1 - x - y - z}^{\,s}\,\dd x\,\dd y\,\dd z\,;\
p, q, r, s >0}$ and $\ds{V=\braces{\pars{ x,y,z} \in {\mathbb R}^{3}_{+}:
x + y + z\ \leq\ 1}}$

\begin{align}&\color{#66f}{\large\left.
\int_{V}\ x^{p}\,y^{q}\,z^{r}\pars{1 - x - y - z}^{\, s}\,\dd x\,\dd y\,\dd z\,
\right\vert_{\, x + y + z\ < 1}}
\\[5mm]&=\left.\int_{0}^{\infty}\int_{0}^{\infty}
\int_{0}^{\infty}x^{p}\,y^{q}\,z^{r}
\pars{1 - x - y - z}^{\,s}\,\dd x\,\dd y\,\dd z\,
\right\vert_{\, x + y + z\ < 1}
\\[5mm]&=\int_{0}^{\infty}\int_{0}^{\infty}
\int_{0}^{\infty}x^{p}\,y^{q}\,z^{r}\
\Theta\pars{1 - x - y - z}\pars{1 - x - y - z}^{\,s}\,\dd x\,\dd y\,\dd z
\\[5mm]&=\int_{0}^{\infty}\int_{0}^{\infty}
\int_{0}^{\infty}x^{p}\,y^{q}\,z^{r}\ \overbrace{\int_{0^{-}}^{\infty}
\delta\pars{1 - x - y - z - \xi}\xi^{\,s}\,\dd\xi}
^{\dsc{\Theta\pars{1 - x - y - z}\pars{1 - x - y - z}^{\,s}}}\
\,\dd x\,\dd y\,\dd z
\\[5mm]&=\int_{0}^{\infty}\int_{0}^{\infty}
\int_{0}^{\infty}\int_{0^{-}}^{\infty}x^{p}\,y^{q}\,z^{r}\xi^{\,s}\ \overbrace{%
\int_{0^{+}\ -\ \infty\ic}^{0^{+}\ +\ \infty\ic}
\exp\pars{\tau\bracks{1 - x - y - z - \xi}}\,{\dd\tau \over 2\pi\ic}}
^{\dsc{\delta\pars{1 - x - y - z - \xi}}}\,\,\,\,
\,\dd x\,\dd y\,\dd z\,\dd\xi
\\[5mm]&=\int_{0^{+}\ -\ \infty\ic}^{0^{+}\ +\ \infty\ic}\expo{\tau}\
\overbrace{\int_{0}^{\infty}x^{p}\expo{-\tau x}\,\dd x}
^{\dsc{\tau^{-p - 1}\ \Gamma\pars{p + 1}}}\
\underbrace{\int_{0}^{\infty}y^{q}\expo{-\tau y}\,\dd y}
_{\dsc{\tau^{-q - 1}\ \Gamma\pars{q + 1}}}\
\overbrace{\int_{0}^{\infty}z^{r}\expo{-\tau z}\,\dd z}
^{\dsc{\tau^{-r - 1}\ \Gamma\pars{r + 1}}}\times \\[2mm] &
\underbrace{\int_{0}^{\infty}\xi^{s}\expo{-\tau\xi}\,\dd\xi}
_{\dsc{\tau^{-s - 1}\ \Gamma\pars{s + 1}}}\,\,\,\,
\,{\dd\tau \over 2\pi\ic}
\\[5mm]&=\Gamma\pars{p + 1}\Gamma\pars{q + 1}\Gamma\pars{r + 1}\Gamma\pars{s + 1}\times \\[2mm] &
\underbrace{\int_{0^{+}\ -\ \infty\ic}^{0^{+}\ +\ \infty\ic}
\frac{\expo{\tau}}{\tau^{p + q + r + s + 4}}\,\,\,{\dd\tau \over 2\pi\ic}}
_{\ds{\dsc{1 \over \pars{p + q + r + s + 3}!}\ =\
      \dsc{1 \over \Gamma\pars{p + q + r + s + 4}}}}
\end{align}

Finally,
\begin{align}
&\color{#66f}{\large\left.%
\int_{V}\ x^{p}\,y^{q}\,z^{r}\pars{1 - x - y - z}^{\,s}\,\dd x\,\dd y\,\dd z
\,\right\vert_{\, x + y + z\ <\ 1}}
\\[5mm]&=\color{#66f}{\large{\Gamma\pars{p + 1}\Gamma\pars{q + 1}\Gamma\pars{r + 1}\Gamma\pars{s + 1} \over \Gamma\pars{p + q + r + s + 4}}}
\end{align}
A: The answer and more is also given in Handbook of Applicable Mathematics, Vol IV : Analysis, § 10.2 Special functions, (ISBN 04712770405 or 0471101499, 1982, John Wiley and Sons, Ltd, ed. by Walter Ledermann and Steven Vajda). With n > o the gamma function
$$\Gamma(n)= \int_0^\infty e^{-x}\cdot x^{(n-1)}\,dx$$
and after integrating by parts: $\Gamma(n)=(−1)!\quad$
Thus the gamma function is a generalization of the factorial function to the case where n may take non-integral values. For m, n > 0 the beta function
$$\beta(m,n)= \int_0^1x^{m-1}\cdot(1-x)^{(n-1)}\,dx$$
Usually their relation is in the form of (10.2.8 in this book):
$$\beta(m,n)=\frac{\Gamma(m)\cdot\Gamma(n)}{\Gamma(m+n)}$$
Now consider the integral
$$\iiint x^p  y^q z^r (1-x-y-z)^s\,dx\,dy\,dz$$
taken throughout the interior of the tetrahedron formed bij the planes $x = 0,\; y = 0, \; z = 0, \; x + y + z = 1$. Define new variables $\alpha = x + y + z, \quad \beta = y + z, \quad \gamma =z\quad$ or equivalently
$x = \alpha(1-\beta),\quad y=\alpha\beta(1-\gamma),\quad z = \alpha\beta\gamma \,$. When x, y, z are all positive and x + y + z < 1,
$\alpha, \beta, \gamma$ all lie between 0 and 1. Hence the tetrahedron transforms into a cube. The Jacobian matrix for this transformation to be inserted [see (5.10.3) and Theorem 6.2.5 in the book]
$$\frac{d(x,y,z)}{d(\alpha,\beta,\gamma)}=\alpha^2\beta$$
Substituting:
$$\int_0^1\int_0^1\int_0^1 \alpha^{(p+q+r)} \cdot\beta^{(q+r)} \cdot(1-β)^p \cdot \gamma^r (1-\gamma)^q (1-\alpha)^s \cdot (\alpha^2β)\,d\alpha \,d\beta\,d\gamma)$$ =
$$\left[\int_0^1 \alpha^{(p+q+r+2)}(1-\alpha)^s\,d\alpha\right]\left[\int_0^1 \beta^{(q+r+2)} (1-\beta)^p\,d\beta\right]\left[\int_0^1 \gamma^r (1-\gamma)^q\,d\gamma\right]$$ =
$$\beta(p+q+r+3,s+1)\cdot\beta(q+r+2,p+1)\cdot\beta(r+1,q+1)$$ And by ref. (10.2.8):
$$\frac{\Gamma(p+q+r+3)\Gamma(s+1)Γ(q+r+2)\Gamma(p+1)\Gamma(r+1)\Gamma(q+1)}{\Gamma(p+q+r+4)\Gamma(p+q+r+3)Γ(q+r+2)}$$ = $$\frac{\Gamma(p+1)\Gamma(q+1)\Gamma(r+1)\Gamma(s+1)}{\Gamma(p+q+r+4)}$$
This result, due to Dirichlet, can easily be extended to the case of n variables. In general, the integral which can be evaluated over the same region by a similar method with $f$ any continuous function, is:
$$\iiint x^py^q z^r f(x+y+z)\,dx\,dy\,dz$$
