Sampling a Dirichlet Distribution efficiently

Question

I'm sampling a DirichletDistribution like so:

n = 5000;
par1 = Table[1./n,n];
data1 = RandomVariate[DirichletDistribution[par1],50];  
(* This takes about 1.5 seconds *)
par2 = Table[1.,n];
data2 = RandomVariate[DirichletDistribution[par2],50];  
(* This takes about 0.10 seconds *)

If I play with the value of the parameters, I find that giving values > 1 radically improves the performance. Any clues on how to circumvent this?

EDIT : The reason why I'm comparing par1 and par2 is that I was hoping for a shortcut following this kind of reasoning (although I know that this is false): DirichletDistribution[n*par] / n

Answer 1

Dirichelet random variates are constructed from gamma random variates. Let $$ v_j \stackrel{\text{iid}}{\sim} \textsf{Gamma}(\alpha_j,1) $$ for $j = 1, \ldots, n$ and set $$ x_j = \frac{v_j}{\sum_{j=1}^n v_j}. $$ Then $x \sim \textsf{Dirichlet}(\alpha)$, where $x = (x_1,\ldots,x_n)$ and $\alpha = (\alpha_1,\ldots,\alpha_n)$.

Gamma random variates themselves are constructed using a complicated rejection method, the efficiency of which depends on the magnitude of $\alpha_j$. That is what you are seeing. As far as I know, there is no way around it.

As an aside, let me point out that Mathematica actually returns $(x_1,\ldots,x_{n-1})$ instead of $(x_1, \ldots, x_n)$. This is consistent with the rest of the usage for the Dirichlet distribution, but for some purposes it is unfortunate. For example, I use draws from the Dirichlet distribution to represent random discrete probability distributions and I require all $n$ components. When I try to compute $x_n = 1-\sum_{j=1}^{n-1} x_j$ for what Mathematica returns, I often get a small negative number (on the order of $10^{-16}$), which causes my code to blow up. As a consequence I have to generate Dirichlet draws from scratch.