Dirichlet distribution

Properties

Moments

Let ${\displaystyle X=(X_1,\dots ,X_K)\sim \mathrm{Dir}(\mathit{\alpha })}$.

Let

$${\displaystyle {\alpha }_0=\sum _{i=1}^K{\alpha }_i.}$$

Then

$${\displaystyle \mathrm{E}[X_i]=\frac{{\alpha }_i}{{\alpha }_0},}$$ $${\displaystyle \mathrm{Var}[X_i]=\frac{{\alpha }_i({\alpha }_0-{\alpha }_i)}{{\alpha }_0^2({\alpha }_0+1)}.}$$

Furthermore, if ${\displaystyle i\ne j}$

$${\displaystyle \mathrm{Cov}[X_i,X_j]=\frac{-{\alpha }_i{\alpha }_j}{{\alpha }_0^2({\alpha }_0+1)}.}$$

The covariance matrix is singular.

More generally, moments of Dirichlet-distributed random variables can be expressed in the following way. For ${\displaystyle \mathit{t}=(t_1,\dots ,t_K)\in {\mathbb{R}}^K}$, denote by ${\displaystyle {\mathit{t}}^{\circ i}=(t_1^i,\dots ,t_K^i)}$ its i-th Hadamard power. Then,

${\displaystyle \mathrm{E}\left[(\mathit{t}\cdot \mathit{X}{)}^n\right]=\frac{n!\mathrm{\Gamma }({\alpha }_0)}{\mathrm{\Gamma }({\alpha }_0+n)}\sum \frac{{t_1}^{k_1}\cdots {t_K}^{k_K}}{k_1!\cdots k_K!}\prod _{i=1}^K\frac{\mathrm{\Gamma }({\alpha }_i+k_i)}{\mathrm{\Gamma }({\alpha }_i)}=\frac{n!\mathrm{\Gamma }({\alpha }_0)}{\mathrm{\Gamma }({\alpha }_0+n)}{Z}_n({\mathit{t}}^{\circ 1}\cdot \mathit{\alpha },\cdots ,{\mathit{t}}^{\circ n}\cdot \mathit{\alpha }),}$

where the sum is over non-negative integers ${\displaystyle k_1,\dots ,k_K}$ with ${\displaystyle n=k_1+\cdots +k_K}$, and ${\displaystyle Z_n}$ is the cycle index polynomial of the Symmetric group of degree n.

We have the special case ${\displaystyle \mathrm{E}\left[\mathit{t}\cdot \mathit{X}\right]=\frac{\mathit{t}\cdot \mathit{\alpha }}{{\alpha }_0}.}$

The multivariate analogue $\mathrm{E}\left[({\mathit{t}}_1\cdot \mathit{X}{)}^{n_1}\cdots ({\mathit{t}}_q\cdot \mathit{X}{)}^{n_q}\right]$ for vectors ${\displaystyle {\mathit{t}}_1,\dots ,{\mathit{t}}_q\in {\mathbb{R}}^K}$ can be expressed in terms of a color pattern of the exponents ${\displaystyle n_1,\dots ,n_q}$ in the sense of the Pólya enumeration theorem.

Particular cases include the simple computation

$${\displaystyle \mathrm{E}\left[\prod _{i=1}^K{X}_i^{{\beta }_i}\right]=\frac{B\left(\mathit{\alpha }+\mathit{\beta }\right)}{B\left(\mathit{\alpha }\right)}=\frac{\mathrm{\Gamma }\left(\sum _{i=1}^K{\alpha }_i\right)}{\mathrm{\Gamma }\left[\sum _{i=1}^K({\alpha }_i+{\beta }_i)\right]}\times \prod _{i=1}^K\frac{\mathrm{\Gamma }({\alpha }_i+{\beta }_i)}{\mathrm{\Gamma }({\alpha }_i)}.}$$

Mode

The mode of the distribution is the vector (x₁, ..., x_K) with

$${\displaystyle x_i=\frac{{\alpha }_i-1}{{\alpha }_0-K},{\alpha }_i>1.}$$

Marginal distributions

The marginal distributions are beta distributions:

$${\displaystyle X_i\sim \mathrm{Beta}({\alpha }_i,{\alpha }_0-{\alpha }_i).}$$

Also see § Related distributions below.

Conjugate to categorical or multinomial

