Multinomial distribution

Properties

Normalization

The multinomial distribution is normalized according to:

$${\displaystyle \sum _{\sum _{j=1}^k{x}_j=n}f(x_1,\dots ,x_k;n,p_1,\dots ,p_k)=1}$$

where the sum is over all permutations of ${\displaystyle x_j}$ such that $\sum _{j=1}^k{x}_j=n$.

Expected value and variance

The expected number of times the outcome i was observed over n trials is

$${\displaystyle \mathrm{E}(X_i)=n{p}_i.}$$

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

$${\displaystyle \mathrm{Var}(X_i)=n{p}_i(1-p_i).}$$

The off-diagonal entries are the covariances:

$${\displaystyle \mathrm{Cov}(X_i,X_j)=-n{p}_i{p}_j}$$

for i, j distinct.

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

When these expressions are combined into a matrix with i, j element ${\displaystyle \mathrm{cov}(X_i,X_j),}$ the result is a k × k positive-semidefinite covariance matrix of rank k − 1. In the special case where k = n and where the p_i are all equal, the covariance matrix is the centering matrix.

The entries of the corresponding correlation matrix are

$${\displaystyle \begin{array}{cc}\rho (X_i,X_i)& =1\\ \rho (X_i,X_j)& =\frac{\mathrm{Cov}(X_i,X_j)}{\sqrt{\mathrm{Var}(X_i)\mathrm{Var}(X_j)}}\\ & =\frac{-p_i{p}_j}{\sqrt{p_i(1-p_i)p_j(1-p_j)}}\\ & =-\sqrt{\frac{p_i{p}_j}{(1-p_i)(1-p_j)}}.\end{array}}$$

Note that the number of trials n drops out of this expression.

Each of the k components separately has a binomial distribution with parameters n and p_i, for the appropriate value of the subscript i.

The support of the multinomial distribution is the set

$${\displaystyle \left\{(n_1,\dots ,n_k)\in {\mathbb{N}}^k\mid n_1+\cdots +n_k=n\right\}.}$$

Its number of elements is

$${\displaystyle (\genfrac{}{}{0ex}{}{n+k-1}{k-1}).}$$

Matrix notation

In matrix notation, $${\displaystyle \mathrm{E}(\mathbf{X})=n\mathbf{p},}$$

and $${\displaystyle \mathrm{Var}(\mathbf{X})=n\{\mathrm{diag}(\mathbf{p})-\mathbf{p}{\mathbf{p}}^{\mathrm{T}}\},}$$

with p^T = the row vector transpose of the column vector p.

Visualization

As slices of generalized Pascal's triangle

Just like one can interpret the binomial distribution as (normalized) one-dimensional (1D) slices of Pascal's triangle, so too can one interpret the multinomial distribution as 2D (triangular) slices of Pascal's pyramid, or 3D/4D/+ (pyramid-shaped) slices of higher-dimensional analogs of Pascal's triangle. This reveals an interpretation of the range of the distribution: discretized equilateral "pyramids" in arbitrary dimension—i.e. a simplex with a grid.

As polynomial coefficients

Similarly, just like one can interpret the binomial distribution as the polynomial coefficients of ${\displaystyle (p+q{)}^n}$ when expanded, one can interpret the multinomial distribution as the coefficients of ${\displaystyle (p_1+p_2+p_3+\cdots +p_k{)}^n}$ when expanded, noting that just the coefficients must sum up to 1.

Large deviation theory

Asymptotics

By Stirling's formula, at the limit of ${\displaystyle n,x_1,\dots ,x_k\to \mathrm{\infty }}$, we have$${\displaystyle \ln (\genfrac{}{}{0ex}{}{n}{x_1,\dots ,x_k})+\sum _{i=1}^k{x}_i\ln p_i=-n{D}_{\text{KL}}(\hat{p}\Vert p)-\frac{k-1}{2}\ln (2\pi n)-\frac{1}{2}\sum _{i=1}^k\ln ({\hat{p}}_i)+o(1)}$$where relative frequencies ${\displaystyle {\hat{p}}_i=x_i/n}$ in the data can be interpreted as probabilities from the empirical distribution ${\displaystyle \hat{p}}$, and ${\displaystyle D_{\text{KL}}}$ is the Kullback–Leibler divergence.

