View source for Dirichlet-multinomial distribution

{{Short description|Distributions in probability theory}}
{{Probability distribution|
  pdf_image  =|
  cdf_image  =|
  name       =Dirichlet-Multinomial|
  type       =mass|

  parameters =<math>n \in \{0, 1, 2, \ldots\}</math> number of trials<br /><math>\alpha_1, \ldots, \alpha_{K}>0, \alpha_0 = \sum \alpha_k </math>  |
  notation = <math>\mathrm{DirMult}(n,\boldsymbol{\alpha})</math>|
  support    =<math>x_i \in \{0,\dots,n\}</math><br><math>\Sigma x_i = n\!, i=1...K</math>|
  pdf        =<math>\frac{\Gamma \left( \alpha_0 \right) \Gamma \left(n+1\right)}
{\Gamma\left(n+\alpha_0 \right)}\prod_{k=1}^K\frac{\Gamma(x_{k}+\alpha_{k})}{\Gamma(\alpha_{k}) \Gamma \left(x_k+1\right)} </math> <ref name=Gluesenkamp2018 />|
  cdf        =|
  mean       =<math>\operatorname{E}(X_i) = n\frac{\alpha_i}{\alpha_0}</math>|
  median     =|
  mode       =|
  variance   =<math>\operatorname{Var}(X_i) = n\frac{\alpha_i}{ \alpha_0}\left(1-\frac{\alpha_i}{ \alpha_0}\right)\left(\frac{n+\alpha_0}{1+\alpha_0}\right)</math><br><math>\textstyle {\mathrm{Cov}}(X_i,X_k) = - n \frac{\alpha_i}{\alpha_0} \frac{\alpha_k}{\alpha_0} \left(\frac{n+\alpha_0}{1+\alpha_0}\right)~~(i\neq k)</math>||
||
  skewness   =|
  kurtosis   =|
  entropy    =|
  mgf        =<math>\operatorname{E}(\prod\limits_{k=1}^{K} {e}^{t_k\cdot x_k}) = \frac{\Gamma(\alpha_0) \Gamma(n+1)}{\Gamma(n+\alpha_0)} \cdot D_n(\boldsymbol{\alpha},(e^{t_1}, ..., e^{t_K}))</math><br>with<br><math> D_{n}= \frac{1}{n}\sum\limits_{u=1}^{n} \left[\left(\sum\limits_{k=1}^{K} \alpha_{k} \cdot {e}^{t_k\cdot u} \right) D_{n-u}\right]  , D_0=1 </math> <ref name=Gluesenkamp2018>{{cite journal |last1=Glüsenkamp |first=T.|year=2018|title=Probabilistic treatment of the uncertainty from the finite size of weighted Monte Carlo data |journal=EPJ Plus |volume=133 |issue=6|pages=218|doi=10.1140/epjp/i2018-12042-x |arxiv=1712.01293|bibcode=2018EPJP..133..218G|s2cid=125665629}}</ref>|
  char       =<math>\operatorname{E}(\prod\limits_{k=1}^{K} {e}^{i t_k\cdot x_k}) = \frac{\Gamma(\alpha_0) \Gamma(n+1)}{\Gamma(n+\alpha_0)} \cdot D_n(\boldsymbol{\alpha},(e^{it_1}, ..., e^{it_K}))</math> 
<br>with<br>
<math> D_{n}= \frac{1}{n}\sum\limits_{u=1}^{n} \left[\left(\sum\limits_{k=1}^{K} \alpha_{k} \cdot {e}^{it_k\cdot u} \right) D_{n-u}\right]  , D_0=1 </math><ref name=Gluesenkamp2018 />|
  pgf =<math>\operatorname{E}(\prod\limits_{k=1}^{K} {z_k}^{x_k}) = \frac{\Gamma(\alpha_0) \Gamma(n+1)}{\Gamma(n+\alpha_0)} \cdot D_n(\boldsymbol{\alpha},\mathbf{z})</math> 
<br>with<br>
<math> D_{n}= \frac{1}{n}\sum\limits_{u=1}^{n} \left[\left(\sum\limits_{k=1}^{K} \alpha_{k} \cdot {z_{k}}^u \right) D_{n-u}\right] , D_0=1  </math> <ref name=Gluesenkamp2018 />|
  conjugate  =|
}}
In [[probability theory]] and [[statistics]], the '''Dirichlet-multinomial distribution''' is a family of discrete multivariate [[probability distribution]]s on a finite support of non-negative integers. It is also called the '''Dirichlet compound multinomial distribution''' ('''DCM''') or '''multivariate Pólya distribution''' (after [[George Pólya]]). It is a [[compound probability distribution]], where a probability vector '''p''' is drawn from a [[Dirichlet distribution]] with parameter vector <math>\boldsymbol{\alpha}</math>, and an observation drawn from a [[multinomial distribution]] with probability vector '''p''' and number of trials ''n''. The Dirichlet parameter vector captures the prior belief about the situation and can be seen as a pseudocount: observations of each outcome that occur before the actual data is collected. The compounding corresponds to a [[Pólya urn model|Pólya urn scheme]]. It is frequently encountered in [[Bayesian statistics]], [[machine learning]], [[empirical Bayes methods]] and [[classical statistics]] as an [[overdispersion|overdispersed]] [[multinomial distribution]].
 
It reduces to the [[categorical distribution]] as a special case when ''n''&nbsp;=&nbsp;1.  It also approximates the [[multinomial distribution]] arbitrarily well for large ''α''. The Dirichlet-multinomial is a multivariate extension of the [[beta-binomial distribution]], as the multinomial and Dirichlet distributions are multivariate versions of the [[binomial distribution]] and [[beta distribution]]s, respectively.

