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<d 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<d\text{ and }i \neq k\text{ , }\frac{\partial z_k}{\partial x_i} &= 0 \\ \text{for }k=d\text{, }\frac{\partial z_d}{\partial x_i} &= \frac{\partial v(1-\sum_{k=1}^{d-1}x_k)}{\partial x_i} \\ &=-v \\ \frac{\partial z_d}{\partial v} &= \frac{\partial v(1-\sum_{k=1}^{d-1}x_k)}{\partial v} \\ &=1-\sum_{k=1}^{d-1}x_k \end{aligned}

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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s