This formula can be interpreted as follows.

Consider ${\displaystyle {\mathrm{\Delta }}_k}$, the space of all possible distributions over the categories ${\displaystyle \{1,2,\dots ,k\}}$. It is a simplex. After ${\displaystyle n}$ independent samples from the categorical distribution ${\displaystyle p}$ (which is how we construct the multinomial distribution), we obtain an empirical distribution ${\displaystyle \hat{p}}$.

By the asymptotic formula, the probability that empirical distribution ${\displaystyle \hat{p}}$ deviates from the actual distribution ${\displaystyle p}$ decays exponentially as we sample more data, at a rate of ${\displaystyle D_{\text{KL}}(\hat{p}\Vert p)}$. The more experiments and the more different ${\displaystyle \hat{p}}$ is from ${\displaystyle p}$, the less likely it is to see such an empirical distribution.

If ${\displaystyle A}$ is a closed subset of ${\displaystyle {\mathrm{\Delta }}_k}$, then by dividing up ${\displaystyle A}$ into pieces, and reasoning about the growth rate of ${\displaystyle Pr(\hat{p}\in A_{ϵ})}$ on each piece ${\displaystyle A_{ϵ}}$, we obtain Sanov's theorem, which states that $${\displaystyle \underset{n\to \mathrm{\infty }}{lim}\frac{1}{n}\ln Pr(\hat{p}\in A)=-\underset{\hat{p}\in A}{inf}{D}_{\text{KL}}(\hat{p}\Vert p)}$$

Concentration at large n

Due to the exponential decay, at large ${\displaystyle n}$, almost all the probability mass is concentrated in a small neighborhood of ${\displaystyle p}$. In this small neighborhood, we can take the first nonzero term in the Taylor expansion of ${\displaystyle D_{KL}}$, to obtain$${\displaystyle \begin{array}{cc}\ln (\genfrac{}{}{0ex}{}{n}{x_1,\cdots ,x_k})p_1^{x_1}\cdots p_k^{x_k}& \approx -\frac{n}{2}\sum _{i=1}^k\frac{({\hat{p}}_i-p_i{)}^2}{p_i}\\ & =-\frac{1}{2}\sum _{i=1}^k\frac{(x_i-n{p}_i{)}^2}{n{p}_i}\end{array}}$$ This resembles the Gaussian distribution, which suggests the following theorem:

Theorem. At the ${\displaystyle n\to \mathrm{\infty }}$ limit, ${\displaystyle n\sum _{i=1}^k\frac{({\hat{p}}_i-p_i{)}^2}{p_i}=\sum _{i=1}^k\frac{(x_i-n{p}_i{)}^2}{n{p}_i}}$ converges in distribution to the chi-squared distribution ${\displaystyle {\chi }^2(k-1)}$.

Conditional concentration at large n

The above concentration phenomenon can be easily generalized to the case where we condition upon independent constraints. This is the theoretical justification for Pearson's chi-squared test.

Theorem.

• Given functions ${\displaystyle f_1,\dots ,f_{\ell }}$, such that they are continuously differentiable in a neighborhood of ${\displaystyle p}$, and the vectors ${\displaystyle (1,1,\dots ,1),\mathrm{\nabla }{f}_1(p),\dots ,\mathrm{\nabla }{f}_{\ell }(p)}$ are linearly independent;
• given sequences ${\displaystyle {ϵ}_1(n),\dots ,{ϵ}_{\ell }(n)}$, such that asymptotically ${\displaystyle \frac{1}{n}\ll {ϵ}_i(n)\ll \frac{1}{\sqrt{n}}}$ for each ${\displaystyle i\in \{1,\dots ,\ell \}}$;
• then for the multinomial distribution conditional on constraints ${\displaystyle f_1(\hat{p})\in [f_1(p)-{ϵ}_1(n),f_1(p)+{ϵ}_1(n)],\dots ,f_{\ell }(\hat{p})\in [f_{\ell }(p)-{ϵ}_{\ell }(n),f_{\ell }(p)+{ϵ}_{\ell }(n)]}$, we have the quantity ${\displaystyle n\sum _i\frac{({\hat{p}}_i-p_i{)}^2}{p_i}=\sum _i\frac{(x_i-n{p}_i{)}^2}{n{p}_i}}$ converging in distribution to ${\displaystyle {\chi }^2(k-1-\ell )}$ at the ${\displaystyle n\to \mathrm{\infty }}$ limit.

