Editing Dirichlet distribution (section)

===From gamma distribution===
With a source of Gamma-distributed random variates, one can easily sample a random vector <math>x=(x_1, \ldots, x_K)</math> from the {{mvar|K}}-dimensional Dirichlet distribution with parameters <math>(\alpha_1, \ldots, \alpha_K)</math> . First, draw {{mvar|K}} independent random samples <math>y_1, \ldots, y_K</math> from [[Gamma distribution]]s each with density

<math display=block>\operatorname{Gamma}(\alpha_i, 1) = \frac{y_i^{\alpha_i-1} \; e^{-y_i}}{\Gamma (\alpha_i)}, \!</math>

and then set

<math display=block>x_i = \frac{y_i}{\sum_{j=1}^K y_j}.</math>

{{hidden begin|style=width:60%|ta1=center|border=1px #aaa solid|title=[Proof]}}
The joint distribution of the independently sampled gamma variates, <math>\{y_{i}\}</math>, is given by the product:

<math display=block>e^{-\sum_{i}y_{i}} \prod _{i=1}^{K} \frac{y_{i}^{\alpha _{i}-1}}{\Gamma (\alpha _{i})} </math>

Next, one uses a change of variables, parametrising <math> \{y_{i}\}</math> in terms of <math> y_{1}, y_{2}, \ldots , y_{K-1} </math> and <math> \sum _{i=1}^{K}y_{i}</math> , and performs a change of variables from  <math> y \to x </math> such that <math>\bar x = \textstyle\sum_{i=1}^{K}y_{i}, x_{1} = \frac{y_{1}}{\bar x}, x_{2} = \frac{y_{2}}{\bar x}, \ldots , x_{K-1} = \frac{y_{K-1}}{\bar x}</math>. Each of the variables <math>0 \leq x_{1}, x_{2}, \ldots , x_{k-1} \leq 1 </math> and likewise <math>0 \leq \textstyle\sum _{i=1}^{K-1}x_{i} \leq 1 </math>. One must then use the change of variables formula, <math> P(x) = P(y(x))\bigg|\frac{\partial y}{\partial x}\bigg| </math> in which <math>\bigg|\frac{\partial y}{\partial x}\bigg|</math> is the transformation Jacobian. Writing y explicitly as a function of x, one obtains 
<math>y_{1} = \bar xx_{1}, y_{2} = \bar xx_{2} \ldots   y_{K-1} = \bar xx_{K-1}, y_{K} = \bar x(1-\textstyle\sum_{i=1}^{K-1}x_{i})  </math>
The Jacobian now looks like
<math display=block>\begin{vmatrix}\bar x & 0 & \ldots &   x_{1} \\  0 & \bar x & \ldots & x_{2} \\ \vdots & \vdots & \ddots & \vdots \\ -\bar x & -\bar x & \ldots & 1-\sum_{i=1}^{K-1}x_{i} \end{vmatrix}</math>

The determinant can be evaluated by noting that it remains unchanged if multiples of a row are added to another row, and adding each of the first K-1 rows to the bottom row to obtain

<math display=block>\begin{vmatrix}\bar x & 0 & \ldots &   x_{1} \\  0 & \bar x & \ldots & x_{2} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & 1 \end{vmatrix} </math>

which can be expanded about the bottom row to obtain the determinant value <math>\bar x^{K-1}</math>. Substituting for x in the joint pdf and including the Jacobian determinant, one obtains:

<math display=block>
\begin{align}
&\frac{\left[\prod _{i=1}^{K-1}(\bar xx_{i})^{\alpha _{i}-1} \right] \left[\bar x(1-\sum_{i=1}^{K-1}x_{i})\right]^{\alpha_{K}-1}}{\prod _{i=1}^{K}\Gamma (\alpha _{i})}\bar x^{K-1}e^{-\bar x} \\
=&\frac{\Gamma(\bar\alpha)\left[\prod _{i=1}^{K-1}(x_{i})^{\alpha _{i}-1} \right] \left[1-\sum_{i=1}^{K-1}x_{i}\right]^{\alpha_{K}-1}}{\prod _{i=1}^{K}\Gamma (\alpha _{i})}\times\frac{\bar x^{\bar\alpha-1}e^{-\bar x}}{\Gamma(\bar\alpha)}
\end{align}
</math>
where <math>\bar\alpha=\textstyle\sum_{i=1}^K\alpha_i</math>. The right-hand side can be recognized as the product of a Dirichlet pdf for the <math>x_i</math> and a gamma pdf for <math>\bar x</math>. The product form shows the Dirichlet and gamma variables are independent, so the latter can be integrated out by simply omitting it, to obtain:
<math display=block>x_{1}, x_{2}, \ldots, x_{K-1} \sim \frac{(1-\sum_{i=1}^{K-1}x_{i})^{\alpha _{K}-1}\prod _{i=1}^{K-1}x_{i}^{\alpha _{i} -1}}{B(\boldsymbol{\alpha})} </math>

Which is equivalent to

<math display=block>\frac{\prod _{i=1}^{K} x_{i}^{\alpha_{i}-1}}{B(\boldsymbol{\alpha})} </math> with support <math> \sum_{i=1}^{K}x_{i}=1 </math>

{{hidden end}}

Below is example Python code to draw the sample:

<syntaxhighlight lang="python">
params = [a1, a2, ..., ak]
sample = [random.gammavariate(a, 1) for a in params]
sample = [v / sum(sample) for v in sample]
</syntaxhighlight>

This formulation is correct regardless of how the Gamma distributions are parameterized (shape/scale vs. shape/rate) because they are equivalent when scale and rate equal 1.0.