==Specification==
===Dirichlet-multinomial as a compound distribution===
The Dirichlet distribution is a [[conjugate distribution]] to the multinomial distribution. This fact leads to an analytically tractable [[compound distribution]].
For a random vector of category counts <math>\mathbf{x}=(x_1,\dots,x_K)</math>, distributed according to a [[multinomial distribution]], the [[marginal distribution]] is obtained by integrating on the distribution for '''p''' which can be thought of as a [[random vector]] following a Dirichlet distribution:
 
:<math>\Pr(\mathbf{x}\mid n,\boldsymbol{\alpha})=\int_{\mathbf{p}}\mathrm{Mult}(\mathbf{x}\mid  n,\mathbf{p})\mathrm{Dir}(\mathbf{p}\mid\boldsymbol{\alpha})\textrm{d}\mathbf{p}</math>

which results in the following explicit formula:

:<math>\Pr(\mathbf{x}\mid n, \boldsymbol{\alpha})=\frac{\Gamma\left(\alpha_0\right)\Gamma\left(n+1\right)}
{\Gamma\left(n+\alpha_0\right)}\prod_{k=1}^K\frac{\Gamma(x_{k}+\alpha_{k})}{\Gamma(\alpha_{k})\Gamma\left(x_{k}+1\right)}</math>

where <math>\alpha_0</math> is defined as the sum <math>\alpha_0 = \sum \alpha_k</math>. Another form for this same compound distribution, written more compactly in terms of the  [[beta function]], ''B'', is as follows:

<math>\Pr(\mathbf{x}\mid n,\boldsymbol{\alpha})=\frac{n B\left(\alpha_0,n\right)}
{\prod_{k:x_k>0} x_k B\left(\alpha_k,x_k \right)} .
</math>

The latter form emphasizes the fact that zero count categories can be ignored in the calculation - a useful fact when the number of categories is very large and [[Sparse Matrix|sparse]] (e.g. word counts in documents).

Observe that the pdf is the Beta-binomial distribution when <math>K=2</math>. It can also be shown that it approaches the multinomial distribution as <math>\alpha_{0}</math> approaches infinity. The parameter <math>\alpha_{0}</math> governs the degree of overdispersion or [[burstiness]] relative to the multinomial. Alternative choices to denote <math>\alpha_{0}</math> found in the literature are S and A.

===Dirichlet-multinomial as an urn model===
The Dirichlet-multinomial distribution can also be motivated via an [[urn model]] for positive [[integer]] values of the vector <math>\boldsymbol{\alpha}</math>, known as the [[Polya urn model]]. Specifically, imagine an urn containing balls of <math>K</math> colors numbering <math>\alpha_{i}</math> for the ith color, where random draws are made. When a ball is randomly drawn and observed, then two balls of the same color are returned to the urn. If this is performed <math>n</math>  times, then the probability of observing the random vector <math>x</math> of color counts is a Dirichlet-multinomial with parameters <math>n</math> and <math>\boldsymbol\alpha</math>.
If the random draws are with simple replacement (no balls over and above the observed ball are added to the urn), then the distribution follows a multinomial distribution and if the random draws are made without replacement, the distribution follows a [[multivariate hypergeometric distribution]].

==Properties==
===Moments===
Once again, let <math>\alpha_0 = \sum \alpha_k</math> and let <math>p_i =\frac{\alpha_i}{\sum \alpha_k}=\frac{\alpha_i}{\alpha_0}</math>, then the [[Expected value|expected]] number of times the outcome ''i'' was observed over ''n'' trials is

:<math>\operatorname{E}(X_i) = n p_i=n\frac{\alpha_i}{\alpha_0}.\,</math>

The [[covariance matrix]] is as follows.  Each diagonal entry is the [[variance]] of a beta-binomially distributed random variable, and is therefore

:<math>\operatorname{var}(X_i)=np_i(1-p_i)\left(\frac{n+\sum \alpha_k}{1+\sum \alpha_k}\right)=n\frac{\alpha_i}{\alpha_0}\left(1-\frac{\alpha_i}{\alpha_0}\right)\left(\frac{n+\alpha_0}{1+\alpha_0}\right).\,</math>

The off-diagonal entries are the [[covariance]]s:

:<math>\operatorname{cov}(X_i,X_j)=-np_i p_j\left(\frac{n+\sum \alpha_k}{1+\sum \alpha_k}\right)=-n\frac{\alpha_i\alpha_j}{\alpha_0^2}\left(\frac{n+\alpha_0}{1+\alpha_0}\right)\,</math>

for ''i'', ''j'' distinct.

All covariances are negative because for fixed ''n'', an increase in one component of a Dirichlet-multinomial vector requires a decrease in another component.

