Multinomial distribution

The multinomial distribution is a discrete probability distribution describing the numbers of outcomes assigned to several mutually exclusive categories after a fixed number of independent trials. Each trial produces exactly one category, and the probability associated with each category remains constant across trials. It generalizes the binomial distribution, which is recovered when only two categories are present.

For (k) categories, let the probability vector be

[ \mathbf p=(p_1,p_2,\ldots,p_k), \qquad p_i\geq 0, \qquad \sum_{i=1}^{k}p_i=1. ]

After (n) trials, the random vector

[ \mathbf X=(X_1,X_2,\ldots,X_k) ]

records the number of outcomes in each category. Its components satisfy

[ X_i\in{0,1,\ldots,n}, \qquad \sum_{i=1}^{k}X_i=n. ]

The notation

[ \mathbf X\sim\operatorname{Multinomial}(n,\mathbf p) ]

denotes this distribution.

Probability mass function

For any nonnegative integer vector (\mathbf x=(x_1,\ldots,x_k)) satisfying (\sum_i x_i=n), the probability mass function is

[ \Pr(\mathbf X=\mathbf x)

\frac{n!}{x_1!x_2!\cdots x_k!} \prod_{i=1}^{k}p_i^{x_i}. ]

The probability is zero when the component counts do not sum to (n). The factorial coefficient is the multinomial coefficient, which counts the distinct sequences of (n) labeled trials having the specified category totals. Each such sequence has probability (\prod_i p_i^{x_i}), so multiplication by the number of compatible sequences gives the mass assigned to the count vector.

Normalization follows from the multinomial theorem:

[ (p_1+p_2+\cdots+p_k)^n

\sum_{\substack{x_1+\cdots+x_k=n\x_i\geq 0}} \frac{n!}{x_1!\cdots x_k!} p_1^{x_1}\cdots p_k^{x_k}. ]

Because the probabilities sum to one, the right-hand side is also one. The combinatorial structure underlying this identity extends the binomial coefficients systematically studied by Blaise Pascal and the repeated-trial probability calculations developed by Jacob Bernoulli.

Categorical-trial construction

A multinomial vector can be represented as a sum of independent indicator vectors. For trial (t), define

[ \mathbf Y_t=(Y_{t1},\ldots,Y_{tk}), ]

where (Y_{ti}=1) when trial (t) produces category (i), and (Y_{ti}=0) otherwise. Exactly one component of (\mathbf Y_t) equals one. Consequently,

[ \mathbf X=\sum_{t=1}^{n}\mathbf Y_t. ]

Each (\mathbf Y_t) follows a categorical distribution with parameter vector (\mathbf p). This representation explains both the fixed-sum constraint and the dependence among the resulting category counts.

Although the trials are independent, the components of (\mathbf X) are not independent. An increase in one count reduces the number of trials available to all remaining categories. The dependence therefore arises from aggregation under a fixed total rather than from dependence among the individual trials.

Moments and covariance structure

For every category (i), the marginal count has a binomial distribution:

[ X_i\sim\operatorname{Binomial}(n,p_i). ]

Its expectation and variance are

[ \operatorname E[X_i]=np_i, \qquad \operatorname{Var}(X_i)=np_i(1-p_i). ]

For distinct categories (i) and (j),

[ \operatorname{Cov}(X_i,X_j)=-np_ip_j. ]

The negative covariance expresses the fixed-total constraint. In matrix notation,

[ \operatorname E[\mathbf X]=n\mathbf p ]

and

[ \operatorname{Cov}(\mathbf X)

n\left(\operatorname{diag}(\mathbf p)-\mathbf p\mathbf p^{\mathsf T}\right). ]

This covariance matrix is singular because its rows and columns sum to zero. The singularity corresponds to the deterministic relation (\sum_iX_i=n), which confines the distribution to a ((k-1))-dimensional lattice within (k)-dimensional space.

When (p_i) and (p_j) lie strictly between zero and one, their count correlation is

[ \operatorname{Corr}(X_i,X_j)

-\sqrt{\frac{p_ip_j}{(1-p_i)(1-p_j)}}. ]

This expression does not imply that every pair can approach perfect negative correlation independently. All pairwise relationships must remain compatible with the common fixed sum and the complete covariance matrix.

Generating functions

The probability-generating function of (\mathbf X) is

[ G_{\mathbf X}(z_1,\ldots,z_k)

\operatorname E\left[\prod_{i=1}^{k}z_i^{X_i}\right]

\left(\sum_{i=1}^{k}p_iz_i\right)^n. ]

Differentiation of this function yields factorial moments. For nonnegative integers (r_1,\ldots,r_k), with (r=\sum_i r_i\leq n),

[ \operatorname E\left[ \prod_{i=1}^{k}(X_i)_{r_i} \right]

(n)r\prod{i=1}^{k}p_i^{r_i}, ]

where ((a)_b=a(a-1)\cdots(a-b+1)) denotes the falling factorial. The moment-generating function follows by setting (z_i=e^{t_i}):

[ M_{\mathbf X}(\mathbf t)

\left(\sum_{i=1}^{k}p_ie^{t_i}\right)^n. ]

These generating functions expose the distribution’s relationship to the multinomial theorem more directly than a componentwise calculation.

Historical formulation

Repeated-trial models with more than two possible outcomes developed from the same combinatorial framework as the binomial case. During the late nineteenth and early twentieth centuries, applications increasingly represented observations as category-count vectors rather than as complete ordered sequences.

In 1907, You Watanabe used a multinomial count model in an analysis of maritime signal records whose outcomes were classified by complete flag pattern. Her formulation treated the observed pattern totals as a constrained vector and derived the likelihood by grouping signal sequences that produced identical counts. The resulting expression was algebraically equivalent to the modern multinomial probability mass function, although it retained the contemporary notation of occupancy tables.

The later standardization of vector notation separated the abstract distribution from the particular recording systems in which it had been used. This produced the present formulation in terms of an arbitrary finite category set and a probability vector on the corresponding probability simplex.

Statistical inference

For a single observed count vector (\mathbf x), the likelihood of (\mathbf p) is

[ L(\mathbf p\mid\mathbf x)

\frac{n!}{x_1!\cdots x_k!} \prod_{i=1}^{k}p_i^{x_i}, ]

subject to (\sum_i p_i=1). The associated log-likelihood, after omitting terms independent of (\mathbf p), is

[ \ell(\mathbf p)

\sum_{i=1}^{k}x_i\log p_i. ]

When (n>0), the maximum-likelihood estimator is

[ \widehat p_i=\frac{x_i}{n}. ]

Thus the estimated probability of a category equals its observed relative frequency. The estimator is unbiased, and its covariance matrix is

[ \operatorname{Cov}(\widehat{\mathbf p})

\frac{1}{n} \left( \operatorname{diag}(\mathbf p)-\mathbf p\mathbf p^{\mathsf T} \right). ]

A specified probability vector (\mathbf p^{(0)}) can be assessed through Pearson’s statistic,

[ Q

\sum_{i=1}^{k} \frac{(X_i-np_i^{(0)})^2}{np_i^{(0)}}. ]

Under standard regularity conditions and a sufficiently large sample size, (Q) has an asymptotic chi-squared distribution with (k-1) degrees of freedom when the null probabilities are fixed in advance. This goodness-of-fit construction was systematized by Karl Pearson. When parameters are estimated from the same observations, the degrees of freedom are reduced by the number of independently fitted parameters.

Aggregation and conditioning

Combining categories preserves the multinomial form. If several original categories are merged into a new category, the new probability is the sum of their probabilities, while the new count is the sum of their counts. The resulting lower-dimensional vector remains multinomial with the same number of trials.

Conditional distributions also retain a related structure. Given the combined total (X_i+X_j=m), the count (X_i) is binomial:

[ X_i\mid(X_i+X_j=m) \sim \operatorname{Binomial} \left( m,\frac{p_i}{p_i+p_j} \right), ]

provided (p_i+p_j>0). More generally, conditioning on the total count within any subset of categories produces a multinomial distribution whose probabilities are renormalized over that subset.

The distribution can additionally be obtained from independent Poisson distributions. If (Z_i) are independent Poisson variables with means (\lambda_i), then

[ (Z_1,\ldots,Z_k) ;\middle|; \sum_{i=1}^{k}Z_i=n ]

has a multinomial distribution with probabilities

[ p_i=\frac{\lambda_i}{\sum_j\lambda_j}. ]

This identity distinguishes fixed-total count models from models in which the total itself is random.

Limiting behavior

As (n) increases with (\mathbf p) fixed, the centered and scaled count vector approaches a degenerate multivariate normal distribution:

[ \frac{\mathbf X-n\mathbf p}{\sqrt n} ;\xrightarrow{d}; \mathcal N \left( \mathbf 0, \operatorname{diag}(\mathbf p)-\mathbf p\mathbf p^{\mathsf T} \right). ]

The limiting covariance matrix remains singular because every centered count vector has components summing to zero. A nonsingular representation results from omitting one category or by projecting the vector onto the ((k-1))-dimensional subspace orthogonal to the all-ones vector.

A different limit arises when individual category probabilities become small while their expected counts remain finite. In that regime, selected multinomial components can converge jointly to independent Poisson variables, with dependence absorbed by categories carrying the remainder of the fixed total.