In the case that all ${\displaystyle {\hat{p}}_i}$ are equal, this reduces to the concentration of entropies around the maximum entropy.

This theorem can be shown by starting with the previous case, then taking the conditional on the constraints.

Related distributions

In some fields such as natural language processing, categorical and multinomial distributions are synonymous and it is common to speak of a multinomial distribution when a categorical distribution is actually meant. This stems from the fact that it is sometimes convenient to express the outcome of a categorical distribution as a "1-of-k" vector (a vector with one element containing a 1 and all other elements containing a 0) rather than as an integer in the range ${\displaystyle 1\dots k}$; in this form, a categorical distribution is equivalent to a multinomial distribution over a single trial.

• When k = 2, the multinomial distribution is the binomial distribution.
• Categorical distribution, the distribution of each trial; for k = 2, this is the Bernoulli distribution.
• The Dirichlet distribution is the conjugate prior of the multinomial in Bayesian statistics.
• Dirichlet-multinomial distribution.
• Beta-binomial distribution.
• Negative multinomial distribution
• Hardy–Weinberg principle ( a trinomial distribution with probabilities ${\displaystyle ({\theta }^2,2\theta (1-\theta ),(1-\theta {)}^2)}$)

Statistical inference

Equivalence tests for multinomial distributions

The goal of equivalence testing is to establish the agreement between a theoretical multinomial distribution and observed counting frequencies. The theoretical distribution may be a fully specified multinomial distribution or a parametric family of multinomial distributions.

Let ${\displaystyle q}$ denote a theoretical multinomial distribution and let ${\displaystyle p}$ be a true underlying distribution. The distributions ${\displaystyle p}$ and ${\displaystyle q}$ are considered equivalent if ${\displaystyle d(p,q)<\epsilon }$ for a distance ${\displaystyle d}$ and a tolerance parameter ${\displaystyle \epsilon >0}$. The equivalence test problem is ${\displaystyle H_0=\{d(p,q)\ge \epsilon \}}$ versus ${\displaystyle H_1=\{d(p,q)<\epsilon \}}$. The true underlying distribution ${\displaystyle p}$ is unknown. Instead, the counting frequencies ${\displaystyle p_n}$ are observed, where ${\displaystyle n}$ is a sample size. An equivalence test uses ${\displaystyle p_n}$ to reject ${\displaystyle H_0}$. If ${\displaystyle H_0}$ can be rejected then the equivalence between ${\displaystyle p}$ and ${\displaystyle q}$ is shown at a given significance level. The equivalence test for Euclidean distance can be found in text book of Wellek (2010). The equivalence test for the total variation distance is developed in Ostrovski (2017). The exact equivalence test for the specific cumulative distance is proposed in Frey (2009).

The distance between the true underlying distribution ${\displaystyle p}$ and a family of the multinomial distributions ${\displaystyle \mathcal{M}}$ is defined by ${\displaystyle d(p,\mathcal{M})=\underset{h\in \mathcal{M}}{min}d(p,h)}$. Then the equivalence test problem is given by ${\displaystyle H_0=\{d(p,\mathcal{M})\ge \epsilon \}}$ and ${\displaystyle H_1=\{d(p,\mathcal{M})<\epsilon \}}$. The distance ${\displaystyle d(p,\mathcal{M})}$ is usually computed using numerical optimization. The tests for this case are developed recently in Ostrovski (2018).