The Dirichlet distribution is the conjugate prior distribution of the categorical distribution (a generic discrete probability distribution with a given number of possible outcomes) and multinomial distribution (the distribution over observed counts of each possible category in a set of categorically distributed observations). This means that if a data point has either a categorical or multinomial distribution, and the prior distribution of the distribution's parameter (the vector of probabilities that generates the data point) is distributed as a Dirichlet, then the posterior distribution of the parameter is also a Dirichlet. Intuitively, in such a case, starting from what we know about the parameter prior to observing the data point, we then can update our knowledge based on the data point and end up with a new distribution of the same form as the old one. This means that we can successively update our knowledge of a parameter by incorporating new observations one at a time, without running into mathematical difficulties.

Formally, this can be expressed as follows. Given a model

$${\displaystyle \begin{array}{rcccl}\mathit{\alpha }& =& \left({\alpha }_1,\dots ,{\alpha }_K\right)& =& \text{concentration hyperparameter}\\ \mathbf{p}\mid \mathit{\alpha }& =& \left(p_1,\dots ,p_K\right)& \sim & \mathrm{Dir}(K,\mathit{\alpha })\\ \mathbb{X}\mid \mathbf{p}& =& \left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_K\right)& \sim & \mathrm{Cat}(K,\mathbf{p})\end{array}}$$

then the following holds:

$${\displaystyle \begin{array}{rcccl}\mathbf{c}& =& \left(c_1,\dots ,c_K\right)& =& \text{number of occurrences of category }i\\ \mathbf{p}\mid \mathbb{X},\mathit{\alpha }& \sim & \mathrm{Dir}(K,\mathbf{c}+\mathit{\alpha })& =& \mathrm{Dir}\left(K,c_1+{\alpha }_1,\dots ,c_K+{\alpha }_K\right)\end{array}}$$

This relationship is used in Bayesian statistics to estimate the underlying parameter p of a categorical distribution given a collection of N samples. Intuitively, we can view the hyperprior vector α as pseudocounts, i.e. as representing the number of observations in each category that we have already seen. Then we simply add in the counts for all the new observations (the vector c) in order to derive the posterior distribution.

In Bayesian mixture models and other hierarchical Bayesian models with mixture components, Dirichlet distributions are commonly used as the prior distributions for the categorical variables appearing in the models. See the section on applications below for more information.

Relation to Dirichlet-multinomial distribution

In a model where a Dirichlet prior distribution is placed over a set of categorical-valued observations, the marginal joint distribution of the observations (i.e. the joint distribution of the observations, with the prior parameter marginalized out) is a Dirichlet-multinomial distribution. This distribution plays an important role in hierarchical Bayesian models, because when doing inference over such models using methods such as Gibbs sampling or variational Bayes, Dirichlet prior distributions are often marginalized out. See the article on this distribution for more details.

Entropy

If X is a ${\displaystyle \mathrm{Dir}(\mathit{\alpha })}$ random variable, the differential entropy of X (in nat units) is

$${\displaystyle h(\mathit{X})=\mathrm{E}[-\ln f(\mathit{X})]=\ln \mathrm{B}(\mathit{\alpha })+({\alpha }_0-K)\psi ({\alpha }_0)-\sum _{j=1}^K({\alpha }_j-1)\psi ({\alpha }_j)}$$

where ${\displaystyle \psi }$ is the digamma function.

The following formula for ${\displaystyle \mathrm{E}[\ln (X_i)]}$ can be used to derive the differential entropy above. Since the functions ${\displaystyle \ln (X_i)}$ are the sufficient statistics of the Dirichlet distribution, the exponential family differential identities can be used to get an analytic expression for the expectation of ${\displaystyle \ln (X_i)}$ (see equation (2.62) in ) and its associated covariance matrix:

$${\displaystyle \mathrm{E}[\ln (X_i)]=\psi ({\alpha }_i)-\psi ({\alpha }_0)}$$

and

$${\displaystyle \mathrm{Cov}[\ln (X_i),\ln (X_j)]={\psi }^{\prime }({\alpha }_i){\delta }_{ij}-{\psi }^{\prime }({\alpha }_0)}$$

where ${\displaystyle \psi }$ is the digamma function, ${\displaystyle {\psi }^{\prime }}$ is the trigamma function, and ${\displaystyle {\delta }_{ij}}$ is the Kronecker delta.

The spectrum of Rényi information for values other than ${\displaystyle \lambda =1}$ is given by

