What happens when merging random variables in Dirichlet distribution?

Question

Imagine that

$$ X_1,\dots,X_k \sim \mathrm{Dirichlet}(\alpha_1,\dots,\alpha_k) $$

Since $x_i \in (0,1)$ for all $x_i$ and $\sum_{i=1}^k x_i = 1$, then $x_i$'s follow the first two axioms of probability and Dirichlet can be (and is) used as "distribution over distributions". Intuitively it should follow that

$$ X_1,\dots,X_{k-2},X_{k-1}+X_k \sim \mathrm{Dirichlet}(\alpha_1,\dots,\alpha_{k-2}, \alpha_{k-1}+\alpha_k) $$

since the properties of $x_i$'s would not change and the total "mass" of $\alpha_i$'s would not change.

But it's probability density function is

$$ f(x_1,\dots,x_k) \propto \prod_{i=1}^k x_i^{\alpha_i - 1}$$

and

$$ x_{k-1}^{\alpha_{k-1} - 1} \times x_k^{\alpha_k - 1} \ne (x_{k-1} + x_k)^{\alpha_{k-1} + \alpha_k - 1} $$

So merging of random variables in Dirichlet distribution does not seem lead to Dirichlet distribution over $k-1$ variables. What does it lead to?

Answer 1

It is a Dirichlet distribution having the expected parameters.

To see this, note that the vector-valued random variable $\mathbf{X}=(X_1, X_2, \ldots, X_k)$ has the same distribution as the variable

$$\frac{1}{\sum_i^k Y_i}\left(Y_1, Y_2, \ldots, Y_k\right)$$

where $Y_i \sim \Gamma(\alpha_i)$ are independently Gamma distributed. Write $Y_i^\prime=Y_i$ for $i=1, 2, \ldots, k-2$ and $Y_{k-1}^\prime = Y_{k-1}+Y_k$. The sum of all the $Y_i$ equals the sum of all the $Y_i^\prime$ and the distribution of $Y_{k-1}^\prime=Y_{k-1}+Y_k$ is $\Gamma(\alpha_{k-1}+ \alpha_k)$. Thus

$$X_{k-1} + X_k = \frac{1}{\sum_i^k Y_i} Y_{k-1} + \frac{1}{\sum_i^k Y_i} Y_{k} = \frac{1}{\sum_i^{k-1} Y_i^\prime} Y_{k-1}^\prime$$

and, for $i < k-1$,

$$X_i = \frac{1}{\sum_i^k Y_i} Y_{k-1} = \frac{1}{\sum_i^{k-1} Y_i^\prime} Y_{k-1}^\prime.$$

Therefore $\mathbf{X}^\prime=(X_1, X_2, \ldots, X_{k-2}, X_{k-1}+X_k)$ has the same distribution as

$$\frac{1}{\sum_i^{k-1} Y_i^\prime}\left(Y_1^\prime, Y_2^\prime, \ldots, Y_k^\prime\right).$$

This demonstrates that $\mathbf{X}^\prime$ has a Dirichlet$(\alpha_1, \alpha_2, \ldots, \alpha_{k-2}, \alpha_{k-1}+\alpha_k)$ distribution, QED.

---

The fault in the argument in the question lies in confusing the arithmetic sum of values $x_{k-1}+x_k$ with the sum of random variables $X_{k-1}+X_k$. The latter is performed with a convolution, of course.

Answer 2

Like I mentioned in the comments, all I am trying suggests that what you are suggesting actually does work.

As you mentioned, intuitively it makes sense if this would work, if $X_i$ represents the posterior draw for some probability $p_i$ for event $i$ happening, you should indeed be able to sum multiple $X_i$ to get the probability of multiple events $i$, if there is no possibility of both events happening at the same time. Since this is a multinomial setting, this is not the case, so we're good.

---

So let's show my simulations:

library('gtools')

K <- 10 
alpha <- c(rpois(K, 50)) #randomly generated alphas, just cause
k <- 2 # the number of alphas we are summing together

sim <- rdirichlet(10000, alpha)
plot(density(rowSums(sim[, 1:k]))) # the density of the summed variable

lines(density(rdirichlet(10000, c(sum(alpha[1:k]), alpha[-(1:k)]))[,1]), col = 'blue') 
# the density of the variable drawn from the Dirichlet distribution with summed alphas

Let's start with $\alpha = \{10, 10, 10\}$. Summing $\alpha_1$ and $\alpha_2$ should get us $Dir(2, \{20, 10\})$:

These marginal densities look pretty similar. According to wikipedia and this random lecture I found on the internet (through wikipedia), the marginal distribution of $X_i$ to the Dirichlet distribution is as follows:

$$X_i = Beta(\alpha_i, \sum_{k=1}^K\left[\alpha_k\right] - \alpha_i)$$

This relies on actually the same principle: summing all the $\alpha$ that are not $\alpha_i$ together, turning the corresponding multinomial to a binomial distributions with the outcomes $i$ and $not\textrm{-}i$. And indeed, if we fit the marginal distribution we would expect from the sum over the density in the previous picture, we see that it looks the same:

So theoretically, we should be able to take a Dirichlet distribution with a high $K$, sum all but one together and end up with a Beta distribution. Heck, let's try:

(99 times $\alpha \sim Pois(50)$ and 1 $\alpha = 1000$, summing together the random $\alpha\mathrm{s}$.)

To show that it also works for the joint densities, this is an example with $K=4$ and $\alpha = \{10,10,10,10\}$:

And this is with $K=3$ and $\alpha=\{20,10,10\}$:

$$\ddot{\smile}$$

---

So where does the confusion about $x_1^{\alpha_1-1}x_2^{\alpha_2-1}\neq (x_1 + x_2)^{\alpha_1+\alpha_2 -1 }$ come from?

Because when we sum $x_1$ and $x_2$ together, we don't just care about the density at $x_1$ and $x_2$, but for any combination of two $x$ that sums to $x_1+x_2$. I am not that strong in integration, so I'll not try and burn myself with that, but I really suggest reading this (page 3-4) for more information.

EDIT:

As @whuber correctly remarked, here is an example with low alphas, $K=4$, summing the first two $X_1$ and $X_2$: