I've already seen https://math.stackexchange.com/questions/1064995/marginal-of-dirichlet-distribution-is-beta-integral, but need to extend this to the $K$-variate case.
We have $\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_K\end{bmatrix}$ following a Dirichlet distribution with parameter vector $\mathbf{a} = \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_K\end{bmatrix}$ such that $\sum_{k=1}^{K}x_k = 1$ and $x_k \in [0, 1]$ for $k = 1, \dots, K$ with density function $$p(\mathbf{x}) = \dfrac{\Gamma(\sum_{i=1}^{K}a_i)}{\prod_{i=1}^{K}\Gamma(a_i)}x_1^{a_1-1}\cdots x_{K-1}^{a_{K-1}-1}\left(1-\sum_{\ell = 1}^{K-1}x_\ell\right)^{a_{K}-1}\text{.}$$ Given $j \in \{1, \dots, K-1\}$, we have $$p(x_j)=x_j^{a_j - 1}\underbrace{\int_{0}^{1}\cdots\int_{0}^{1}}_{K-2\text{ times }}\left(\prod_{p \neq j}x_p^{a_p - 1}\right)\left(1-\sum_{\ell = 1}^{K-1}x_\ell\right)^{a_{K}-1}\text{ d}\mathbf{x}_{-j}\tag{1}$$ where $\mathbf{x}_{-j}$ is $\mathbf{x}$ without $x_j$.
How does one evaluate the integral given in $(1)$?
Edit: I know that the integral is wrong, since it's not integrating over the simplex. But I'm not sure how the limits would be formed. No reference I've found has shown how to integrate this to find the marginals. Please DO NOT use the method found at http://www.mas.ncl.ac.uk/~nmf16/teaching/mas3301/week6.pdf; I've already seen this. I would like to see this done using integration.
Also, given that $x_K$ doesn't appear anywhere in the density above... I'm concerned that the way in which the $x_i$ are ordered matters. How would one obtain the density for $x_K$ in this case since $x_K$ doesn't explicitly appear in $p(\mathbf{x})$? Would one just find the PDF of $1-\sum_{\ell=1}^{K-1}x_{\ell}$?