Probability that a Linear Combination of Dirichlet Random Variables is a Distribution I've been putting a lot of thought on this problem, but it seems I ran out of ideas. Any help would be appreciated!  Suppose we generate two probability vectors $\boldsymbol{\theta}_1, \boldsymbol{\theta}_2 \sim \operatorname{IID Dirichlet}(\boldsymbol{\alpha}) \in \mathbb{R}^L$  where $\boldsymbol{\alpha} = (\alpha_1, \ldots, \alpha_L)$ with $\alpha_1 = \ldots = \alpha_L = 1$.
Now, suppose we take a scalar value $a \geq 1$ and we define the vector $\mathbf{X}$ as:
\begin{equation}
\mathbf{X} = a \boldsymbol{\theta}_1 + (1-a) \boldsymbol{\theta}_2.
\end{equation}
It is clear that $\sum X_i = 1$, so the norming requirement of probability is met.  In order for $\mathbf{X}$ to be a probability vector we also need $\mathbf{X} \geqslant \mathbf{0}$.  How can I compute the probability that this is true?

What I have tried: I've tried looking at Chen (2013) to solve this problem but the $z$ in there seems too restrictive.  (I may also not be using their results correctly.)  I have also tried using another trick using Hoeffding's inequality, but none of my methods have worked.  I'm also new on the ideas of concentration bounds, so that may have been a problem too.  Any help would be appreciated!
 A: To facilitate analysis, we will define the scalar value $\phi \equiv (a-1)/a$, which will serve as an alternative parameterisation of the problem.  Using this value you can rewrite the probability of interest as follows:
$$\begin{align}
\mathbb{P}(\mathbf{X} \geqslant \mathbf{0})
&= \mathbb{P}(a \boldsymbol{\theta}_1 + (1-a) \boldsymbol{\theta}_2 \geqslant \mathbf{0}) \\[6pt]
&= \mathbb{P}(a \boldsymbol{\theta}_1 \geqslant (a-1) \boldsymbol{\theta}_2) \\[6pt]
&= \mathbb{P} \Big( \boldsymbol{\theta}_1 \geqslant \frac{a-1}{a} \cdot \boldsymbol{\theta}_2 \Big) \\[6pt]
&= \mathbb{P}( \boldsymbol{\theta}_1 \geqslant \phi\cdot  \boldsymbol{\theta}_2), \\[6pt]
\end{align}$$
Since $0 \leqslant \phi < 1$, your question boils down to finding the probability that one Dirichlet random variable is at least as large (componentwise) as some positive proportion (less than one) of another identically distributed Dirichlet random variable.
So far as I am aware, there is no closed form expression for this probability.  It can be written as an integral over the joint density function of the random variables:
$$\begin{align}
\mathbb{P}(\mathbf{X} \geqslant \mathbf{0})
&= \mathbb{P}( \boldsymbol{\theta}_1 \geqslant \phi \cdot \boldsymbol{\theta}_2), \\[6pt]
&= \iint \limits_{\boldsymbol{\theta}_1 \geqslant \phi \cdot \boldsymbol{\theta}_2} \ \text{Dirichlet}(\boldsymbol{\theta}_1|\boldsymbol{\alpha}) \text{Dirichlet}(\boldsymbol{\theta}_2|\boldsymbol{\alpha}) \ d\boldsymbol{\theta}_1 d \boldsymbol{\theta}_2. \\[6pt]
&= \frac{1}{\text{B}(\boldsymbol{\alpha})^2} \iint \limits_{\boldsymbol{\theta}_1 \geqslant \phi\cdot  \boldsymbol{\theta}_2} \prod_{i=1}^L (\theta_{1,i} \theta_{2,i})^{\alpha_i-1} \ d\boldsymbol{\theta}_1 d \boldsymbol{\theta}_2. \\[6pt]
\end{align}$$
In practice you would need to compute this integral using numerical methods, or approximate it using Monte Carlo methods or importance sampling.  So long as $L$ is not too large, it should be possible to compute this probability to a reasonable level of accuracy using Monte Carlo sampling.