$${\displaystyle F_R(\lambda )=(1-\lambda {)}^{-1}\left(-\lambda \log \mathrm{B}(\mathit{\alpha })+\sum _{i=1}^K\log \mathrm{\Gamma }(\lambda ({\alpha }_i-1)+1)-\log \mathrm{\Gamma }(\lambda ({\alpha }_0-K)+K)\right)}$$

and the information entropy is the limit as ${\displaystyle \lambda }$ goes to 1.

Another related interesting measure is the entropy of a discrete categorical (one-of-K binary) vector Z with probability-mass distribution X, i.e., ${\displaystyle P(Z_i=1,Z_{j\ne i}=0|\mathit{X})=X_i}$. The conditional information entropy of Z, given X is

$${\displaystyle S(\mathit{X})=H(\mathit{Z}|\mathit{X})={\mathrm{E}}_{\mathit{Z}}[-\log P(\mathit{Z}|\mathit{X})]=\sum _{i=1}^K-X_i\log X_i}$$

This function of X is a scalar random variable. If X has a symmetric Dirichlet distribution with all ${\displaystyle {\alpha }_i=\alpha }$, the expected value of the entropy (in nat units) is

$${\displaystyle \mathrm{E}[S(\mathit{X})]=\sum _{i=1}^K\mathrm{E}[-X_i\ln X_i]=\psi (K\alpha +1)-\psi (\alpha +1)}$$

Kullback–Leibler divergence

The Kullback–Leibler (KL) divergence between two Dirichlet distributions, ${\displaystyle \text{Dir}(\mathit{\alpha })}$ and ${\displaystyle \text{Dir}(\mathit{\beta })}$, over the same simplex is:

${\displaystyle \begin{array}{cc}{D}_{\mathrm{K}\mathrm{L}}(\mathrm{D}\mathrm{i}\mathrm{r}(\mathit{\alpha })\Vert \mathrm{D}\mathrm{i}\mathrm{r}(\mathit{\beta }))& =\log \frac{\mathrm{\Gamma }\left(\sum _{i=1}^K{\alpha }_i\right)}{\mathrm{\Gamma }\left(\sum _{i=1}^K{\beta }_i\right)}+\sum _{i=1}^K\left[\log \frac{\mathrm{\Gamma }({\beta }_i)}{\mathrm{\Gamma }({\alpha }_i)}+({\alpha }_i-{\beta }_i)\left(\psi ({\alpha }_i)-\psi \left(\sum _{j=1}^K{\alpha }_j\right)\right)\right]\end{array}}$

Aggregation

If

$${\displaystyle X=(X_1,\dots ,X_K)\sim \mathrm{Dir}({\alpha }_1,\dots ,{\alpha }_K)}$$

then, if the random variables with subscripts i and j are dropped from the vector and replaced by their sum,

$${\displaystyle X^{\prime }=(X_1,\dots ,X_i+X_j,\dots ,X_K)\sim \mathrm{Dir}({\alpha }_1,\dots ,{\alpha }_i+{\alpha }_j,\dots ,{\alpha }_K).}$$

This aggregation property may be used to derive the marginal distribution of ${\displaystyle X_i}$ mentioned above.

Neutrality

If ${\displaystyle X=(X_1,\dots ,X_K)\sim \mathrm{Dir}(\mathit{\alpha })}$, then the vector X is said to be neutral in the sense that X_K is independent of ${\displaystyle X^{(-K)}}$ where

$${\displaystyle X^{(-K)}=\left(\frac{X_1}{1-X_K},\frac{X_2}{1-X_K},\dots ,\frac{X_{K-1}}{1-X_K}\right),}$$

and similarly for removing any of ${\displaystyle X_2,\dots ,X_{K-1}}$. Observe that any permutation of X is also neutral (a property not possessed by samples drawn from a generalized Dirichlet distribution).

Combining this with the property of aggregation it follows that X_j + ... + X_K is independent of ${\displaystyle \left(\frac{X_1}{X_1+\cdots +X_{j-1}},\frac{X_2}{X_1+\cdots +X_{j-1}},\dots ,\frac{X_{j-1}}{X_1+\cdots +X_{j-1}}\right)}$. In fact it is true, further, for the Dirichlet distribution, that for ${\displaystyle 3\le j\le K-1}$, the pair ${\displaystyle \left(X_1+\cdots +X_{j-1},X_j+\cdots +X_K\right)}$, and the two vectors ${\displaystyle \left(\frac{X_1}{X_1+\cdots +X_{j-1}},\frac{X_2}{X_1+\cdots +X_{j-1}},\dots ,\frac{X_{j-1}}{X_1+\cdots +X_{j-1}}\right)}$ and ${\displaystyle \left(\frac{X_j}{X_j+\cdots +X_K},\frac{X_{j+1}}{X_j+\cdots +X_K},\dots ,\frac{X_K}{X_j+\cdots +X_K}\right)}$, viewed as triple of normalised random vectors, are mutually independent. The analogous result is true for partition of the indices {1, 2, ..., K} into any other pair of non-singleton subsets.

Characteristic function

The characteristic function of the Dirichlet distribution is a confluent form of the Lauricella hypergeometric series. It is given by Phillips as

$${\displaystyle CF\left(s_1,\dots ,s_{K-1}\right)=\mathrm{E}\left(e^{i\left(s_1{X}_1+\cdots +s_{K-1}{X}_{K-1}\right)}\right)={\mathrm{\Psi }}^{\left[K-1\right]}({\alpha }_1,\dots ,{\alpha }_{K-1};{\alpha }_0;i{s}_1,\dots ,i{s}_{K-1})}$$

where

$${\displaystyle {\mathrm{\Psi }}^{[m]}(a_1,\dots ,a_m;c;z_1,\dots z_m)=\sum \frac{(a_1{)}_{k_1}\cdots (a_m{)}_{k_m}{z}_1^{k_1}\cdots z_m^{k_m}}{(c{)}_k{k}_1!\cdots k_m!}.}$$

The sum is over non-negative integers ${\displaystyle k_1,\dots ,k_m}$ and ${\displaystyle k=k_1+\cdots +k_m}$. Phillips goes on to state that this form is "inconvenient for numerical calculation" and gives an alternative in terms of a complex path integral:

$${\displaystyle {\mathrm{\Psi }}^{[m]}=\frac{\mathrm{\Gamma }(c)}{2\pi i}{\int }_L{e}^t{t}^{a_1+\cdots +a_m-c}\prod _{j=1}^m(t-z_j{)}^{-a_j}dt}$$

where L denotes any path in the complex plane originating at ${\displaystyle -\mathrm{\infty }}$, encircling in the positive direction all the singularities of the integrand and returning to ${\displaystyle -\mathrm{\infty }}$.

Inequality

Probability density function ${\displaystyle f\left(x_1,\dots ,x_{K-1};{\alpha }_1,\dots ,{\alpha }_K\right)}$ plays a key role in a multifunctional inequality which implies various bounds for the Dirichlet distribution.

Another inequality relates the moment-generating function of the Dirichlet distribution to the convex conjugate of the scaled reversed Kullback-Leibler divergence:

$${\displaystyle \log \mathrm{E}\left(\exp \sum _{i=1}^K{s}_i{X}_i\right)\le \underset{p}{sup}\sum _{i=1}^K\left(p_i{s}_i-{\alpha }_i\log \left(\frac{{\alpha }_i}{{\alpha }_0{p}_i}\right)\right),}$$ where the supremum is taken over p spanning the (K − 1)-simplex.

Related distributions

When ${\displaystyle \mathit{X}=(X_1,\dots ,X_K)\sim \mathrm{Dir}\left({\alpha }_1,\dots ,{\alpha }_K\right)}$, the marginal distribution of each component ${\displaystyle X_i\sim \mathrm{Beta}({\alpha }_i,{\alpha }_0-{\alpha }_i)}$, a Beta distribution. In particular, if K = 2 then ${\displaystyle X_1\sim \mathrm{Beta}({\alpha }_1,{\alpha }_2)}$ is equivalent to ${\displaystyle \mathit{X}=(X_1,1-X_1)\sim \mathrm{Dir}\left({\alpha }_1,{\alpha }_2\right)}$.

For K independently distributed Gamma distributions:

$${\displaystyle Y_1\sim \mathrm{Gamma}({\alpha }_1,1),\dots ,Y_K\sim \mathrm{Gamma}({\alpha }_K,1)}$$

we have writing ${\displaystyle {\alpha }_0:=\sum _{i=1}^K{\alpha }_i}$:

$${\displaystyle V=\sum _{i=1}^K{Y}_i\sim \mathrm{Gamma}\left({\alpha }_0,1\right),}$$ $${\displaystyle X=(X_1,\dots ,X_K)=\left(\frac{Y_1}{V},\dots ,\frac{Y_K}{V}\right)\sim \mathrm{Dir}\left({\alpha }_1,\dots ,{\alpha }_K\right).}$$

Although the X_is are not independent from one another, they can be seen to be generated from a set of K independent gamma random variables. Unfortunately, since the sum V is lost in forming X (in fact it can be shown that V is stochastically independent of X), it is not possible to recover the original gamma random variables from these values alone. Nevertheless, because independent random variables are simpler to work with, this reparametrization can still be useful for proofs about properties of the Dirichlet distribution.

Conjugate prior of the Dirichlet distribution

Because the Dirichlet distribution is an exponential family distribution it has a conjugate prior. The conjugate prior is of the form:

$${\displaystyle \mathrm{CD}(\mathit{\alpha }\mid \mathit{v},\eta )\propto {\left(\frac{1}{\mathrm{B}(\mathit{\alpha })}\right)}^{\eta }\exp \left(-\sum _k{v}_k{\alpha }_k\right).}$$

Here ${\displaystyle \mathit{v}}$ is a K-dimensional real vector and ${\displaystyle \eta }$ is a scalar parameter. The domain of ${\displaystyle (\mathit{v},\eta )}$ is restricted to the set of parameters for which the above unnormalized density function can be normalized. The (necessary and sufficient) condition is:

$${\displaystyle \mathrm{\forall }k{v}_k>0\text{ and }\eta >-1\text{ and }(\eta \le 0\text{ or }\sum _k\exp -\frac{v_k}{\eta }<1)}$$

The conjugation property can be expressed as

if [prior: ${\displaystyle \mathit{\alpha }\sim \mathrm{CD}(\cdot \mid \mathit{v},\eta )}$] and [observation: ${\displaystyle \mathit{x}\mid \mathit{\alpha }\sim \mathrm{Dirichlet}(\cdot \mid \mathit{\alpha })}$] then [posterior: ${\displaystyle \mathit{\alpha }\mid \mathit{x}\sim \mathrm{CD}(\cdot \mid \mathit{v}-\log \mathit{x},\eta +1)}$].

In the published literature there is no practical algorithm to efficiently generate samples from ${\displaystyle \mathrm{CD}(\mathit{\alpha }\mid \mathit{v},\eta )}$.

Generalization by scaling and translation of log-probabilities

As noted above, Dirichlet variates can be generated by normalizing independent gamma variates. If instead one normalizes generalized gamma variates, one obtains variates from the simplicial generalized beta distribution (SGB). On the other hand, SGB variates can also be obtained by applying the softmax function to scaled and translated logarithms of Dirichlet variates. Specifically, let ${\displaystyle \mathbf{x}=(x_1,\dots ,x_K)\sim \mathrm{Dir}(\mathit{\alpha })}$ and let ${\displaystyle \mathbf{y}=(y_1,\dots ,y_K)}$, where applying the logarithm elementwise: $${\displaystyle \mathbf{y}=\mathrm{softmax}(a^{-1}\log \mathbf{x}+\log \mathbf{b})⟺\mathbf{x}=\mathrm{softmax}(a\log \mathbf{y}-a\log \mathbf{b})}$$ or $${\displaystyle y_k=\frac{b_k{x}_k^{1/a}}{\sum _{i=1}^K{b}_i{x}_i^{1/a}}⟺x_k=\frac{(y_k/b_k{)}^a}{\sum _{i=1}^K(y_i/b_i{)}^a}}$$ where ${\displaystyle a>0}$ and ${\displaystyle \mathbf{b}=(b_1,\dots ,b_K)}$, with all ${\displaystyle b_k>0}$, then ${\displaystyle \mathbf{y}\sim \mathrm{SGB}(a,\mathbf{b},\mathit{\alpha })}$. The SGB density function can be derived by noting that the transformation ${\displaystyle \mathbf{x}\mapsto \mathbf{y}}$, which is a bijection from the simplex to itself, induces a differential volume change factor of: $${\displaystyle R(\mathbf{y},a,\mathbf{b})=a^{1-K}\prod _{k=1}^K\frac{y_k}{x_k}}$$ where it is understood that ${\displaystyle \mathbf{x}}$ is recovered as a function of ${\displaystyle \mathbf{y}}$, as shown above. This facilitates writing the SGB density in terms of the Dirichlet density, as: $${\displaystyle f_{\text{SGB}}(\mathbf{y}\mid a,\mathbf{b},\mathit{\alpha })=\frac{f_{\text{Dir}}(\mathbf{x}\mid \mathit{\alpha })}{R(\mathbf{y},a,\mathbf{b})}}$$ This generalization of the Dirichlet density, via a change of variables, is closely related to a normalizing flow, while the differential volume change is not given by the Jacobian determinant of ${\displaystyle \mathbf{x}\mapsto \mathbf{y}:{\mathbb{R}}^K\to {\mathbb{R}}^K}$ which is zero, but by the Jacobian determinant of ${\displaystyle (x_1,\dots ,x_{K-1})\mapsto \mathbf{(}{y}_1,\dots ,y_{K-1})}$, as explained in more detail at Normalizing flow § Simplex flow.