Confidence intervals for the difference of two proportions

In the setting of a multinomial distribution, constructing confidence intervals for the difference between the proportions of observations from two events, ${\displaystyle p_i-p_j}$, requires the incorporation of the negative covariance between the sample estimators ${\displaystyle {\hat{p}}_i=\frac{X_i}{n}}$ and ${\displaystyle {\hat{p}}_j=\frac{X_j}{n}}$.

Some of the literature on the subject focused on the use-case of matched-pairs binary data, which requires careful attention when translating the formulas to the general case of ${\displaystyle p_i-p_j}$ for any multinomial distribution. Formulas in the current section will be generalized, while formulas in the next section will focus on the matched-pairs binary data use-case.

Wald's standard error (SE) of the difference of proportion can be estimated using:

$${\displaystyle \widehat{\mathrm{SE}}({\hat{p}}_i-{\hat{p}}_j)=\sqrt{\frac{\left({\hat{p}}_i+{\hat{p}}_j\right)-{\left({\hat{p}}_i-{\hat{p}}_j\right)}^2}{n}}}$$

For a ${\displaystyle 100(1-\alpha )\mathrm{\%}}$ approximate confidence interval, the margin of error may incorporate the appropriate quantile from the standard normal distribution, as follows:

$${\displaystyle ({\hat{p}}_i-{\hat{p}}_j)\pm z_{\alpha /2}\cdot \widehat{\mathrm{SE}}({\hat{p}}_i-{\hat{p}}_j)}$$

A modification which includes a continuity correction adds ${\displaystyle \frac{1}{n}}$ to the margin of error as follows:

$${\displaystyle ({\hat{p}}_i-{\hat{p}}_j)\pm \left(z_{\alpha /2}\cdot \widehat{\mathrm{SE}}({\hat{p}}_i-{\hat{p}}_j)+\frac{1}{n}\right)}$$

Another alternative is to rely on a Bayesian estimator using Jeffreys prior which leads to using a dirichlet distribution, with all parameters being equal to 0.5, as a prior. The posterior will be the calculations from above, but after adding 1/2 to each of the k elements, leading to an overall increase of the sample size by ${\displaystyle \frac{k}{2}}$. This was originally developed for a multinomial distribution with four events, and is known as wald+2, for analyzing matched pairs data (see the next section for more details).

This leads to the following SE:

$${\displaystyle \widehat{\mathrm{SE}}{({\hat{p}}_i-{\hat{p}}_j)}_{wald+\frac{k}{2}}=\sqrt{\frac{\left({\hat{p}}_i+{\hat{p}}_j+\frac{1}{n}\right)\frac{n}{n+\frac{k}{2}}-{\left({\hat{p}}_i-{\hat{p}}_j\right)}^2{\left(\frac{n}{n+\frac{k}{2}}\right)}^2}{n+\frac{k}{2}}}}$$

Which can just be plugged into the original Wald formula as follows:

$${\displaystyle \left(p_i-p_j\right)\frac{n}{n+\frac{k}{2}}\pm z_{\alpha /2}\cdot \widehat{\mathrm{SE}}{({\hat{p}}_i-{\hat{p}}_j)}_{wald+\frac{k}{2}}}$$

Occurrence and applications

Confidence intervals for the difference in matched-pairs binary data (using multinomial with k=4)

For the case of matched-pairs binary data, a common task is to build the confidence interval of the difference of the proportion of the matched events. For example, we might have a test for some disease, and we may want to check the results of it for some population at two points in time (1 and 2), to check if there was a change in the proportion of the positives for the disease during that time.

Such scenarios can be represented using a two-by-two contingency table with the number of elements that had each of the combination of events. We can use small f for sampling frequencies: ${\displaystyle f_{11},f_{10},f_{01},f_{00}}$, and capital F for population frequencies: ${\displaystyle F_{11},F_{10},F_{01},F_{00}}$. These four combinations could be modeled as coming from a multinomial distribution (with four potential outcomes). The sizes of the sample and population can be n and N respectively. And in such a case, there is an interest in building a confidence interval for the difference of proportions from the marginals of the following (sampled) contingency table:

