# Sampling a Dirichlet Distribution efficiently

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

• Possibly related: (75303) – dr.blochwave May 4 '16 at 17:59
• Also, how are you timing this? With AbsoluteTiming[]? – dr.blochwave May 4 '16 at 18:00
• But while the Dirichlet distributions from par1 and par2 have the same means, they have different variances and are therefore different distributions. Why would you choose one over the other based on timing? (If they were the same distribution, then I would understand the concern about timing.) – JimB May 4 '16 at 18:39
• Also, both of these only take 0.1 seconds on my machine... – dr.blochwave May 4 '16 at 19:15
• I can confirm OP's problem, the first one takes 1.25 for me and the second 0.12. Mathematica 10.3, OS X. – C. E. May 4 '16 at 22:13

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.