For further insight into the interaction between the Dirichlet shape parameters ${\displaystyle \mathit{\alpha }}$, and the transformation parameters ${\displaystyle a,\mathbf{b}}$, it may be helpful to consider the logarithmic marginals, ${\displaystyle \log \frac{x_k}{1-x_k}}$, which follow the logistic-beta distribution, ${\displaystyle B_{\sigma }({\alpha }_k,\sum _{i\ne k}{\alpha }_i)}$. See in particular the sections on tail behaviour and generalization with location and scale parameters.

Application

When ${\displaystyle b_1=b_2=\cdots =b_K}$, then the transformation simplifies to ${\displaystyle \mathbf{x}\mapsto \mathrm{softmax}(a^{-1}\log \mathbf{x})}$, which is known as temperature scaling in machine learning, where it is used as a calibration transform for multiclass probabilistic classifiers. Traditionally the temperature parameter (${\displaystyle a}$ here) is learnt discriminatively by minimizing multiclass cross-entropy over a supervised calibration data set with known class labels. But the above PDF transformation mechanism can be used to facilitate also the design of generatively trained calibration models with a temperature scaling component.

Occurrence and applications

Bayesian models

Dirichlet distributions are most commonly used as the prior distribution of categorical variables or multinomial variables in Bayesian mixture models and other hierarchical Bayesian models. (In many fields, such as in natural language processing, categorical variables are often imprecisely called "multinomial variables". Such a usage is unlikely to cause confusion, just as when Bernoulli distributions and binomial distributions are commonly conflated.)

Inference over hierarchical Bayesian models is often done using Gibbs sampling, and in such a case, instances of the Dirichlet distribution are typically marginalized out of the model by integrating out the Dirichlet random variable. This causes the various categorical variables drawn from the same Dirichlet random variable to become correlated, and the joint distribution over them assumes a Dirichlet-multinomial distribution, conditioned on the hyperparameters of the Dirichlet distribution (the concentration parameters). One of the reasons for doing this is that Gibbs sampling of the Dirichlet-multinomial distribution is extremely easy; see that article for more information.

Intuitive interpretations of the parameters

The concentration parameter

Dirichlet distributions are very often used as prior distributions in Bayesian inference. The simplest and perhaps most common type of Dirichlet prior is the symmetric Dirichlet distribution, where all parameters are equal. This corresponds to the case where you have no prior information to favor one component over any other. As described above, the single value α to which all parameters are set is called the concentration parameter. If the sample space of the Dirichlet distribution is interpreted as a discrete probability distribution, then intuitively the concentration parameter can be thought of as determining how "concentrated" the probability mass of the Dirichlet distribution to its center, leading to samples with mass dispersed almost equally among all components, i.e., with a value much less than 1, the mass will be highly concentrated in a few components, and all the rest will have almost no mass, and with a value much greater than 1, the mass will be dispersed almost equally among all the components. See the article on the concentration parameter for further discussion.

String cutting

One example use of the Dirichlet distribution is if one wanted to cut strings (each of initial length 1.0) into K pieces with different lengths, where each piece had a designated average length, but allowing some variation in the relative sizes of the pieces. Recall that ${\displaystyle {\alpha }_0=\sum _{i=1}^K{\alpha }_i.}$ The ${\displaystyle {\alpha }_i/{\alpha }_0}$ values specify the mean lengths of the cut pieces of string resulting from the distribution. The variance around this mean varies inversely with ${\displaystyle {\alpha }_0}$.

Pólya's urn

Consider an urn containing balls of K different colors. Initially, the urn contains α₁ balls of color 1, α₂ balls of color 2, and so on. Now perform N draws from the urn, where after each draw, the ball is placed back into the urn with an additional ball of the same color. In the limit as N approaches infinity, the proportions of different colored balls in the urn will be distributed as Dir(α₁, ..., α_K).