In this case, checking the difference in marginal proportions means we are interested in using the following definitions: ${\displaystyle p_{1\ast }=\frac{F_{1\ast }}{N}=\frac{F_{11}+F_{10}}{N}}$, ${\displaystyle p_{\ast 1}=\frac{F_{\ast 1}}{N}=\frac{F_{11}+F_{01}}{N}}$. And the difference we want to build confidence intervals for is:

$${\displaystyle p_{\ast 1}-p_{1\ast }=\frac{F_{11}+F_{01}}{N}-\frac{F_{11}+F_{10}}{N}=\frac{F_{01}}{N}-\frac{F_{10}}{N}=p_{01}-p_{10}}$$

Hence, a confidence intervals for the marginal positive proportions (${\displaystyle p_{\ast 1}-p_{1\ast }}$) is the same as building a confidence interval for the difference of the proportions from the secondary diagonal of the two-by-two contingency table (${\displaystyle p_{01}-p_{10}}$).

Calculating a p-value for such a difference is known as McNemar's test. Building confidence interval around it can be constructed using methods described above for Confidence intervals for the difference of two proportions.

The Wald confidence intervals from the previous section can be applied to this setting, and appears in the literature using alternative notations. Specifically, the SE often presented is based on the contingency table frequencies instead of the sample proportions. For example, the Wald confidence intervals, provided above, can be written as:

$${\displaystyle \begin{array}{cc}\widehat{\mathrm{SE}}(p_{\ast 1}-p_{1\ast })& =\widehat{\mathrm{SE}}(p_{01}-p_{10})\\ & =\frac{\sqrt{n\left(f_{10}+f_{01}\right)-{\left(f_{10}-f_{01}\right)}^2}}{n\sqrt{n}}\end{array}}$$

Further research in the literature has identified several shortcomings in both the Wald and the Wald with continuity correction methods, and other methods have been proposed for practical application.

One such modification includes Agresti and Min’s Wald+2 (similar to some of their other works) in which each cell frequency had an extra ${\displaystyle \frac{1}{2}}$ added to it. This leads to the Wald+2 confidence intervals. In a Bayesian interpretation, this is like building the estimators taking as prior a dirichlet distribution with all parameters being equal to 0.5 (which is, in fact, the Jeffreys prior). The +2 in the name wald+2 can now be taken to mean that in the context of a two-by-two contingency table, which is a multinomial distribution with four possible events, then since we add 1/2 an observation to each of them, then this translates to an overall addition of 2 observations (due to the prior).

This leads to the following modified SE for the case of matched pairs data:

$${\displaystyle \widehat{\mathrm{SE}}(p_{\ast 1}-p_{1\ast })=\frac{\sqrt{\left(n+2\right)\left(f_{10}+f_{01}+1\right)-{\left(f_{10}-f_{01}\right)}^2}}{\left(n+2\right)\sqrt{n+2}}}$$

Which can just be plugged into the original Wald formula as follows:

$${\displaystyle \left(p_{\ast 1}-p_{1\ast }\right)\frac{n}{n+2}\pm z_{\alpha /2}\cdot \widehat{\mathrm{SE}}({\hat{p}}_i-{\hat{p}}_j{)}_{wald+2}}$$

Other modifications include Bonett and Price’s Adjusted Wald, and Newcombe’s Score.

Computational methods

Random variate generation

First, reorder the parameters ${\displaystyle p_1,\dots ,p_k}$ such that they are sorted in descending order (this is only to speed up computation and not strictly necessary). Now, for each trial, draw an auxiliary variable X from a uniform (0, 1) distribution. The resulting outcome is the component

$${\displaystyle j=min\left\{j^{\prime }\in \{1,\dots ,k\}:\left(\sum _{i=1}^{j^{\prime }}{p}_i\right)-X\ge 0\right\}.}$$

