Multivariate Normal Distribution
The covariance matrix and the
correlation coefficient give us tools for
quantifying linear relationships among random variables. For a single variable, the
Gaussian \(\mathcal{N}(\mu, \sigma^2)\) is
characterized by its mean and variance. The multivariate normal distribution (MVN) is the natural extension of the
Gaussian to random vectors, and it is arguably the most important joint distribution in all of statistics and machine learning. Its
shape is entirely determined by a mean vector and a covariance matrix. Everything about its geometry therefore rests on the covariance
and correlation machinery already in hand.
Definition: Multivariate Normal Distribution
A \(D\)-dimensional random vector \(\boldsymbol{x} \in \mathbb{R}^D\) follows a multivariate normal distribution,
written \(\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)\), if its p.d.f. is
\[
f(\boldsymbol{x}) = \frac{1}{\sqrt{(2\pi)^D \det(\Sigma)}}
\exp\!\left[
-\frac{1}{2}(\boldsymbol{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{\mu})
\right],
\]
where \(\boldsymbol{\mu} = \mathbb{E}[\boldsymbol{x}] \in \mathbb{R}^D\) is the mean vector and
\(\Sigma = \operatorname{Cov}[\boldsymbol{x}] \in \mathbb{R}^{D \times D}\) is the covariance matrix, here required to be
strictly positive definite so that \(\Sigma^{-1}\) exists and \(\det(\Sigma) \gt 0\).
When \(\Sigma\) is only positive semi-definite, the distribution becomes degenerate and is supported on a lower-dimensional
affine subspace, which calls for a different formulation that we do not pursue here. Positive semi-definiteness itself is guaranteed
in general by the
PSD theorem for covariance matrices.
The expression inside the exponential (ignoring the factor of \(-\tfrac{1}{2}\)) is the squared
Mahalanobis distance between \(\boldsymbol{x}\) and \(\boldsymbol{\mu}\).
Definition: Mahalanobis Distance
Let \(\Sigma \in \mathbb{R}^{D \times D}\) be a strictly positive definite matrix. The
Mahalanobis distance between \(\boldsymbol{x}, \boldsymbol{y} \in \mathbb{R}^D\) (with respect to \(\Sigma\)) is
\[
d_{\Sigma}(\boldsymbol{x}, \boldsymbol{y})
= \sqrt{(\boldsymbol{x} - \boldsymbol{y})^\top \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{y})}.
\]
Unlike the Euclidean distance, the Mahalanobis distance accounts for the correlation
structure among variables. The contours of constant density of \(\mathcal{N}(\boldsymbol{\mu}, \Sigma)\) are
ellipsoids defined by
\[
(\boldsymbol{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) = c,
\]
whose axes are aligned with the eigenvectors of \(\Sigma\) and whose lengths are
proportional to the square roots of the corresponding eigenvalues.
To build concrete intuition, let us work out the special case \(D = 2\) in full detail. When \(\boldsymbol{x} \in \mathbb{R}^2\), the
MVN is called the bivariate normal distribution. In this case,
\[
\begin{align*}
\Sigma &= \begin{bmatrix}
\operatorname{Var}[X_1] & \operatorname{Cov}[X_1, X_2] \\\\
\operatorname{Cov}[X_2, X_1] & \operatorname{Var}[X_2]
\end{bmatrix} \\\\
&= \begin{bmatrix}
\sigma_1^2 & \rho \sigma_1 \sigma_2 \\\\
\rho \sigma_1 \sigma_2 & \sigma_2^2
\end{bmatrix}
\end{align*}
\]
where \(\rho\) is the correlation coefficient
of \(X_1\) and \(X_2\),
\[
\begin{align*}
\rho
&= \operatorname{Corr}[X_1, X_2] \\\\
&= \frac{\operatorname{Cov}[X_1, X_2]}{\sqrt{\operatorname{Var}[X_1]\,\operatorname{Var}[X_2]}}.
\end{align*}
\]
Then
\[
\begin{align*}
\det (\Sigma) &= \sigma_1^2 \sigma_2^2 - \rho^2 \sigma_1^2 \sigma_2^2 \\\\
&= \sigma_1^2 \sigma_2^2 (1 - \rho^2)
\end{align*}
\]
and
\[
\begin{align*}
\Sigma^{-1} &= \frac{1}{\det (\Sigma )}
\begin{bmatrix}
\sigma_2^2 & -\rho \sigma_1 \sigma_2 \\\\
-\rho \sigma_1 \sigma_2 & \sigma_1^2
\end{bmatrix} \\\\
&= \frac{1}{1 - \rho^2}
\begin{bmatrix}
\frac{1}{\sigma_1^2 } & \frac{-\rho} {\sigma_1 \sigma_2} \\\\
\frac{-\rho} {\sigma_1 \sigma_2} & \frac{1}{\sigma_2^2 }
\end{bmatrix}
\end{align*}
\]
Note that the exponent in the bivariate density is a quadratic form:
\[
(\boldsymbol{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\boldsymbol{x} -\boldsymbol{\mu}).
\]
Expanding,
\[
\begin{align*}
(\boldsymbol{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\boldsymbol{x} -\boldsymbol{\mu})
&= \frac{1}{1 - \rho^2} \begin{bmatrix} X_1 - \mu_1 & X_2 - \mu_2 \end{bmatrix}
\begin{bmatrix}
\frac{1}{\sigma_1^2 } & \frac{-\rho} {\sigma_1 \sigma_2} \\\\
\frac{-\rho} {\sigma_1 \sigma_2} & \frac{1}{\sigma_2^2 }
\end{bmatrix}
\begin{bmatrix} X_1 - \mu_1 \\ X_2 - \mu_2 \end{bmatrix} \\\\
&= \frac{1}{1 - \rho^2}\left[\frac{1}{\sigma_1^2 }(X_1 - \mu_1)^2
-\frac{2\rho} {\sigma_1 \sigma_2}(X_1 - \mu_1)(X_2 - \mu_2)
+\frac{1}{\sigma_2^2 }(X_2 - \mu_2)^2 \right].
\end{align*}
\]
Therefore, the p.d.f. of the bivariate normal distribution becomes
\[
f(\boldsymbol{x}) = \frac{1}{2\pi \sigma_1 \sigma_2 \sqrt{1 - \rho^2}}
\exp\!\left\{-\frac{1}{2(1 - \rho^2)}
\left[\left(\frac{X_1 - \mu_1}{\sigma_1}\right)^2
-2\rho \left(\frac{X_1 - \mu_1} {\sigma_1}\right) \left(\frac{X_2 - \mu_2} {\sigma_2}\right)
+\left(\frac{X_2 - \mu_2}{\sigma_2}\right)^2
\right]
\right\}.
\]
When \(\rho = -1\) or \(\rho = 1\), this density is undefined and the distribution is degenerate.
Cholesky Decomposition
The previous section defined the MVN in terms of its density function, but in practice we frequently need to sample from
\[
\mathcal{N}(\boldsymbol{\mu}, \Sigma).
\]
Sampling is needed when training variational autoencoders, when implementing the reparameterization trick, and when running
Monte Carlo simulations.
Generating independent standard normal samples is straightforward, but we need a way to introduce the correlation structure encoded in
\(\Sigma\). The Cholesky decomposition provides a numerically stable factorization that transforms uncorrelated samples
into correlated ones.
Theorem: Cholesky Decomposition
Let \(\boldsymbol{A} \in \mathbb{R}^{D \times D}\) be a symmetric,
positive definite
matrix. Then there exists a unique lower triangular matrix \(\boldsymbol{L}\) with positive diagonal entries such that
\[
\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^\top.
\]
(The factorization can equivalently be written as \(\boldsymbol{A} = \boldsymbol{R}^\top \boldsymbol{R}\) where
\(\boldsymbol{R} = \boldsymbol{L}^\top\) is upper triangular.)
Conversely, if \(\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^\top\) for some lower triangular matrix \(\boldsymbol{L}\)
with positive diagonal entries, then \(\boldsymbol{A}\) is symmetric and positive definite.
Proof:
Existence. We argue by induction on \(D\). For \(D = 1\), write \(\boldsymbol{A} = [a]\). Positive definiteness
gives \(a = [1]^\top \boldsymbol{A}\,[1] \gt 0\), so \(\boldsymbol{L} = [\sqrt{a}]\) is the required factor.
Now let \(D \geq 2\) and assume the result for symmetric positive definite matrices of size \(D - 1\). Partition
\[
\boldsymbol{A} = \begin{bmatrix} a & \boldsymbol{b}^\top \\ \boldsymbol{b} & \boldsymbol{C} \end{bmatrix},
\]
where \(a \in \mathbb{R}\), \(\boldsymbol{b} \in \mathbb{R}^{D-1}\), and \(\boldsymbol{C} \in \mathbb{R}^{(D-1) \times (D-1)}\)
is symmetric because \(\boldsymbol{A}\) is. Evaluating the quadratic form of \(\boldsymbol{A}\) at the first standard basis vector
\(\boldsymbol{e}_1\) gives \(a = \boldsymbol{e}_1^\top \boldsymbol{A} \boldsymbol{e}_1 \gt 0\). We may therefore form the
Schur complement
of the entry \(a\) in \(\boldsymbol{A}\),
\[
\boldsymbol{S} = \boldsymbol{C} - \frac{1}{a}\boldsymbol{b}\boldsymbol{b}^\top,
\]
which is symmetric because \(\boldsymbol{C}\) and \(\boldsymbol{b}\boldsymbol{b}^\top\) are.
The matrix \(\boldsymbol{S}\) is also positive definite. Given a nonzero \(\boldsymbol{w} \in \mathbb{R}^{D-1}\), set
\(t = -\boldsymbol{b}^\top \boldsymbol{w} / a\) and let \(\boldsymbol{v} \in \mathbb{R}^D\) be the vector whose first entry is
\(t\) and whose remaining entries are those of \(\boldsymbol{w}\). Then \(\boldsymbol{v} \neq \boldsymbol{0}\) because \(\boldsymbol{w} \neq \boldsymbol{0}\), and
block multiplication
gives
\[
\begin{align*}
\boldsymbol{v}^\top \boldsymbol{A} \boldsymbol{v}
&= a t^2 + 2t\,\boldsymbol{b}^\top \boldsymbol{w} + \boldsymbol{w}^\top \boldsymbol{C} \boldsymbol{w} \\\\
&= \boldsymbol{w}^\top \boldsymbol{C} \boldsymbol{w} - \frac{1}{a}\left(\boldsymbol{b}^\top \boldsymbol{w}\right)^2 \\\\
&= \boldsymbol{w}^\top \boldsymbol{S} \boldsymbol{w}.
\end{align*}
\]
The first equality also uses \(\boldsymbol{w}^\top \boldsymbol{b} = \boldsymbol{b}^\top \boldsymbol{w}\), the second substitutes
the value of \(t\), and the last uses
\(\left(\boldsymbol{b}^\top \boldsymbol{w}\right)^2 = \boldsymbol{w}^\top \boldsymbol{b}\boldsymbol{b}^\top \boldsymbol{w}\)
together with the definition of \(\boldsymbol{S}\).
Since \(\boldsymbol{A}\) is positive definite, \(\boldsymbol{w}^\top \boldsymbol{S} \boldsymbol{w} \gt 0\).
By the induction hypothesis, \(\boldsymbol{S} = \boldsymbol{L}'\boldsymbol{L}'^\top\) for a lower triangular matrix
\(\boldsymbol{L}'\) with positive diagonal entries. Define
\[
\boldsymbol{L} = \begin{bmatrix} \sqrt{a} & \boldsymbol{0}^\top \\ \boldsymbol{b}/\sqrt{a} & \boldsymbol{L}' \end{bmatrix},
\]
which is lower triangular with positive diagonal entries. Transposing block by block, we have
\[
\boldsymbol{L}^\top = \begin{bmatrix} \sqrt{a} & \boldsymbol{b}^\top/\sqrt{a} \\ \boldsymbol{0} & \boldsymbol{L}'^\top \end{bmatrix},
\]
and block multiplication, followed by \(\boldsymbol{L}'\boldsymbol{L}'^\top = \boldsymbol{S}\), gives
\[
\begin{align*}
\boldsymbol{L}\boldsymbol{L}^\top
&= \begin{bmatrix} a & \boldsymbol{b}^\top \\ \boldsymbol{b} & \frac{1}{a}\boldsymbol{b}\boldsymbol{b}^\top + \boldsymbol{L}'\boldsymbol{L}'^\top \end{bmatrix} \\\\
&= \begin{bmatrix} a & \boldsymbol{b}^\top \\ \boldsymbol{b} & \boldsymbol{C} \end{bmatrix} \\\\
&= \boldsymbol{A}.
\end{align*}
\]
Uniqueness. We again argue by induction on \(D\). For \(D = 1\), the single entry \(\ell\) of
\(\boldsymbol{L}\) must satisfy \(\ell^2 = a\) and \(\ell \gt 0\), so \(\ell = \sqrt{a}\). For \(D \geq 2\), keep the
partition of \(\boldsymbol{A}\) above and let \(\boldsymbol{L}\) be any lower triangular matrix with positive diagonal
entries such that \(\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^\top\). Partition it conformally as
\[
\boldsymbol{L} = \begin{bmatrix} \ell & \boldsymbol{0}^\top \\ \boldsymbol{m} & \boldsymbol{M} \end{bmatrix},
\]
where \(\ell \gt 0\), \(\boldsymbol{m} \in \mathbb{R}^{D-1}\), and \(\boldsymbol{M}\) is lower triangular with positive
diagonal entries. Taking the transpose blockwise and multiplying out, we find
\[
\boldsymbol{L}\boldsymbol{L}^\top
= \begin{bmatrix} \ell^2 & \ell\,\boldsymbol{m}^\top \\ \ell\,\boldsymbol{m} & \boldsymbol{m}\boldsymbol{m}^\top + \boldsymbol{M}\boldsymbol{M}^\top \end{bmatrix}.
\]
Comparing blocks with the partition of \(\boldsymbol{A}\), we see that \(\ell^2 = a\), which forces \(\ell = \sqrt{a}\)
because \(\ell \gt 0\). The equation \(\ell\,\boldsymbol{m} = \boldsymbol{b}\) then gives \(\boldsymbol{m} = \boldsymbol{b}/\sqrt{a}\),
and so \(\boldsymbol{M}\boldsymbol{M}^\top = \boldsymbol{C} - \boldsymbol{m}\boldsymbol{m}^\top = \boldsymbol{S}\). The existence
part showed that \(\boldsymbol{S}\) is symmetric positive definite, so the induction hypothesis determines \(\boldsymbol{M}\)
uniquely. Hence \(\boldsymbol{L}\) is unique.
Converse. Suppose \(\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^\top\) with \(\boldsymbol{L}\) as in the
statement. Then
\(\boldsymbol{A}^\top = (\boldsymbol{L}^\top)^\top \boldsymbol{L}^\top = \boldsymbol{L}\boldsymbol{L}^\top = \boldsymbol{A}\),
and for every \(\boldsymbol{u} \in \mathbb{R}^D\),
\[
\boldsymbol{u}^\top \boldsymbol{A} \boldsymbol{u} = \left(\boldsymbol{L}^\top \boldsymbol{u}\right)^\top \left(\boldsymbol{L}^\top \boldsymbol{u}\right),
\]
the sum of the squares of the entries of \(\boldsymbol{L}^\top \boldsymbol{u}\). It remains to show that
\(\boldsymbol{L}^\top \boldsymbol{u} \neq \boldsymbol{0}\) whenever \(\boldsymbol{u} \neq \boldsymbol{0}\). The \(i\)-th entry
of \(\boldsymbol{L}^\top \boldsymbol{u}\) is \(L_{ii} u_i + \sum_{j \gt i} L_{ji} u_j\), because \(L_{ji} = 0\) for
\(j \lt i\). Suppose \(\boldsymbol{L}^\top \boldsymbol{u} = \boldsymbol{0}\). The last entry reads \(L_{DD} u_D = 0\), and
\(L_{DD} \gt 0\) gives \(u_D = 0\). Working upward through \(i = D-1, \ldots, 1\), the entries \(u_j\) with \(j \gt i\) are already
known to vanish, so we obtain \(L_{ii} u_i = 0\) at each step, and \(u_i = 0\) because \(L_{ii} \gt 0\). Thus \(\boldsymbol{u} = \boldsymbol{0}\), and
\(\boldsymbol{u}^\top \boldsymbol{A} \boldsymbol{u} \gt 0\) for every nonzero \(\boldsymbol{u}\).
Returning to sampling, the Cholesky factor \(\boldsymbol{L}\) of \(\Sigma\) provides a recipe for converting standard normal
samples into MVN samples with arbitrary mean and covariance.
Theorem: MVN Sampling via Cholesky Reparametrization
Let \(\Sigma \in \mathbb{R}^{D \times D}\) be strictly positive definite with Cholesky decomposition
\(\Sigma = \boldsymbol{L}\boldsymbol{L}^\top\), and let \(\boldsymbol{\mu} \in \mathbb{R}^D\).
If \(\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\) (obtained by sampling
\(D\) independent standard normal variables), then
\[
\boldsymbol{y} = \boldsymbol{L}\boldsymbol{x} + \boldsymbol{\mu}
\]
is a multivariate normal random vector with mean \(\boldsymbol{\mu}\) and covariance \(\Sigma\):
\[
\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma).
\]
Proof:
By linearity of expectation, applied
component-wise and iterated over the \(D\) terms of each linear combination,
\[
\begin{align*}
\mathbb{E}[\boldsymbol{y}]
&= \mathbb{E}[\boldsymbol{L}\boldsymbol{x} + \boldsymbol{\mu}] \\\\
&= \boldsymbol{L}\,\mathbb{E}[\boldsymbol{x}] + \boldsymbol{\mu} \\\\
&= \boldsymbol{\mu},
\end{align*}
\]
since \(\mathbb{E}[\boldsymbol{x}] = \boldsymbol{0}\).
For the covariance, by the
definition of the covariance matrix,
\[
\operatorname{Cov}[\boldsymbol{y}]
= \mathbb{E}\!\left[(\boldsymbol{y} - \mathbb{E}[\boldsymbol{y}])
(\boldsymbol{y} - \mathbb{E}[\boldsymbol{y}])^\top \right].
\]
Substituting \(\boldsymbol{y} = \boldsymbol{L}\boldsymbol{x} + \boldsymbol{\mu}\) and using
\(\mathbb{E}[\boldsymbol{x}] = \boldsymbol{0}\) (so that
\(\operatorname{Cov}[\boldsymbol{x}] = \mathbb{E}[\boldsymbol{x}\boldsymbol{x}^\top]\)), the matrix-valued linearity of expectation
(entry-wise application of scalar linearity) gives
\[
\begin{align*}
\operatorname{Cov}[\boldsymbol{y}]
&= \mathbb{E}\!\left[ (\boldsymbol{L}\boldsymbol{x})
(\boldsymbol{L}\boldsymbol{x})^\top \right] \\\\
&= \boldsymbol{L}\,\mathbb{E}\!\left[\boldsymbol{x}\boldsymbol{x}^\top\right]\boldsymbol{L}^\top \\\\
&= \boldsymbol{L}\,\operatorname{Cov}[\boldsymbol{x}]\,\boldsymbol{L}^\top \\\\
&= \boldsymbol{L}\,\boldsymbol{I}\,\boldsymbol{L}^\top \\\\
&= \boldsymbol{L}\boldsymbol{L}^\top \\\\
&= \Sigma.
\end{align*}
\]
Finally, \(\boldsymbol{y}\) is multivariate normal because it is an invertible affine transformation of a standard MVN
(\(\boldsymbol{L}\) has positive diagonal entries). Each component is a linear combination of independent normals plus a
constant, and invertible affine transformations preserve normality (a property of the Gaussian family that we take for granted
here). Therefore \(\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)\).
Dirichlet Distribution
The multivariate normal distribution models random vectors taking values in all of \(\mathbb{R}^D\).
However, many quantities of interest in machine learning are probability vectors, whose
entries are non-negative and sum to one. For example, the class probabilities output by a softmax layer,
the topic proportions in a document, or the mixing weights in a mixture model all live on a
probability simplex. We need a distribution that respects this constraint.
The Dirichlet distribution is a multivariate generalization of the
beta distribution. It has support over the
\((K - 1)\)-dimensional probability simplex, defined by
\[
S_K = \Bigl\{(x_1, x_2, \ldots, x_K) \in \mathbb{R}^K: x_k \ge 0,\ \sum_{k=1}^K x_k = 1 \Bigr\}.
\]
Definition: Dirichlet Distribution
Let \(K \geq 2\). A random vector \(\boldsymbol{x} \in \mathbb{R}^K\) taking values in \(S_K\) has a Dirichlet
distribution with parameters \(\boldsymbol{\alpha} = (\alpha_1, \alpha_2, \ldots, \alpha_K)\) (with each \(\alpha_k \gt 0\)),
written \(\boldsymbol{x} \sim \operatorname{Dir}(\boldsymbol{\alpha})\), if, at the points of \(S_K\) where every
\(x_k \gt 0\), its probability density function is given by
\[
f(x_1, \ldots, x_K; \boldsymbol{\alpha}) = \frac{1}{B(\boldsymbol{\alpha})} \prod_{k=1}^K x_k^{\alpha_k - 1}, \quad (x_1, \ldots, x_K) \in S_K,
\]
or equivalently
\[
\operatorname{Dir}(\boldsymbol{x} \mid \boldsymbol{\alpha})
= \frac{1}{B(\boldsymbol{\alpha})} \prod_{k=1}^K x_k^{\alpha_k - 1}\, \mathbb{1}\{\boldsymbol{x} \in S_K,\ \min_k x_k \gt 0\},
\]
where the multivariate beta function \(B(\boldsymbol{\alpha})\) is defined as
\[
B(\boldsymbol{\alpha}) = \frac{\prod_{k=1}^K \Gamma(\alpha_k)}{\Gamma\Bigl(\sum_{k=1}^K \alpha_k\Bigr)}.
\]
On \(S_K\) the last coordinate is determined by the others, \(x_K = 1 - \sum_{k=1}^{K-1} x_k\), so the density is a
function of \((x_1, \ldots, x_{K-1})\) and is integrated against Lebesgue measure on the region
\(\{(x_1, \ldots, x_{K-1}) : x_1, \ldots, x_{K-1} \geq 0,\ \sum_{k=1}^{K-1} x_k \leq 1\}\). That \(B(\boldsymbol{\alpha})\) is the correct normalizing constant for this
integral is a fact we take for granted here. For \(K = 2\) it is the
Beta-Gamma identity.
Theorem: Moments of the Dirichlet Distribution
Let \(\boldsymbol{x} \sim \operatorname{Dir}(\boldsymbol{\alpha})\) and let \(\alpha_0 = \sum_{k=1}^K \alpha_k\). Then
for each \(k\),
- Mean:
\[
\mathbb{E}[x_k] = \frac{\alpha_k}{\alpha_0}.
\]
- Variance:
\[
\operatorname{Var}[x_k] = \frac{\alpha_k (\alpha_0 - \alpha_k)}{\alpha_0^2 (\alpha_0+1)}.
\]
For the symmetric Dirichlet prior with \(\alpha_k = \frac{\alpha}{K}\), these reduce to
\[
\mathbb{E}[x_k] = \frac{1}{K}, \quad \operatorname{Var}[x_k] = \frac{K-1}{K^2 (\alpha +1)},
\]
so increasing \(\alpha\) increases the precision (decreases the variance) of the distribution.
- Covariance: for \(i \neq j\),
\[
\operatorname{Cov}[x_i, x_j] = \frac{-\alpha_i \alpha_j}{\alpha_0^2 (\alpha_0+1)}.
\]
Proof:
The argument rests on one fact about the Dirichlet family that we do not derive here. If \(A \subseteq \{1, \ldots, K\}\) is
non-empty and proper, then
\[
\sum_{k \in A} x_k \sim \operatorname{Beta}\!\left(\sum_{k \in A} \alpha_k,\ \alpha_0 - \sum_{k \in A} \alpha_k\right).
\]
Taking \(A = \{k\}\) leaves \(x_k \sim \operatorname{Beta}(\alpha_k, \alpha_0 - \alpha_k)\), whose two parameters sum to
\(\alpha_0\). The mean and variance recorded in the
beta distribution therefore give
\[
\begin{align*}
\mathbb{E}[x_k] &= \frac{\alpha_k}{\alpha_0}, \\\\
\operatorname{Var}[x_k] &= \frac{\alpha_k (\alpha_0 - \alpha_k)}{\alpha_0^2 (\alpha_0 + 1)}.
\end{align*}
\]
Setting \(\alpha_k = \alpha/K\) makes \(\alpha_0 = \alpha\), and the two expressions collapse to
\(1/K\) and \((K-1)/(K^2 (\alpha + 1))\).
For \(i \neq j\) the covariance follows the same route. Taking \(A = \{i, j\}\) gives
\[
\operatorname{Var}[x_i + x_j]
= \frac{(\alpha_i + \alpha_j)(\alpha_0 - \alpha_i - \alpha_j)}{\alpha_0^2 (\alpha_0 + 1)},
\]
while the variance of a sum expands the left-hand
side into \(\operatorname{Var}[x_i] + \operatorname{Var}[x_j] + 2\operatorname{Cov}[x_i, x_j]\). The finite-second-moment hypothesis
holds here because every coordinate lies in \([0, 1]\). Subtracting the two individual variances and halving yields
\[
\operatorname{Cov}[x_i, x_j] = \frac{-\alpha_i \alpha_j}{\alpha_0^2 (\alpha_0 + 1)}.
\]
The parameters \(\alpha_k\) can be thought of as "pseudocounts" or prior observations of each category. When all \(\alpha_k\) are equal
(that is, \(\boldsymbol{\alpha} = \alpha\, \boldsymbol{1}\)), the distribution is said to be symmetric, and reduces to
the uniform distribution over the simplex when \(\alpha = 1\). This symmetry makes the Dirichlet distribution a natural
prior in Bayesian models where no category is favored a priori.
The Dirichlet distribution is widely used as a conjugate prior for the parameters
of a multinomial distribution in Bayesian statistics.
Wishart Distribution
The Dirichlet distribution placed a prior on probability vectors. In Bayesian multivariate analysis, we also need priors on
covariance matrices themselves. Inferring the parameters of a multivariate normal model is one such case. The
Wishart distribution fills this role. It generalizes the
gamma distribution from positive scalars to
positive definite matrices.
In this section we follow the standard machine-learning notation. Here \(D\) denotes the matrix dimension and \(N\) the sample size.
Definition: Multivariate Gamma Function
For an integer \(D \geq 1\) and real \(a \gt (D-1)/2\), the multivariate gamma function is
\[
\Gamma_D(a) = \pi^{D(D-1)/4} \prod_{j=1}^{D} \Gamma\!\left(a + \frac{1-j}{2}\right),
\]
where \(\Gamma\) is the gamma function. The
condition \(a \gt (D-1)/2\) makes every argument \(a + (1-j)/2\) positive. For \(D = 1\) the product has the single
factor \(\Gamma(a)\), so \(\Gamma_1 = \Gamma\). The function enters the normalizing constant of the Wishart density below,
and that this constant makes the density integrate to one is a fact we take for granted here.
Definition: Wishart Distribution
A \(D \times D\) symmetric positive definite random matrix \(\Sigma\) follows a
Wishart distribution with parameters \(\boldsymbol{S}\) (scale matrix) and
\(\nu\) (degrees of freedom) if its p.d.f. is
\[
\operatorname{Wi}(\Sigma \mid \boldsymbol{S}, \nu)
= \frac{1}{Z}\, |\Sigma|^{(\nu - D - 1)/2} \exp\!\left(-\frac{1}{2}\operatorname{tr}(\boldsymbol{S}^{-1} \Sigma)\right),
\]
where the normalization constant is
\[
Z = |\boldsymbol{S}|^{\nu/2}\, 2^{\nu D/2}\, \Gamma_D\!\left(\frac{\nu}{2}\right),
\]
and \(\Gamma_D\) is the multivariate gamma function. The parameters are constrained as follows:
- \(\nu\): the degrees of freedom, which must satisfy \(\nu \gt D - 1\) for the normalization to exist.
- \(\boldsymbol{S}\): the scale matrix, a \(D \times D\) symmetric positive definite matrix.
With these parameters, the distribution has mean \(\nu \boldsymbol{S}\) and (when
\(\nu \gt D + 1\)) mode \((\nu - D - 1)\boldsymbol{S}\).
If \(D = 1\), the Wishart reduces to the gamma distribution:
\[
\operatorname{Wi}(\lambda \mid s, \nu) = \operatorname{Gamma}\!\left(\lambda \,\Big|\, \text{shape} = \frac{\nu}{2},\ \text{rate} = \frac{1}{2s}\right),
\]
so that the mean \(\nu s\) agrees with the general formula \(\nu \boldsymbol{S}\) specialized to \(\boldsymbol{S} = s\).
Setting \(s = 1\) further reduces this to the chi-squared distribution \(\chi^2_\nu\),
whose density is \(\operatorname{Gamma}(\text{shape} = \nu/2,\ \text{rate} = 1/2)\).
Applications
In the applications that follow, \(\Sigma\) denotes the covariance matrix of a multivariate normal model rather than the
Wishart-distributed matrix of the definition above.
- Covariance Matrix Estimation:
In multivariate statistics, the Wishart distribution arises naturally in the sampling distribution
of empirical
covariance matrices.
Specifically, given \(N \geq D\) i.i.d. samples
\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N \sim \mathcal{N}(\boldsymbol{0}, \Sigma)\), the scatter matrix
\[
\boldsymbol{S}_0 = \sum_{n=1}^N \boldsymbol{x}_n \boldsymbol{x}_n^\top
\]
(the subscript \(0\) emphasizes that this is the \(\boldsymbol{\mu} = \boldsymbol{0}\) centered scatter, distinct from the
Wishart scale parameter \(\boldsymbol{S}\) introduced above) follows a Wishart distribution with scale \(\Sigma\) and \(N\)
degrees of freedom:
\[
\boldsymbol{S}_0 \sim \operatorname{Wi}(\Sigma, N).
\]
- Bayesian Inference:
The Wishart serves as a conjugate prior on the precision matrix \(\Sigma^{-1}\) in multivariate normal models with
known mean, enabling closed-form updates of the posterior distribution.
For a prior on the covariance matrix \(\Sigma\) directly, we use the inverse Wishart distribution: for
\(\nu \gt D - 1\) and symmetric positive definite \(\boldsymbol{S}\),
\[
\operatorname{IW}(\Sigma \mid \boldsymbol{S}, \nu)
= \frac{1}{Z_{IW}}\, |\Sigma|^{-(\nu + D + 1)/2} \exp\!\left(-\frac{1}{2}\operatorname{tr}(\boldsymbol{S}\, \Sigma^{-1})\right),
\]
where
\[
Z_{IW} = |\boldsymbol{S}|^{-\nu/2}\, 2^{\nu D/2}\, \Gamma_D\!\left(\frac{\nu}{2}\right).
\]
The mean of this distribution is \(\boldsymbol{S} / (\nu - D - 1)\) for \(\nu \gt D + 1\). Its relationship to the Wishart is
the matrix analogue of the scalar fact
\(\lambda \sim \operatorname{Gamma}(a, b) \Rightarrow 1/\lambda \sim \operatorname{IG}(a, b)\). In particular,
\(\Sigma \sim \operatorname{IW}(\boldsymbol{S}, \nu)\) if and only if
\(\Sigma^{-1} \sim \operatorname{Wi}(\boldsymbol{S}^{-1}, \nu)\), so the inverse Wishart is the multivariate generalization of
the inverse gamma. Specializing to \(D = 1\) confirms this:
\[
\operatorname{IW}(\sigma^2 \mid s, \nu) = \operatorname{IG}\!\left(\sigma^2 \,\Big|\, \text{shape} = \frac{\nu}{2},\ \text{scale} = \frac{s}{2}\right),
\]
with mean \(s/(\nu - 2)\) for \(\nu \gt 2\), agreeing with the general formula \(\boldsymbol{S}/(\nu - D - 1)\) at
\(\boldsymbol{S} = s\), \(D = 1\). This mirrors the Wishart-to-gamma reduction shown earlier.
Note that the first application concerns the sampling distribution of an estimator (a frequentist usage), while the
second concerns prior beliefs and posterior inference (a Bayesian usage). Both deal with uncertainty in covariance
matrices, and the same distribution arises in both roles, with different interpretations.
Connections to Machine Learning
The multivariate distributions introduced here form the backbone of probabilistic machine learning.
Gaussian processes use the MVN to define a prior over
functions, and it reappears in the E-step of Gaussian mixture models and in the reparameterization
trick for variational autoencoders.
The Dirichlet distribution is the standard prior for topic models such as latent Dirichlet allocation (LDA) and for categorical
distributions in Bayesian inference. The Wishart and inverse Wishart distributions are essential whenever the covariance matrix
itself is a parameter to be inferred, as in Bayesian multivariate regression and Bayesian Gaussian mixture models.
Covariance, correlation, and the principal multivariate distributions are now in place. With that framework we have the
probabilistic foundations needed for statistical inference. These foundations are put to work in
maximum likelihood estimation, the principal method for fitting probabilistic models
to observed data.