For a formal proof, note that the proportions of the different colored balls form a bounded [0,1]^K-valued martingale, hence by the martingale convergence theorem, these proportions converge almost surely and in mean to a limiting random vector. To see that this limiting vector has the above Dirichlet distribution, check that all mixed moments agree.

Each draw from the urn modifies the probability of drawing a ball of any one color from the urn in the future. This modification diminishes with the number of draws, since the relative effect of adding a new ball to the urn diminishes as the urn accumulates increasing numbers of balls.

Random variate generation

From gamma distribution

With a source of Gamma-distributed random variates, one can easily sample a random vector ${\displaystyle x=(x_1,\dots ,x_K)}$ from the K-dimensional Dirichlet distribution with parameters ${\displaystyle ({\alpha }_1,\dots ,{\alpha }_K)}$ . First, draw K independent random samples ${\displaystyle y_1,\dots ,y_K}$ from Gamma distributions each with density

$${\displaystyle \mathrm{Gamma}({\alpha }_i,1)=\frac{y_i^{{\alpha }_i-1}{e}^{-y_i}}{\mathrm{\Gamma }({\alpha }_i)},}$$

and then set

$${\displaystyle x_i=\frac{y_i}{\sum _{j=1}^K{y}_j}.}$$

Below is example Python code to draw the sample:

<span id="mwAo4"></span><span class="n" id="mwAo8">params</span> <span class="o" id="mwApA">=</span> <span class="p" id="mwApE">[</span><span class="n" id="mwApI">a1</span><span class="p" id="mwApM">,</span> <span class="n" id="mwApQ">a2</span><span class="p" id="mwApU">,</span> <span class="o" id="mwApY">...</span><span class="p" id="mwApc">,</span> <span class="n" id="mwApg">ak</span><span class="p" id="mwApk">]</span> <span class="n" id="mwApo">sample</span> <span class="o" id="mwAps">=</span> <span class="p" id="mwApw">[</span><span class="n" id="mwAp0">random</span><span class="o" id="mwAp4">.</span><span class="n" id="mwAp8">gammavariate</span><span class="p" id="mwAqA">(</span><span class="n" id="mwAqE">a</span><span class="p" id="mwAqI">,</span> <span class="mi" id="mwAqM">1</span><span class="p" id="mwAqQ">)</span> <span class="k" id="mwAqU">for</span> <span class="n" id="mwAqY">a</span> <span class="ow" id="mwAqc">in</span> <span class="n" id="mwAqg">params</span><span class="p" id="mwAqk">]</span> <span class="n" id="mwAqo">sample</span> <span class="o" id="mwAqs">=</span> <span class="p" id="mwAqw">[</span><span class="n" id="mwAq0">v</span> <span class="o" id="mwAq4">/</span> <span class="nb" id="mwAq8">sum</span><span class="p" id="mwArA">(</span><span class="n" id="mwArE">sample</span><span class="p" id="mwArI">)</span> <span class="k" id="mwArM">for</span> <span class="n" id="mwArQ">v</span> <span class="ow" id="mwArU">in</span> <span class="n" id="mwArY">sample</span><span class="p" id="mwArc">]</span>

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.

From marginal beta distributions

A less efficient algorithm relies on the univariate marginal and conditional distributions being beta and proceeds as follows. Simulate ${\displaystyle x_1}$ from

$${\displaystyle \text{Beta}\left({\alpha }_1,\sum _{i=2}^K{\alpha }_i\right)}$$

Then simulate ${\displaystyle x_2,\dots ,x_{K-1}}$ in order, as follows. For ${\displaystyle j=2,\dots ,K-1}$, simulate ${\displaystyle {\varphi }_j}$ from

$${\displaystyle \text{Beta}\left({\alpha }_j,\sum _{i=j+1}^K{\alpha }_i\right),}$$

and let

$${\displaystyle x_j=\left(1-\sum _{i=1}^{j-1}{x}_i\right){\varphi }_j.}$$

Finally, set

$${\displaystyle x_K=1-\sum _{i=1}^{K-1}{x}_i.}$$

This iterative procedure corresponds closely to the "string cutting" intuition described above.

Below is example Python code to draw the sample:

<span id="mwAtI"></span><span class="n" id="mwAtM">params</span> <span class="o" id="mwAtQ">=</span> <span class="p" id="mwAtU">[</span><span class="n" id="mwAtY">a1</span><span class="p" id="mwAtc">,</span> <span class="n" id="mwAtg">a2</span><span class="p" id="mwAtk">,</span> <span class="o" id="mwAto">...</span><span class="p" id="mwAts">,</span> <span class="n" id="mwAtw">ak</span><span class="p" id="mwAt0">]</span> <span class="n" id="mwAt4">xs</span> <span class="o" id="mwAt8">=</span> <span class="p" id="mwAuA">[</span><span class="n" id="mwAuE">random</span><span class="o" id="mwAuI">.</span><span class="n" id="mwAuM">betavariate</span><span class="p" id="mwAuQ">(</span><span class="n" id="mwAuU">params</span><span class="p" id="mwAuY">[</span><span class="mi" id="mwAuc">0</span><span class="p" id="mwAug">],</span> <span class="nb" id="mwAuk">sum</span><span class="p" id="mwAuo">(</span><span class="n" id="mwAus">params</span><span class="p" id="mwAuw">[</span><span class="mi" id="mwAu0">1</span><span class="p" id="mwAu4">:]))]</span> <span class="k" id="mwAu8">for</span> <span class="n" id="mwAvA">j</span> <span class="ow" id="mwAvE">in</span> <span class="nb" id="mwAvI">range</span><span class="p" id="mwAvM">(</span><span class="mi" id="mwAvQ">1</span><span class="p" id="mwAvU">,</span> <span class="nb" id="mwAvY">len</span><span class="p" id="mwAvc">(</span><span class="n" id="mwAvg">params</span><span class="p" id="mwAvk">)</span> <span class="o" id="mwAvo">-</span> <span class="mi" id="mwAvs">1</span><span class="p" id="mwAvw">):</span> <span class="n" id="mwAv0">phi</span> <span class="o" id="mwAv4">=</span> <span class="n" id="mwAv8">random</span><span class="o" id="mwAwA">.</span><span class="n" id="mwAwE">betavariate</span><span class="p" id="mwAwI">(</span><span class="n" id="mwAwM">params</span><span class="p" id="mwAwQ">[</span><span class="n" id="mwAwU">j</span><span class="p" id="mwAwY">],</span> <span class="nb" id="mwAwc">sum</span><span class="p" id="mwAwg">(</span><span class="n" id="mwAwk">params</span><span class="p" id="mwAwo">[</span><span class="n" id="mwAws">j</span> <span class="o" id="mwAww">+</span> <span class="mi" id="mwAw0">1</span> <span class="p" id="mwAw4">:]))</span> <span class="n" id="mwAw8">xs</span><span class="o" id="mwAxA">.</span><span class="n" id="mwAxE">append</span><span class="p" id="mwAxI">((</span><span class="mi" id="mwAxM">1</span> <span class="o" id="mwAxQ">-</span> <span class="nb" id="mwAxU">sum</span><span class="p" id="mwAxY">(</span><span class="n" id="mwAxc">xs</span><span class="p" id="mwAxg">))</span> <span class="o" id="mwAxk">*</span> <span class="n" id="mwAxo">phi</span><span class="p" id="mwAxs">)</span> <span class="n" id="mwAxw">xs</span><span class="o" id="mwAx0">.</span><span class="n" id="mwAx4">append</span><span class="p" id="mwAx8">(</span><span class="mi" id="mwAyA">1</span> <span class="o" id="mwAyE">-</span> <span class="nb" id="mwAyI">sum</span><span class="p" id="mwAyM">(</span><span class="n" id="mwAyQ">xs</span><span class="p" id="mwAyU">))</span>

When each alpha is 1

When α₁ = ... = α_K = 1, a sample from the distribution can be found by randomly drawing a set of K − 1 values independently and uniformly from the interval [0, 1], adding the values 0 and 1 to the set to make it have K + 1 values, sorting the set, and computing the difference between each pair of order-adjacent values, to give x₁, ..., x_K.

When each alpha is 1/2 and relationship to the hypersphere

When α₁ = ... = α_K = 1/2, a sample from the distribution can be found by randomly drawing K values independently from the standard normal distribution, squaring these values, and normalizing them by dividing by their sum, to give x₁, ..., x_K.

A point (x₁, ..., x_K) can be drawn uniformly at random from the (K−1)-dimensional unit hypersphere (which is the surface of a K-dimensional hyperball) via a similar procedure. Randomly draw K values independently from the standard normal distribution and normalize these coordinate values by dividing each by the constant that is the square root of the sum of their squares.