{X_j = 1, X_k = 0 for k ≠ j } is one observation from the multinomial distribution with ${\displaystyle p_1,\dots ,p_k}$ and n = 1. A sum of independent repetitions of this experiment is an observation from a multinomial distribution with n equal to the number of such repetitions.

Sampling using repeated conditional binomial samples

Given the parameters ${\displaystyle p_1,p_2,\dots ,p_k}$ and a total for the sample ${\displaystyle n}$ such that $\sum _{i=1}^k{X}_i=n$, it is possible to sample sequentially for the number in an arbitrary state ${\displaystyle X_i}$, by partitioning the state space into ${\displaystyle i}$ and not-${\displaystyle i}$, conditioned on any prior samples already taken, repeatedly.

Algorithm: Sequential conditional binomial sampling

<span id="mwAkg"></span><span class="nv" id="mwAkk">S</span><span class="w" id="mwAko"> </span><span class="o" id="mwAks">=</span><span class="w" id="mwAkw"> </span>n <span class="nv" id="mwAk0">rho</span><span class="w" id="mwAk4"> </span><span class="o" id="mwAk8">=</span><span class="w" id="mwAlA"> </span><span class="m" id="mwAlE">1</span> <span class="k" id="mwAlI">for</span><span class="w" id="mwAlM"> </span>i<span class="w" id="mwAlQ"> </span><span class="k" id="mwAlU">in</span><span class="w" id="mwAlY"> </span><span class="o" id="mwAlc">[</span><span class="m" id="mwAlg">1</span>,k-1<span class="o" id="mwAlk">]</span>: <span class="w" id="mwAlo"> </span><span class="k" id="mwAls">if</span><span class="w" id="mwAlw"> </span>rho<span class="w" id="mwAl0"> </span>!<span class="o" id="mwAl4">=</span><span class="w" id="mwAl8"> </span><span class="m" id="mwAmA">0</span>: <span class="w" id="mwAmE"> </span>X<span class="o" id="mwAmI">[</span>i<span class="o" id="mwAmM">]</span><span class="w" id="mwAmQ"> </span>~<span class="w" id="mwAmU"> </span>Binom<span class="o" id="mwAmY">(</span>S,p<span class="o" id="mwAmc">[</span>i<span class="o" id="mwAmg">]</span>/rho<span class="o" id="mwAmk">)</span> <span class="w" id="mwAmo"> </span><span class="k" id="mwAms">else</span><span class="w" id="mwAmw"> </span> <span class="w" id="mwAm0"> </span>X<span class="o" id="mwAm4">[</span>i<span class="o" id="mwAm8">]</span><span class="w" id="mwAnA"> </span><span class="o" id="mwAnE">=</span><span class="w" id="mwAnI"> </span><span class="m" id="mwAnM">0</span> <span class="w" id="mwAnQ"> </span><span class="nv" id="mwAnU">S</span><span class="w" id="mwAnY"> </span><span class="o" id="mwAnc">=</span><span class="w" id="mwAng"> </span>S<span class="w" id="mwAnk"> </span>-<span class="w" id="mwAno"> </span>X<span class="o" id="mwAns">[</span>i<span class="o" id="mwAnw">]</span> <span class="w" id="mwAn0"> </span><span class="nv" id="mwAn4">rho</span><span class="w" id="mwAn8"> </span><span class="o" id="mwAoA">=</span><span class="w" id="mwAoE"> </span>rho<span class="w" id="mwAoI"> </span>-<span class="w" id="mwAoM"> </span>p<span class="o" id="mwAoQ">[</span>i<span class="o" id="mwAoU">]</span> X<span class="o" id="mwAoY">[</span>k<span class="o" id="mwAoc">]</span><span class="w" id="mwAog"> </span><span class="o" id="mwAok">=</span><span class="w" id="mwAoo"> </span>S

Heuristically, each application of the binomial sample reduces the available number to sample from and the conditional probabilities are likewise updated to ensure logical consistency.

Software implementations

• The MultinomialCI R package allows the computation of simultaneous confidence intervals for the probabilities of a multinomial distribution given a set of observations.