# Sampling from Dirichlet Distribution using Gamma distributed samples

There is an algorithm to generate Dirichlet samples using a sampler for Gamma distribution for any $\alpha > 0$ and $\beta > 0$. We will generate Gamma distributed variables $z_k \sim \text{gamma}(\alpha_k,1)$, for $k \in {1,\cdots,d}$, and do the following variable transformation to get Dirichlet samples $x_k = \frac{z_k}{\sum_k z_k}$. First we should demonstrate that this transformation results in Dirichlet distributed samples.

Consider the following tranformation $(z_1,\cdots,z_d) \leftarrow (x_1,\cdots,x_d,v)$, where $x_k = \frac{z_k}{\sum_k z_k}$ and $v = {\sum_k z_k}$. We can rewrite this transformation as $(x_1,\cdots,x_d,v)=h(z_1,\cdots,z_d)$, where $x_k = \frac{z_k}{v}$ and $v = {\sum_k z_k}$. Also we can imediatly calculate the inverse transformation $(z_1,\cdots,z_d)=h^{-1}(x_1,\cdots,x_d,v)$, with $z_k=v x_k$. From the transformation definition we know that ${\sum_{k=1}^d x_k=1}$, implying that $x_d = 1-\sum_{k=1}^{d-1} x_k$ and $z_d=v(1-\sum_{k=1}^{d-1}x_k)$.

Now we can apply the formula for transformation of variables to $z_{1:d}$,$x_{1:(d-1)}$ and $v$ (using the notation that $f_X$ is the pdf of the random variables $(x_{1:(d-1)},v)$ and  $f_Z$ is the pdf of the random variables $z_{1:d}$):

\begin{aligned} f_{Z}(z_{1:d})&=\prod_{k=1}^d \frac{z_k^{\alpha_k-1} e^{-z_k}}{\Gamma(\alpha_k)} \\ f_{X}(x_{1:(d-1)},v)&=f_{Z}(z_{1:d})J(z_{1:d}) \\ &=f_{Z}(h^{-1}(x_{1:(d-1)},v))\left|\frac{\partial z_{1:d}}{\partial (x_{1:(d-1)},v)}\right| \end{aligned}

So we need to compute the determinant of the Jacobian  $J(z_{1:d})$ of the transformation.

\begin{aligned} J(z_{1:d})&=\det \left( \begin{bmatrix} \frac{\partial z_1}{\partial x_1} & \cdots & \frac{\partial z_1}{\partial x_{d-1}} & \frac{\partial z_1}{\partial v} \\ \frac{\partial z_2}{\partial x_1} & \cdots & \frac{\partial z_2}{\partial x_{d-1}} & \frac{\partial z_2}{\partial v} \\ \vdots & \ddots & \cdots & \vdots \\ \frac{\partial z_d}{\partial x_1} & \cdots & \frac{\partial z_d}{\partial x_{d-1}} & \frac{\partial z_d}{\partial v} \end{bmatrix} \right) \end{aligned}

Computing the individual terms, for  $k and  $i=k$,

\begin{aligned} \frac{\partial z_k}{\partial x_i} &= \frac{\partial vx_k}{\partial x_k}\\ &=v \\ \frac{\partial z_k}{\partial v} &= \frac{\partial vx_k}{\partial v} \\ &=x_k \\ \text{for }k

In matricial forms we have:

\begin{aligned} J(z_{1:d})&=\det \left( \begin{bmatrix} v & 0 & 0 & \cdots & x_1 \\ 0 & v & 0 & \cdots & x_2 \\ 0 & 0 & v & \cdots & x_2 \\ \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\ddots}{v} & \underset{\vdots}{x_k} \\ - v & -v & -v & \cdots & 1-\sum_{k=1}^{d-1}x_k \end{bmatrix} \right) \\ \text{summing first line }&\text{with the last line}\\ &=\det \left( \begin{bmatrix} v & 0 & 0 & \cdots & x_1 \\ 0 & v & 0 & \cdots & x_2 \\ 0 & 0 & v & \cdots & x_2 \\ \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\ddots}{v} & \underset{\vdots}{x_k} \\ 0 & -v & -v & \cdots & 1+x_1-\sum_{k=1}^{d-1}x_k \end{bmatrix} \right) \\ \text{summing all the other lines }&\text{with the last line}\\ &=\det \left( \begin{bmatrix} v & 0 & 0 & \cdots & x_1 \\ 0 & v & 0 & \cdots & x_2 \\ 0 & 0 & v & \cdots & x_2 \\ \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\vdots}{0} & \underset{\ddots}{v} & \underset{\vdots}{x_k} \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix} \right) \\ &=v^{d-1} \\ \end{aligned}
Now with this result from the determinant of the Jacobian, we will apply it in the transformation of random variables formula for the pdf and marginilize  $v$ to get the pdf of  $x_{1:(d-1)}$.

\begin{aligned} \Rightarrow f_{X}(x_{1:(d-1)},v)&=f_{Z}(z_{1:d}) v^{d-1} \\ &= v^{d-1} \prod_{k=1}^d \frac{{(vx_k)}^{\alpha_k-1} e^{-vx_k}}{\Gamma(\alpha_k)} \\ \Rightarrow f_{X}(x_{1:(d-1)})&= \int_v f_{X}(x_{1:(d-1)},v) dv \\ &=\prod_{k=1}^d \frac{x_k^{\alpha_k-1}}{\Gamma(\alpha_k) } \int_v v^{d-1} \prod_{k=1}^d v^{\alpha_k-1} e^{-vx_k} dv \\ &=\prod_{k=1}^d \frac{x_k^{\alpha_k-1}}{\Gamma(\alpha_k) } \int_v v^{d-1-d+\sum_k \alpha_k} e^{-v \sum_k x_k} dv \\ &=\prod_{k=1}^d \frac{x_k^{\alpha_k-1}}{\Gamma(\alpha_k) } \underbrace{\int_v v^{-1+\sum_k \alpha_k} e^{-v} dv}_{\Gamma(\sum_k \alpha_k)} \\ \Rightarrow f_{X}(x_{1:(d-1)}) &=\Gamma(\sum_k \alpha_k) \prod_{k=1}^d \frac{x_k^{\alpha_k-1}}{\Gamma(\alpha_k) } \square \\ \end{aligned}
This finish our demonstration that the given transformation of Gamma distributed random variables results in Dirichlet distributed random variables. Notice that $x_d$ is totally determined by  $x_{1:(d-1)}$ ($\sum_k x_k = 1$), so we can actually express the right side of the equation as a formula of only $x_{1:(d-1)}$.

If drawing $d$ samples from a Gamma distribution can be performed with complexity $O(g(d))$, than $d$ samples from a Dirichlet distribution can be computed with the same $O(g(d))$ complexity.