Marginal distribution of the difference of two elements of a Dirichlet distributed vector Let $(X_1, X_2, \ldots X_n) \sim \text{Dirichlet}(\alpha_1, \alpha_2, \ldots \alpha_n)$, what is the marginal distribution of $X_i - X_j$?
 A: Analysis
The thread at Construction of Dirichlet distribution with Gamma distribution shows that the Dirichlet distribution with parameters $(\alpha_1, \alpha_2, \ldots, \alpha_{n+1})$ arises as the distribution of the ratios
$$X_i=\frac{Y_i}{Y_1+Y_2+\cdots + Y_{n+1}},$$
$i=1, 2, \ldots, n$ where the $Y_j$ are independently distributed with Gamma$(\alpha_j)$ distributions.  This permits two simplifications, because (a) we may choose the order of the $\alpha_i,$ reducing the question to the difference $X_2 - X_1$ and (b) since the sum $Y_3 + \cdots +Y_{n+1}$ has a Gamma$(\alpha_3+\cdots+\alpha_{n+1})$ distribution, we only have to consider the case $n=2.$
Solution
As usual for computing marginal distributions, change the variables in the Dirichlet distribution from $(X_1,X_2)$ to $(X,X+Y)$ where $Y$ represents the difference $X_2-X_1.$  This transformation has unit Jacobian, so the integrand remains otherwise unchanged.  The (unnormalized) integrand is proportional to
$$X^{\alpha_1-1}  (X + Y)^{\alpha_2-1}  (1 - X - (X+Y))^{\alpha_3-1}$$
and we have to integrate out $X$ to find the marginal distribution of $Y$.
Since none of the three terms can be negative, the integral breaks into two parts: one from $x=0$ to $x=(1-y)/2$ when $y \ge 0$ and the other from $x=-y$ to $x=(1-y)/2$ when $y \lt 0.$  These are well-known integrals--Mathematica or Maple or even tables of integrals should provide answers for you.  The results are
$$f_Y(y) = \frac{\Gamma (\alpha_1+\alpha_2+\alpha_3)2^{-\alpha_1}  }{\Gamma (\alpha_2)} y^{\alpha_2-1} (1-y)^{\alpha_1+\alpha_3-1} \,
   _2\tilde{F}_1\left(\alpha_1,1-\alpha_2;\alpha_1+\alpha_3;\frac{y-1}{2 y}\right)$$
for $0 \le y \lt 1$ and
$$ \eqalign{f_Y(y) &=
\csc \left(\pi  \left(\alpha _1+\alpha _2\right)\right) (1-y)^{\alpha _3-2} \\ &\left(\frac{(y-1) \sin \left(\pi  \alpha
   _1\right) \Gamma \left(\alpha\right) (-y)^{\alpha _1+\alpha _2} 
}{y \Gamma \left(\alpha _3\right)} \,_2\tilde{F}_1\left(\alpha
   _1,1-\alpha _3;\alpha _1+\alpha _2;\frac{2 y}{y-1}\right)\\
-\frac{\pi  2^{-\alpha_1-\alpha_2+1} \left(\alpha-1\right) (1-y)^{\alpha_1+\alpha_2}}{\Gamma
   \left(\alpha _1\right) \Gamma \left(\alpha _2\right)} 
\,_2\tilde{F}_1\left(1-\alpha _2,2-\alpha;-\alpha _1-\alpha _2+2;\frac{2 y}{y-1}\right)\right)
}
$$
for $-1 \lt y \lt 0,$ where $\alpha=\alpha_1+\alpha_2+\alpha_3$ and $\, _2\tilde{F}_1$ is the regularized Hypergeometric function.
For integral values of $\alpha_1+\alpha_2$ and negative values of $y$ you also have to take a limit (because the cosecant blows up).  For small integral values of the parameters the function is algebraic (because the Hypergeometric functions that are involved reduce to polynomials).  You can compute these by elementary means if you wish.
Verification
Here are histograms of independent simulations of $X_2-X_1$ using 50,000 iterations each for three combinations of $(\alpha_1,\alpha_2, \alpha_3+\cdots+\alpha_{n+1}).$ On each is superimposed the graph of $f_Y.$  All show close agreement.