This is a ''K'' &times; ''K'' [[Positive-definite matrix#Negative-definite, semidefinite and indefinite matrices|positive-semidefinite]] matrix of [[Rank (linear algebra)|rank]] ''K''&nbsp;&minus;&nbsp;1.

The entries of the corresponding [[Correlation matrix#Correlation matrices|correlation matrix]] are

:<math>\rho(X_i,X_i) = 1.</math>

:<math>\rho(X_i,X_j) = \frac{\operatorname{cov}(X_i,X_j)}{\sqrt{\operatorname{var}(X_i)\operatorname{var}(X_j)}} = \frac{-p_i  p_j(\frac{n+\alpha_0}{1+\alpha_0})}{\sqrt{p_i(1-p_i)(\frac{n+\alpha_0}{1+\alpha_0}) p_j(1-p_j)(\frac{n+\alpha_0}{1+\alpha_0})}} =  -\sqrt{\frac{\alpha_i  \alpha_j}{(\alpha_0-\alpha_i)(\alpha_0-\alpha_j)}}.</math>

The sample size drops out of this expression.

Each of the ''k'' components separately has a beta-binomial distribution.

The [[Support (mathematics)|support]] of the Dirichlet-multinomial distribution is the set

: <math>\{(n_1,\dots,n_k)\in \mathbb{N}^{k}| n_1+\cdots+n_k=n\}.\,</math>

Its number of elements is

: <math>{n+k-1 \choose k-1}.</math>

=== Matrix notation ===

In matrix notation, 
:<math>\operatorname{E}(\mathbf{X}) = n \mathbf{p},\,</math>

and 
:<math>\operatorname{var}(\mathbf{X}) = n \lbrace \operatorname{diag}(\mathbf{p}) - \mathbf{p}\mathbf{p}^{\rm T} \rbrace \left( \frac{n+\alpha_0}{1+ \alpha_0} \right) ,\,</math>

with {{math|'''p'''<sup>T</sup>}} = the row vector transpose of the column vector {{math|'''p'''}}. Letting

:<math>\alpha_0 =  \frac{1-\rho^2}{\rho^2}\,</math>, we can write alternatively

:<math>\operatorname{var}(\mathbf{X}) = n \lbrace \operatorname{diag}(\mathbf{p}) - \mathbf{p}\mathbf{p}^{\rm T} \rbrace (1+\rho^2(n-1)) ,\,</math>

The parameter <math> \rho \!</math> is known as the "intra class" or "intra cluster" correlation. It is this positive correlation which gives rise to overdispersion relative to the multinomial distribution.

===Aggregation===
If

:<math>X = (X_1, \ldots, X_K)\sim\operatorname{DM}(\alpha_1,\cdots,\alpha_K)</math>

then, if the random variables with subscripts ''i'' and ''j'' are dropped from the vector and replaced by their sum{{citation needed|reason=this is a nontrivial calculation|date=December 2021}},

:<math>X' = (X_1, \ldots, X_i + X_j, \ldots, X_K)\sim\operatorname{DM} \left(\alpha_1,\cdots,\alpha_i+\alpha_j,\cdots,\alpha_K \right).</math>

This aggregation property may be used to derive the marginal distribution of <math>X_i</math>.

==Likelihood function==
Conceptually, we are making ''N'' independent draws from a categorical distribution with ''K'' categories.  Let us represent the independent draws as random categorical variables <math>z_n</math> for <math>n = 1 \dots N</math>.  Let us denote the number of times a particular category <math>k</math> has been seen (for <math>k = 1 \dots K</math>) among all the categorical variables as <math>n_k</math>, and <math>\sum_k n_k = N</math>.  Then, we have two separate views onto this problem:
# A set of <math>N</math>  categorical variables <math>z_1,\dots,z_N</math>.
# A single vector-valued variable <math>\mathbf{x}=(n_1,\dots,n_K)</math>, distributed according to a [[multinomial distribution]].
The former case is a set of random variables specifying each ''individual'' outcome, while the latter is a variable specifying the ''number'' of outcomes of each of the ''K'' categories.  The distinction is important, as the two cases have correspondingly different probability distributions.

The parameter of the categorical distribution is <math>\mathbf{p} = (p_1,p_2,\dots,p_K),</math> where <math>p_k</math> is the probability to draw value  <math>k</math>;  <math>\mathbf{p}</math>  is likewise the parameter of the multinomial distribution <math>P(\mathbf{x}|\mathbf{p})</math>.  Rather than specifying  <math>\mathbf{p}</math> directly, we give it a [[conjugate prior distribution]], and hence it is drawn from a Dirichlet distribution with parameter vector <math>\boldsymbol\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_K)</math>.

By integrating out <math>\mathbf{p}</math>, we obtain a compound distribution.  However, the form of the distribution is different depending on which view we take.

===For a set of individual outcomes===
====Joint distribution====
For categorical variables <math>\mathbb{Z}=z_1,\dots,z_N</math>, the [[marginal distribution|marginal]] [[joint distribution]] is obtained by integrating out <math>\mathbf{p}</math>:

:<math>\Pr(\mathbb{Z}\mid\boldsymbol{\alpha})=\int_{\mathbf{p}}\Pr(\mathbb{Z}\mid \mathbf{p})\Pr(\mathbf{p}\mid\boldsymbol{\alpha})\textrm{d}\mathbf{p}</math>

which results in the following explicit formula:

:<math>\Pr(\mathbb{Z}\mid\boldsymbol{\alpha})=\frac{\Gamma\left(\alpha_{0}\right)}
{\Gamma\left(N+\alpha_{0}\right)}\prod_{k=1}^K\frac{\Gamma(n_{k}+\alpha_{k})}{\Gamma(\alpha_{k})}</math>

where <math>\Gamma</math> is the [[gamma function]], with

:<math>\alpha_0=\sum_k \alpha_k\text{ and }N=\sum_k n_k\text{, and where }n_k=\text{number of }z_n\text{'s with the value }k.</math>
Note the absence of the multinomial coefficient due to the formula being about the probability of a sequence of categorical variables instead of a probability on the counts within each category.

Although the variables <math>z_1,\dots,z_N</math> do not appear explicitly in the above formula, they enter in through the <math>n_k</math> values.{{clarification needed|reason=By the same reasoning, the z_1, \dots, z_N also enter into the multinomial distribution too. The multinomial distribution is literally the distribution of the vector of the n_k's. So either the purported distinction between Dirichlet-Multinomial and Multinomial does not exist, or the given reason for the distinction is incorrect.|date=June 2021}}

====Conditional distribution====
Another useful formula, particularly in the context of [[Gibbs sampling]], asks what the conditional density of a given variable <math>z_n</math> is, conditioned on all the other variables (which we will denote <math>\mathbb{Z}^{(-n)}</math>).  It turns out to have an extremely simple form:

:<math>\Pr(z_n=k\mid\mathbb{Z}^{(-n)},\boldsymbol{\alpha}) \propto n_k^{(-n)} + \alpha_k</math>

where <math>n_k^{(-n)}</math> specifies the number of counts of category <math>k</math> seen in all variables other than <math>z_n</math>.

It may be useful to show how to derive this formula.  In general, [[conditional distribution]]s are proportional to the corresponding [[joint distribution]]s, so we simply start with the above formula for the joint distribution of all the <math>z_1,\dots,z_N</math> values and then eliminate any factors not dependent on the particular <math>z_n</math> in question.  To do this, we make use of the notation <math>n_k^{(-n)}</math> defined above, and

:<math>
n_j=
\begin{cases}
 n_j^{(-n)}, & \text{if }j\not=k \\
 n_j^{(-n)}+1, & \text{if }j=k
\end{cases}
</math>

We also use the fact that
:<math>\Gamma(n+1) = n\Gamma(n)</math>

Then:

:<math>
\begin{align}
& \Pr(z_n=k\mid\mathbb{Z}^{(-n)},\boldsymbol{\alpha})\\ 
\propto\ & \Pr(z_n=k,\mathbb{Z}^{(-n)}\mid\boldsymbol{\alpha}) \\
=\ &\ \frac{\Gamma\left(\alpha_{0}\right)}{\Gamma\left(N+\alpha_{0}\right)}\prod_{j=1}^K\frac{\Gamma(n_{j}+\alpha_{j})}{\Gamma(\alpha_{j})} \\
\propto\ & \prod_{j=1}^K\Gamma(n_{j}+\alpha_{j}) \\
=\ & \Gamma(n_{k}+\alpha_{k})\prod_{j\not=k}\Gamma(n_{j}+\alpha_{j}) \\
=\ & \Gamma(n_k^{(-n)}+1+\alpha_{k})\prod_{j\not=k}\Gamma(n_j^{(-n)}+\alpha_{j}) \\
=\ & (n_k^{(-n)}+\alpha_{k}) \Gamma(n_k^{(-n)}+\alpha_{k})\prod_{j\not=k}\Gamma(n_j^{(-n)}+\alpha_{j}) \\
=\ & (n_k^{(-n)}+\alpha_{k}) \prod_{j}\Gamma(n_j^{(-n)}+\alpha_{j}) \\
\propto\ & n_k^{(-n)}+\alpha_{k}\\
\end{align}
</math>

In general, it is not necessary to worry about the [[normalizing constant]] at the time of deriving the equations for conditional distributions.  The normalizing constant will be determined as part of the algorithm for sampling from the distribution (see [[Categorical distribution#Sampling]]).  However, when the conditional distribution is written in the simple form above, it turns out that the normalizing constant assumes a simple form:

:<math>\sum_k \left( n_k^{(-n)} + \alpha_{k} \right) = \alpha_{0} + \sum_k n_k^{(-n)} = \alpha_{0} + N - 1</math>

Hence

:<math>\Pr(z_n=k\mid\mathbb{Z}^{(-n)},\boldsymbol{\alpha}) = \frac{n_k^{(-n)} + \alpha_{k}}{\alpha_{0} + N - 1}</math>

This formula is closely related to the [[Chinese restaurant process]], which results from taking the limit as <math>K \to \infty</math>.

====In a Bayesian network====
In a larger [[Bayesian network]] in which categorical (or so-called "multinomial") distributions occur with [[Dirichlet distribution]] priors as part of a larger network, all Dirichlet priors can be collapsed provided that the only nodes depending on them are categorical distributions.  The collapsing happens for each Dirichlet-distribution node separately from the others, and occurs regardless of any other nodes that may depend on the categorical distributions.  It also occurs regardless of whether the categorical distributions depend on nodes additional to the Dirichlet priors (although in such a case, those other nodes must remain as additional conditioning factors). Essentially, all of the categorical distributions depending on a given Dirichlet-distribution node become connected into a single Dirichlet-multinomial joint distribution defined by the above formula.  The joint distribution as defined this way will depend on the parent(s) of the integrated-out Dirichet prior nodes, as well as any parent(s) of the categorical nodes other than the Dirichlet prior nodes themselves.

In the following sections, we discuss different configurations commonly found in Bayesian networks.  We repeat the probability density from above, and define it using the symbol <math>\operatorname{DirMult}(\mathbb{Z}\mid\boldsymbol{\alpha})</math>:

:<math>\Pr(\mathbb{Z}\mid\boldsymbol{\alpha})=\operatorname{DirMult}(\mathbb{Z}\mid\boldsymbol{\alpha})=\frac{\Gamma\left(\sum_k \alpha_k\right)}
{\Gamma\left(\sum_k n_k+\alpha_k\right)}\prod_{k=1}^K\frac{\Gamma(n_{k}+\alpha_{k})}{\Gamma(\alpha_{k})}</math>

=====Multiple Dirichlet priors with the same hyperprior=====

Imagine we have a hierarchical model as follows:

:<math>
\begin{array}{lcl}
\boldsymbol\alpha &\sim& \text{some distribution} \\
\boldsymbol\theta_{d=1 \dots M} &\sim& \operatorname{Dirichlet}_K(\boldsymbol\alpha) \\
z_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{Categorical}_K(\boldsymbol\theta_d)
\end{array}
</math>

In cases like this, we have multiple Dirichet priors, each of which generates some number of categorical observations (possibly a different number for each prior).  The fact that they are all dependent on the same hyperprior, even if this is a random variable as above, makes no difference.  The effect of integrating out a Dirichlet prior links the categorical variables attached to that prior, whose joint distribution simply inherits any conditioning factors of the Dirichlet prior.  The fact that multiple priors may share a hyperprior makes no difference:

:<math>\Pr(\mathbb{Z}\mid\boldsymbol\alpha) = \prod_d \operatorname{DirMult}(\mathbb{Z}_d\mid\boldsymbol\alpha)</math>

where <math>\mathbb{Z}_d</math> is simply the collection of categorical variables dependent on prior ''d''.

Accordingly, the conditional probability distribution can be written as follows:

:<math>\Pr(z_{dn}=k\mid\mathbb{Z}^{(-dn)},\boldsymbol\alpha)\ \propto\ n_{k,d}^{(-n)} + \alpha_k</math>

where <math>n_{k,d}^{(-n)}</math> specifically means the number of variables ''among the set'' <math>\mathbb{Z}_d</math>, excluding <math>z_{dn}</math> itself, that have the value <math>k</math> .

It is necessary to count '''only''' the variables having the value ''k'' that are tied together to the variable in question through having the same prior.  We do '''not''' want to count any other variables also having the value ''k''.

=====Multiple Dirichlet priors with the same hyperprior, with dependent children=====

Now imagine a slightly more complicated hierarchical model as follows:

:<math>
\begin{array}{lcl}
\boldsymbol\alpha &\sim& \text{some distribution} \\
\boldsymbol\theta_{d=1 \dots M} &\sim& \operatorname{Dirichlet}_K(\boldsymbol\alpha) \\
z_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{Categorical}_K(\boldsymbol\theta_d) \\
\boldsymbol\phi &\sim& \text{some other distribution} \\
w_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{F}(w_{dn}\mid z_{dn},\boldsymbol\phi)
\end{array}
</math>

This model is the same as above, but in addition, each of the categorical variables has a child variable dependent on it.  This is typical of a [[mixture model]].

Again, in the joint distribution, only the categorical variables dependent on the same prior are linked into a single Dirichlet-multinomial:

:<math>\Pr(\mathbb{Z},\mathbb{W}\mid\boldsymbol\alpha,\boldsymbol\phi) = \prod_d \operatorname{DirMult}(\mathbb{Z}_d\mid\boldsymbol\alpha) \prod_{d=1}^{M} \prod_{n=1}^{N_d} \operatorname{F}(w_{dn}\mid z_{dn},\boldsymbol\phi)</math>

The conditional distribution of the categorical variables dependent only on their parents and ancestors would have the identical form as above in the simpler case.  However, in Gibbs sampling it is necessary to determine the conditional distribution of a given node <math>z_{dn}</math> dependent not only on <math>\mathbb{Z}^{(-dn)}</math> and ancestors such as <math>\alpha</math> but on ''all'' the other parameters.

The simplified expression for the conditional distribution is derived above simply by rewriting the expression for the joint probability and removing constant factors.  Hence, the same simplification would apply in a larger joint probability expression such as the one in this model, composed of Dirichlet-multinomial densities plus factors for many other random variables dependent on the values of the categorical variables.

This yields the following:

:<math>\Pr(z_{dn}=k\mid\mathbb{Z}^{(-dn)},\mathbb{W},\boldsymbol\alpha,\boldsymbol\phi)\ \propto\ (n_{k,d}^{(-n)} + \alpha_k) \operatorname{F}(w_{dn}\mid z_{dn},\boldsymbol\phi)</math>

Here the probability density of <math>\operatorname{F}</math> appears directly.  To do [[pseudo-random number sampling|random sampling]] over <math>z_{dn}</math>, we would compute the unnormalized probabilities for all ''K'' possibilities for <math>z_{dn}</math> using the above formula, then normalize them and proceed as normal using the algorithm described in the [[categorical distribution]] article.

Correctly speaking, the additional factor that appears in the conditional distribution is derived not from the model specification but directly from the joint distribution.  This distinction is important when considering models where a given node with Dirichlet-prior parent has multiple dependent children, particularly when those children are dependent on each other (e.g. if they share a parent that is collapsed out).  This is discussed more below.

=====Multiple Dirichlet priors with shifting prior membership=====

Now imagine we have a hierarchical model as follows:

:<math>
\begin{array}{lcl}
\boldsymbol\theta &\sim& \text{some distribution} \\
z_{n=1 \dots N} &\sim& \operatorname{Categorical}_K(\boldsymbol\theta) \\
\boldsymbol\alpha &\sim& \text{some distribution} \\
\boldsymbol\phi_{k=1 \dots K} &\sim& \operatorname{Dirichlet}_V(\boldsymbol\alpha) \\
w_{n=1 \dots N} &\sim& \operatorname{Categorical}_V(\boldsymbol\phi_{z_{n}}) \\
\end{array}
</math>

Here we have a tricky situation where we have multiple Dirichlet priors as before and a set of dependent categorical variables, but the relationship between the priors and dependent variables isn't fixed, unlike before.  Instead, the choice of which prior to use is dependent on another random categorical variable.  This occurs, for example, in topic models, and indeed the names of the variables above are meant to correspond to those in [[latent Dirichlet allocation]].  In this case, the set <math>\mathbb{W}</math> is a set of words, each of which is drawn from one of <math>K</math> possible topics, where each topic is a Dirichlet prior over a vocabulary of <math>V</math> possible words, specifying the frequency of different words in the topic.  However, the topic membership of a given word isn't fixed; rather, it's determined from a set of [[latent variable]]s <math>\mathbb{Z}</math>.  There is one latent variable per word, a <math>K</math> -dimensional [[categorical variable]] specifying the topic the word belongs to.

In this case, all variables dependent on a given prior are tied together (i.e. [[correlated]]) in a group, as before — specifically, all words belonging to a given topic are linked.  In this case, however, the group membership shifts, in that the words are not fixed to a given topic but the topic depends on the value of a latent variable associated with the word.  However, the definition of the Dirichlet-multinomial density doesn't actually depend on the number of categorical variables in a group (i.e. the number of words in the document generated from a given topic), but only on the counts of how many variables in the group have a given value (i.e. among all the word tokens generated from a given topic, how many of them are a given word).  Hence, we can still write an explicit formula for the joint distribution:

:<math>\Pr(\mathbb{W}\mid\boldsymbol\alpha,\mathbb{Z}) = \prod_{k=1}^K \operatorname{DirMult}(\mathbb{W}_k\mid\mathbb{Z},\boldsymbol\alpha) = \prod_{k=1}^K \left[\frac{\Gamma\left(\sum_v \alpha_v\right)}
{\Gamma\left(\sum_v n_v^{k}+\alpha_v\right)}\prod_{v=1}^V\frac{\Gamma(n_v^{k}+\alpha_{v})}{\Gamma(\alpha_{v})} \right]</math>

Here we use the notation <math>n_v^{k}</math> to denote the number of word tokens whose value is word symbol ''v'' and which belong to topic ''k''.

The conditional distribution still has the same form:

:<math>\Pr(w_n=v\mid\mathbb{W}^{(-n)},\mathbb{Z},\boldsymbol\alpha)\ \propto\ n_v^{k,(-n)} + \alpha_v</math>

Here again, '''only''' the categorical variables for words belonging to a given topic are linked (even though this linking will depend on the assignments of the latent variables), and hence the word counts need to be over '''only''' the words generated by a given topic.  Hence the symbol <math>n_v^{k,(-n)}</math>, which is the count of words tokens having the word symbol ''v'', but '''only''' among those generated by topic ''k'', and excluding the word itself whose distribution is being described.

(The reason why excluding the word itself is necessary, and why it even makes sense at all, is that in a [[Gibbs sampling]] context, we repeatedly resample the values of each random variable, after having run through and sampled all previous variables.  Hence the variable will already have a value, and we need to exclude this existing value from the various counts that we make use of.)

=====A combined example: LDA topic models=====

We now show how to combine some of the above scenarios to demonstrate how to [[Gibbs sampling|Gibbs sample]] a real-world model, specifically a smoothed [[latent Dirichlet allocation]] (LDA) [[topic model]].

The model is as follows:

:<math>
\begin{array}{lcl}
\boldsymbol\alpha &\sim& \text{A Dirichlet hyperprior, either a constant or a random variable} \\
\boldsymbol\beta &\sim& \text{A Dirichlet hyperprior, either a constant or a random variable} \\
\boldsymbol\theta_{d=1 \dots M} &\sim& \operatorname{Dirichlet}_K(\boldsymbol\alpha) \\
\boldsymbol\phi_{k=1 \dots K} &\sim& \operatorname{Dirichlet}_V(\boldsymbol\beta) \\
z_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{Categorical}_K(\boldsymbol\theta_d) \\
w_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{Categorical}_V(\boldsymbol\phi_{z_{dn}}) \\
\end{array}
</math>

Essentially we combine the previous three scenarios: We have categorical variables dependent on multiple priors sharing a hyperprior; we have categorical variables with dependent children (the [[latent variable]] topic identities); and we have categorical variables with shifting membership in multiple priors sharing a hyperprior. In the standard LDA model, the words are completely observed, and hence we never need to resample them. (However, Gibbs sampling would equally be possible if only some or none of the words were observed.  In such a case, we would want to initialize the distribution over the words in some reasonable fashion — e.g. from the output of some process that generates sentences, such as a [[machine translation]] model — in order for the resulting [[posterior distribution|posterior]] latent variable distributions to make any sense.)

Using the above formulas, we can write down the conditional probabilities directly:

:<math>
\begin{array}{lcl}
\Pr(w_{dn}=v\mid\mathbb{W}^{(-dn)},\mathbb{Z},\boldsymbol\beta)\ &\propto\ & \#\mathbb{W}_v^{k,(-dn)} + \beta_v \\
\Pr(z_{dn}=k\mid\mathbb{Z}^{(-dn)},w_{dn}=v,\mathbb{W}^{(-dn)},\boldsymbol\alpha)\ &\propto\ &(\#\mathbb{Z}_k^{d,(-dn)} + \alpha_k) \Pr(w_{dn}=v\mid\mathbb{W}^{(-dn)},\mathbb{Z},\boldsymbol\beta) \\
\end{array}
</math>

Here we have defined the counts more explicitly to clearly separate counts of words and counts of topics:

:<math>
\begin{array}{lcl}
\#\mathbb{W}_v^{k,(-dn)} &=& \text{number of words having value }v\text{ among topic }k\text{ excluding }w_{dn} \\
\#\mathbb{Z}_k^{d,(-dn)} &=& \text{number of topics having value }k\text{ among document }d\text{ excluding }z_{dn} \\
\end{array}
</math>

As in the scenario above with categorical variables with dependent children, the conditional probability of those dependent children appears in the definition of the parent's conditional probability.  In this case, each latent variable has only a single dependent child word, so only one such term appears. (If there were multiple dependent children, all would have to appear in the parent's conditional probability, regardless of whether there was overlap between different parents and the same children, i.e. regardless of whether the dependent children of a given parent also have other parents.  In a case where a child has multiple parents, the conditional probability for that child appears in the conditional probability definition of each of its parents.)

The definition above specifies only the ''unnormalized'' conditional probability of the words, while the topic conditional probability requires the ''actual'' (i.e. normalized) probability.  Hence we have to normalize by summing over all word symbols:

:<math>
\begin{array}{rcl}
\Pr(z_{dn}=k\mid\mathbb{Z}^{(-dn)},w_{dn}=v,\mathbb{W}^{(-dn)},\boldsymbol\alpha)\ &\propto\ &\bigl(\#\mathbb{Z}_k^{d,(-dn)} + \alpha_k\bigr) \dfrac{\#\mathbb{W}_v^{k,(-dn)} + \beta_v}{\sum_{v'=1}^{V} (\#\mathbb{W}_{v'}^{k,(-dn)} + \beta_{v'})} \\
&& \\
&=& \bigl(\#\mathbb{Z}_k^{d,(-dn)} + \alpha_k\bigr) \dfrac{\#\mathbb{W}_v^{k,(-dn)} + \beta_v}{\#\mathbb{W}^{k} + B - 1}
\end{array}
</math>

where

:<math>
\begin{array}{lcl}
\#\mathbb{W}^{k} &=& \text{number of words generated by topic }k \\
B &=& \sum_{v=1}^{V} \beta_v \\
\end{array}
</math>

It's also worth making another point in detail, which concerns the second factor above in the conditional probability.  Remember that the conditional distribution in general is derived from the joint distribution, and simplified by removing terms not dependent on the domain of the conditional (the part on the left side of the vertical bar).  When a node <math>z</math> has dependent children, there will be one or more factors <math>\operatorname{F}(\dots\mid z)</math> in the joint distribution that are dependent on <math>z</math>.  ''Usually'' there is one factor for each dependent node, and it has the same density function as the distribution appearing the mathematical definition.  However, if a dependent node has another parent as well (a co-parent), and that co-parent is collapsed out, then the node will become dependent on all other nodes sharing that co-parent, and in place of multiple terms for each such node, the joint distribution will have only one joint term.  We have exactly that situation here.  Even though <math>z_{dn}</math> has only one child <math>w_{dn}</math>, that child has a Dirichlet co-parent that we have collapsed out, which induces a Dirichlet-multinomial over the entire set of nodes <math>\mathbb{W}^{k}</math>.

It happens in this case that this issue does not cause major problems, precisely because of the one-to-one relationship between <math>z_{dn}</math> and <math>w_{dn}</math>.  We can rewrite the joint distribution as follows:

:<math>
\begin{array}{lcl}
p(\mathbb{W}^{k}\mid z_{dn}) &=& p(w_{dn}\mid\mathbb{W}^{k,(-dn)},z_{dn})\,p(\mathbb{W}^{k,(-dn)}\mid z_{dn}) \\
&=& p(w_{dn}\mid\mathbb{W}^{k,(-dn)},z_{dn})\,p(\mathbb{W}^{k,(-dn)}) \\
&\sim& p(w_{dn}\mid\mathbb{W}^{k,(-dn)},z_{dn})
\end{array}
</math>

where in the set <math>\mathbb{W}^{k,(-dn)}</math> (i.e. the set of nodes <math>\mathbb{W}^{k}</math> excluding <math>w_{dn}</math> ), none of the nodes have <math>z_{dn}</math> as a parent.  Hence it can be eliminated as a conditioning factor (line 2), meaning that the entire factor can be eliminated from the conditional distribution (line 3).

=====A second example: Naive Bayes document clustering=====

Here is another model, with a different set of issues.  This is an implementation of an unsupervised [[Naive Bayes]] model for document clustering.  That is, we would like to [[document classification|classify documents]] into multiple categories (e.g. "[[spam (electronic)|spam]]" or "non-spam", or "scientific journal article", "newspaper article about finance", "newspaper article about politics", "love letter") based on textual content.  However, we don't already know the correct category of any documents; instead, we want to [[document clustering|cluster]] them based on mutual similarities. (For example, a set of scientific articles will tend to be similar to each other in word use but very different from a set of love letters.) This is a type of [[unsupervised learning]]. (The same technique can be used for doing [[semi-supervised learning]], i.e. where we know the correct category of some fraction of the documents and would like to use this knowledge to help in clustering the remaining documents.)

The model is as follows:

:<math>
\begin{array}{lcl}
\boldsymbol\alpha &\sim& \text{A Dirichlet hyperprior, either a constant or a random variable} \\
\boldsymbol\beta &\sim& \text{A Dirichlet hyperprior, either a constant or a random variable} \\
\boldsymbol\theta_{d=1 \dots M} &\sim& \operatorname{Dirichlet}_K(\boldsymbol\alpha) \\
\boldsymbol\phi_{k=1 \dots K} &\sim& \operatorname{Dirichlet}_V(\boldsymbol\beta) \\
z_{d=1 \dots M} &\sim& \operatorname{Categorical}_K(\boldsymbol\theta_d) \\
w_{d=1 \dots M,n=1 \dots N_d} &\sim& \operatorname{Categorical}_V(\boldsymbol\phi_{z_{d}}) \\
\end{array}
</math>

In many ways, this model is very similar to the [[latent Dirichlet allocation|LDA]] [[topic model]] described above, but it assumes one topic per document rather than one topic per word, with a document consisting of a mixture of topics.  This can be seen clearly in the above model, which is identical to the LDA model except that there is only one [[latent variable]] per document instead of one per word.  Once again, we assume that we are collapsing all of the Dirichlet priors.

The conditional probability for a given word is almost identical to the LDA case.  Once again, all words generated by the same Dirichlet prior are interdependent.  In this case, this means the words of all documents having a given label — again, this can vary depending on the label assignments, but all we care about is the total counts.  Hence:

:<math>
\begin{array}{lcl}
\Pr(w_{dn}=v\mid\mathbb{W}^{(-dn)},\mathbb{Z},\boldsymbol\beta)\ &\propto\ & \#\mathbb{W}_v^{k,(-dn)} + \beta_v \\
\end{array}
</math>

where

:<math>
\begin{array}{lcl}
\#\mathbb{W}_v^{k,(-dn)} &=& \text{number of words having value }v\text{ among documents with label }k\text{ excluding }w_{dn} \\
\end{array}
</math>

However, there is a critical difference in the conditional distribution of the latent variables for the label assignments, which is that a given label variable has multiple children nodes instead of just one — in particular, the nodes for all the words in the label's document.  This relates closely to the discussion above about the factor <math>\operatorname{F}(\dots\mid z_d)</math> that stems from the joint distribution.  In this case, the joint distribution needs to be taken over all words in all documents containing a label assignment equal to the value of <math>z_d</math>, and has the value of a Dirichlet-multinomial distribution.  Furthermore, we cannot reduce this joint distribution down to a conditional distribution over a single word.  Rather, we can reduce it down only to a smaller joint conditional distribution over the words in the document for the label in question, and hence we cannot simplify it using the trick above that yields a simple sum of expected count and prior.  Although it is in fact possible to rewrite it as a product of such individual sums, the number of factors is very large, and is not clearly more efficient than directly computing the Dirichlet-multinomial distribution probability.

==Related distributions==
The one-dimensional version of the Dirichlet-multinomial distribution is known as the [[Beta-binomial distribution]].

The Dirichlet-multinomial distribution has a relationship with the [[negative binomial]] distribution analogous to the relationship of the [[multinomial distribution]] with the [[Poisson distribution]].<ref name=Zhou2018>Theorem 1 of {{cite journal |last1=Zhou |first=M.|year=2018|title=Nonparametric Bayesian Negative Binomial Factor Analysis |journal=Bayesian Analysis |volume=13 |issue=4|pages=1065–1093|doi=10.1214/17-BA1070 |doi-access=free |arxiv=1604.07464 }}</ref>

==Uses==
The Dirichlet-multinomial distribution is used in automated [[document classification]] and clustering, [[genetics]], [[economy]], combat modeling, and quantitative marketing.

{{more footnotes|date=June 2012}}

==See also==
* [[Beta-binomial distribution]]
* [[Chinese restaurant process]]
* [[Stars and bars (combinatorics)]]
* [[Dirichlet process]]
* [[Generalized Dirichlet distribution]]
* [[Krichevsky–Trofimov estimator]]
* [[Dirichlet negative multinomial distribution]]

== References ==
=== Citations ===
{{Reflist}}

=== Sources ===
{{refbegin}}
* Elkan, C. (2006) [http://www-cse.ucsd.edu/users/elkan/edcm.pdf Clustering documents with an exponential-family approximation of the Dirichlet compound multinomial distribution]. ICML, 289–296.
* Johnson, N. L., Kotz, S. and Balakrishnan, N. (1997) Discrete multivariate distributions (Vol. 165). New York: Wiley.
* Kvam, P. and Day, D. (2001) The multivariate Polya distribution in combat modeling. Naval Research Logistics, 48, 1–17.
* Madsen, R. E., Kauchak, D. and Elkan, C. (2005) [http://www.cse.ucsd.edu/~dkauchak/kauchak05modeling.pdf  Modeling Word Burstiness Using the Dirichlet Distribution]. ICML, 545–552.
* Minka, T. (2003) [http://research.microsoft.com/~minka/papers/dirichlet/ Estimating a Dirichlet distribution]. Technical report Microsoft Research. Includes Matlab code for fitting distributions to data.
* Mosimann, J. E. (1962) [https://www.jstor.org/stable/2333468 On the compound multinomial distribution, the multivariate β-distribution, and correlations among proportions]. Biometrika, 49(1–2), 65–82.
* Wagner, U. and Taudes, A. (1986) A Multivariate Polya Model of Brand Choice and Purchase Incidence. Marketing Science, 5(3), 219–244.
{{refend}}

{{-}}
{{ProbDistributions|multivariate}}
{{Peter Gustav Lejeune Dirichlet}}

[[Category:Multivariate discrete distributions]]
[[Category:Discrete distributions]]
[[Category:Compound probability distributions]]