Find marginal distribution of $K$-variate Dirichlet

Question

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}$?

Answer 1

The marginal distribution of $x_j$ is,

$$ p(x_j) = \frac{1}{B({\bf a})} \int_0^{1 - x_j} \int_0^{1 - x_j - x_1} \cdots \int_0^{1 - \sum_{k =1}^{K-2} x_k} \prod_{p=1}^{K-1} x_p^{a_p - 1} \left( 1 - \sum_{l=1}^{K-1} x_l \right)^{a_K - 1} d x_{K-1} d x_{K-2} \dots d x_1, $$ where $\bf a$ is the vector of all $a_j$ values, $B({\bf a})$ is the multivariate Beta function, and the integration variables do not include $d x_j.$ We can marginalize out $x_{K-1}$ by doing the innermost integral,

$$ \int_0^{1 - \sum_{k =1}^{K-2} x_k} x_{K-1}^{a_{K-1} -1} \left( 1 - \sum_{l=1}^{K-1} x_l \right)^{a_K - 1} d x_{K-1}. $$ Let $z (1 - \sum_{k=1}^{K-2} x_k) = x_{K-1}.$ Then the above integral becomes,

$$ \begin{split} & \left( 1 - \sum_{k=1}^{K-2} x_k \right)^{a_{K-1}-1} \int_0^1 z^{a_{K-1}-1} \left( [1 - z] \left[1 - \sum_{k=1}^{K-2} x_k \right] \right)^{a_K -1} \left( 1 - \sum_{k=1}^{K-2} x_k \right) dz \\ = & \left( 1 - \sum_{k=1}^{K-2} x_k \right)^{a_K + a_{K-1} -1} \int_0^1 z^{a_{K-1} -1} (1-z)^{a_K -1} dz \\ = & \left( 1 - \sum_{k=1}^{K-2} x_k \right)^{a_K + a_{K-1} -1} B(a_{K-1}, a_K). \end{split} $$

Plugging this into $p(x_j)$ we have,

$$ p(x_j) = \frac{B(a_{K-1}, a_K)}{B({\bf a})} \int_0^{1 - x_j} \int_0^{1 - x_j - x_1} \cdots \int_0^{1 - \sum_{k =1}^{K-3} x_k} \prod_{p=1}^{K-2} x_p^{a_p - 1} \left( 1 - \sum_{k=1}^{K-2} x_k \right)^{a_K + a_{K-1} -1} d x_{K-2} d x_{K-3} \dots d x_1. $$

Compare this to the original expression. This is very similar, except that the "determined" value $x_K$ has its role replaced by $x_K + x_{K-1}.$ We can now marginalize out $x_{K-2}$ purely by analogy to how we marginalized $x_{K-1}$ (i.e. replacing $a_K$ with $a_K + a_{K-1},$ etc.):

$$ p(x_j) = \frac{B(a_{K-1}, a_K) B(a_{K-2}, a_{K-1} + a_K)}{B({\bf a})} \int_0^{1 - x_j} \int_0^{1 - x_j - x_1} \cdots \int_0^{1 - \sum_{k =1}^{K-4} x_k} \prod_{p=1}^{K-3} x_p^{a_p - 1} \left( 1 - \sum_{k=1}^{K-3} x_k \right)^{a_K + a_{K-1} + a_{K-2} -1} d x_{K-3} d x_{K-3} \dots d x_1. $$

Note, however, that $B(a_{K-1},a_{K}) B(a_{K-2}, a_{K-1} + a_K) = B(a_{K-2}, a_{K-1}, a_K).$

As we can see, each iteration of this procedure involves taking out the last factor from the product in the integral, the last term from the sum in the integral, adding the last $a$ coefficient to the exponent on the sum, and adding the same coefficient to the list of variables in the multivariate Beta coefficient outside the integral. We need only apply this pattern to all variables in the integral, from the inside out. We get,

$$ p(x_j) = \frac{B({\bf a}_{-j})}{B({\bf a})} x_j^{a_j -1} (1 - x_j)^{\sum_{i \ne j} a_i -1}, $$ where ${\bf a}_{-j}$ is all values of $a_k$ for $k \ne j.$ Note that $\frac{B({\bf a}_{-j})}{B({\bf a})}$ is just $\frac{1}{B(a_j, \sum_{i \ne j} a_i)}.$ Therefore,

$$ p(x_j) = \text{Beta}(x_j; a_j, \sum_{i \ne j} a_i). $$

This was just a general version of the link provided by @marmle. (I even stole the idea of the integration substitution from it.)

EDIT: It's not clear from the notation, but in all of the integrals above, the integration variables do not include $x_j.$

Answer 2

It is possible to derive this integral using the expression for the normalization constant of the dirichlet distribution. The key idea is to recast the integral containing the terms not involving $x_j$ by performing a change of variables into an integral which is similar to the one used to derive the normalization constant of the dirichlet distribution.

I will perform an equivalent integration process by using the following density function.

$$p(\mathbf{x}) = \frac{\Gamma(a_0)}{\prod_{i=1}^{K}\Gamma(a_i)} \prod_{j=1}^{K} x_j^{a_j-1} \\ where \quad a_0 = \sum_{i=1}^{K}a_i \\ and \quad \sum_{i=1}^{K}x_i = 1 \\ Also, \quad x_K = 1 - \sum_{i=1}^{K-1}x_i$$

I aim to find the marginal distribution of $p(x_j)$.Now, in order to incorporate the constraint $\sum_{i=1}^{K}x_i = 1$, I make the observation that $\sum_{i=1}^{j-1}x_i + \sum_{l=j+1}^{K}x_l = 1 - x_j$ and derive the limits of integration as shown in the following equation. Additionally, integration is not performed over the variable $x_K$ which is set equal to $1 - \sum_{i=1}^{K-1}x_i$.

$$p(x_j) = x_j^{a_j-1}\int_{0}^{1-\sum_{p=1}^{K-2}x_p}\int_{0}^{1-\sum_{p=1}^{K-3}x_p}\ldots \int_{0}^{1-\sum_{p=1}^{j+1-1}x_p} \int_{0}^{1-x_j-\sum_{p=1}^{j-2}x_p} \ldots \int_{0}^{1-x_j-x_1} \int_{0}^{1-x_j}\frac{\Gamma(a_0)}{\prod_{i=1}^{K}\Gamma(a_i)} \prod_{q=1,\\ q\neq j}^{K} x_q^{a_q-1} dx_1 dx_2 \ldots dx_{j-1} dx_{j+1}\ldots dx_{K-2} dx_{K-1}$$

Now I perform a change of variables to arrive at an integral which is similar to the one used to calculate the normalization constant of the dirichlet distribution. More specifically, I use the new variables as shown below to perform this change of variables.

$$x_1' = \frac{x_1}{1-x_j}, ~~ x_2' = \frac{x_2}{1-x_j},~~\ldots~~,x_{j-1}' = \frac{x_{j-1}}{1-x_j},~~\ldots~~,x_K' = \frac{x_K}{1-x_j}$$

The limits of integration for the variable $x_1$ now get modified from $x_1 = \{0, x_j-1\}$ to $x_1' = \{0, 1\}$. For the variable $x_2$, these get transformed from $x_2 = \{0, 1-x_j-x_1\}$ to $x_2' = \{0, 1-x_1'\}$. A similar transformation occurs for the other variables. The resulting variable transformation results in a factor of $(1-x_j)^{K-2}$ as a result of the Jacobian term. This arises due to the presence of $K-2$ differential terms $(dx_i)$ in the integral and the observation that $dx_l=(1 - x_j)dx_l'$. The final integral can be simplified in the following manner.

$$p(x_j) = x_j^{a_j-1} (1-x_j)^{K-2}\int_{0}^{1-\sum_{p=1,p \neq j}^{K-2}x_p'}\int_{0}^{1-\sum_{p=1,p \neq j}^{K-3}x_p'}\ldots \int_{0}^{1-\sum_{p=1,p \neq j}^{j+1-1}x_p'} \int_{0}^{1-\sum_{p=1}^{j-2}x_p'} \ldots \int_{0}^{1-x_1'} \int_{0}^{1}\frac{\Gamma(a_0)}{\prod_{i=1}^{K}\Gamma(a_i)} \prod_{q=1,\\ q\neq j}^{K} {x_q'}^{a_q-1} (1-x_j)^{a_q-1} dx_1' dx_2' \ldots dx_{j-1}' dx_{j+1}'\ldots dx_{K-2}' dx_{K-1}'$$

It can be observed that $\prod_{q=1,\\ q\neq j}^{K}(1-x_j)^{a_q-1} = (1-x_j)^{\sum_{q=1,\\ q\neq j}^{K}a_q-1}=(1-x_j)^{a_0-a_j-(K-1)}$. This term can be taken outside the integration and multiplied with $(1-x_j)^{K-2}$ to yield $(1-x_j)^{a_0-a_j-1}$. The simplified integral is shown as follows.

$$p(x_j) = \frac{\Gamma(a_0)}{\prod_{i=1}^{K}\Gamma(a_i)} x_j^{a_j-1} (1-x_j)^{a_0-a_j-1}\int_{0}^{1-\sum_{p=1,p \neq j}^{K-2}x_p'}\int_{0}^{1-\sum_{p=1,p \neq j}^{K-3}x_p'}\ldots \int_{0}^{1-\sum_{p=1,p \neq j}^{j+1-1}x_p'} \int_{0}^{1-\sum_{p=1}^{j-2}x_p'} \ldots \int_{0}^{1-x_1'} \int_{0}^{1} \prod_{q=1,\\ q\neq j}^{K} {x_q'}^{a_q-1} dx_1' dx_2' \ldots dx_{j-1}' dx_{j+1}'\ldots dx_{K-2}' dx_{K-1}'$$

It can be observed that $\sum_{i=1,i \neq j}^{K}x_i' = \frac{\sum_{i=1,i\neq j}^{K} x_i}{1-x_j} = \frac{1-x_j}{1-x_j} = 1$. Hence, the resulting integral with the new variables is the same as the integral used for calculating the normalization constant of the dirichlet distribution except that the $x_j^{a_j-1}$ term is absent. Hence, the integral simplifies to $\frac{\prod_{i=1,i \neq j}^{K}\Gamma(a_i)}{\Gamma(a_0 - a_j)}$. After cancellation with the original gamma function terms present outside the integral, we obtain the final result shown below.

$$p(x_j) = \frac{\Gamma(a_0)}{\Gamma(a_j)\Gamma(a_0 - a_j)} x_j^{a_j-1} (1-x_j)^{a_0-a_j-1} \\ p(x_j) = \frac{1}{B(a_j, a_0-a_j)} x_j^{a_j-1} (1-x_j)^{a_0-a_j-1} \\ p(x_j) = Beta(x_j;a_j, a_0-a_j)$$