Likelihood of a multivariate normal sample

Section of Estimation and inference for the multivariate normal, part of Statistical learning

Type	Section
Taught in	Estimation and inference for the multivariate normal
Part	Statistical learning

A likelihood holds the observed sample fixed and varies the parameters. Parameter values that give this sample greater density have greater likelihood. That comparison does not assign a probability distribution to the parameters. For continuous data the likelihood is built from density values, which may exceed one. An estimator is the rule applied to a random sample. An estimate is the number that rule returns for the sample in hand. Sampling uncertainty belongs to the estimator.

Model the sample as n independent draws \mathbf{x}_1,\ldots,\mathbf{x}_n from \operatorname{N}_p(\boldsymbol{\mu},\Sigma), and stack the realised vectors as the rows of a data matrix,

\mathbf{X}_{\mathrm{all}} =\begin{pmatrix} x_{11}&\cdots&x_{1p}\\ \vdots&&\vdots\\ x_{n1}&\cdots&x_{np} \end{pmatrix} =\begin{pmatrix}\mathbf{x}_1^\top\\ \vdots\\ \mathbf{x}_n^\top\end{pmatrix}.

Independence multiplies the n normal densities. Hold the rows fixed and regard the product as a function of (\boldsymbol{\mu},\Sigma). That function is the likelihood

L(\boldsymbol{\mu},\Sigma;\mathbf{X}_{\mathrm{all}}) =\frac{1}{(2\pi)^{np/2}|\Sigma|^{n/2}} \exp\Bigl\{-\frac12\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu})\Bigr\}.

The logarithm is strictly increasing, so L and \ell=\log L have the same maximiser, and the derivative is easier on the logarithm:

\ell(\boldsymbol{\mu},\Sigma;\mathbf{X}_{\mathrm{all}}) =-\frac12 np\log(2\pi)-\frac12 n\log|\Sigma| -\frac12\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}).

Products become sums. Terms that do not depend on a parameter can be dropped while maximising in that parameter. They stay when an absolute density value, or a comparison that uses the whole likelihood, is required.

Two summaries of the sample carry the maximiser. The sample mean vector is

\bar{\mathbf{x}}=\frac1n\sum_{i=1}^n\mathbf{x}_i,

and the matrix of corrected sums of squares and cross-products is

A=\sum_{i=1}^n(\mathbf{x}_i-\bar{\mathbf{x}})(\mathbf{x}_i-\bar{\mathbf{x}})^\top =\Bigl(\sum_{i=1}^n\mathbf{x}_i\mathbf{x}_i^\top\Bigr)-n\bar{\mathbf{x}}\bar{\mathbf{x}}^\top.

The second expression is the first one expanded. The cross terms cancel because \sum_i\mathbf{x}_i=n\bar{\mathbf{x}}. In particular A does not depend on (\boldsymbol{\mu},\Sigma).

The mean

Main article: Maximum likelihood for a normal mean vector

The log-likelihood depends on \boldsymbol{\mu} only through the sum of quadratic forms, and \Sigma is an arbitrary positive definite matrix. Split each residual into a deviation from the sample mean and a deviation of the sample mean from \boldsymbol{\mu}:

\begin{aligned} \sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}) &=\sum_{i=1}^n\bigl[(\mathbf{x}_i-\bar{\mathbf{x}})+(\bar{\mathbf{x}}-\boldsymbol{\mu})\bigr]^\top \Sigma^{-1} \bigl[(\mathbf{x}_i-\bar{\mathbf{x}})+(\bar{\mathbf{x}}-\boldsymbol{\mu})\bigr]\\ &=\sum_{i=1}^n(\mathbf{x}_i-\bar{\mathbf{x}})^\top\Sigma^{-1}(\mathbf{x}_i-\bar{\mathbf{x}}) +n(\bar{\mathbf{x}}-\boldsymbol{\mu})^\top\Sigma^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}). \end{aligned}

The cross terms contain \sum_i(\mathbf{x}_i-\bar{\mathbf{x}})=\mathbf{0}, so they vanish. The first sum does not involve \boldsymbol{\mu}. The second is n times a quadratic form in the positive definite matrix \Sigma^{-1}, so it is nonnegative and equals zero only at \boldsymbol{\mu}=\bar{\mathbf{x}}. The maximiser in \boldsymbol{\mu} is therefore

\widehat{\boldsymbol{\mu}}=\bar{\mathbf{x}},

for every positive definite \Sigma. The covariance does not have to be known before the mean is estimated. The next step substitutes this maximiser and maximises in \Sigma alone.

The covariance and the invariance lemma

Main article: Maximum likelihood for a covariance matrix

The inverse covariance is a convenient optimisation variable. Set \Theta=\Sigma^{-1}. The quadratic part of the log-likelihood is linear in \Theta, and the determinant term becomes \log|\Theta|. Together they form a concave function on the positive definite matrices. The maximiser can be found and then inverted. The invariance result says why the inverted matrix is the maximum likelihood estimator of \Sigma. The general statement is the invariance of maximum likelihood. The argument used for the covariance is the following.

Proposition 1 (Invariance of the maximum likelihood estimator) Let the distribution of the data depend on \theta=(\theta_1,\ldots,\theta_k), with maximum likelihood estimator \widehat\theta. Let \tau(\theta)=(\tau_1(\theta),\ldots,\tau_r(\theta)), 1\le r\le k, be a transformation of the parameter. Then the maximum likelihood estimator of \tau_l(\theta) is \widehat\tau_l=\tau_l(\widehat\theta).

Article: Invariance of the maximum likelihood estimator →

Let M(\tau;\mathbf{x})=\sup_{\theta:\tau(\theta)=\tau}L(\theta;\mathbf{x}) be the likelihood induced on the transformed parameter. For any value \tau in the range of the transformation,

\sup_{\theta:\tau(\theta)=\tau}L(\theta;\mathbf{x}) \le\sup_{\theta}L(\theta;\mathbf{x}) =L(\widehat\theta;\mathbf{x}) =\sup_{\theta:\tau(\theta)=\tau(\widehat\theta)}L(\theta;\mathbf{x}) =M\bigl(\tau(\widehat\theta);\mathbf{x}\bigr).

The last equality holds because \widehat\theta itself satisfies the constraint \tau(\theta)=\tau(\widehat\theta), so the supremum over that constraint set is attained at \widehat\theta. Hence M(\tau(\widehat\theta))\ge M(\tau) for every \tau, and \tau(\widehat\theta) maximises M.

Apply the lemma with \theta=\Sigma^{-1}=:\Theta. The free coordinates are the k=p(p+1)/2 distinct entries of the precision, and the transformed coordinates are the distinct entries of \Sigma=\Theta^{-1}. Substitute \Sigma^{-1}=\Theta and \log|\Sigma|=-\log|\Theta|:

\ell(\boldsymbol{\mu},\Theta) =-\frac12 np\log(2\pi)+\frac12 n\log|\Theta| -\frac12\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})^\top\Theta(\mathbf{x}_i-\boldsymbol{\mu}).

The derivative with respect to a symmetric matrix uses \partial\log|\Theta|/\partial\Theta=\Theta^{-1} and \partial(\mathbf{v}^\top\Theta\mathbf{v})/\partial\Theta=\mathbf{v}\mathbf{v}^\top. Those two identities give

\frac{\partial\ell}{\partial\Theta} =\frac n2\,\Theta^{-1}-\frac12\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})(\mathbf{x}_i-\boldsymbol{\mu})^\top.

The mean has already been maximised. Set \boldsymbol{\mu}=\widehat{\boldsymbol{\mu}}=\bar{\mathbf{x}} and set the derivative equal to zero:

\frac n2\,\widehat\Theta^{-1}-\frac12\sum_{i=1}^n(\mathbf{x}_i-\bar{\mathbf{x}})(\mathbf{x}_i-\bar{\mathbf{x}})^\top=0.

Multiply by 2. The sum is A, so n\widehat\Theta^{-1}-A=0 and \widehat\Theta^{-1}=A/n. Invariance, applied to \Theta\mapsto\Theta^{-1}, says that the maximum likelihood estimator of \Sigma is

\widehat\Sigma=\frac1n A.

There is a condition behind that answer: A must be positive definite. Then strict concavity in the precision makes the stationary point the unique maximum. If A is singular, the unrestricted positive definite Gaussian covariance has no finite maximum likelihood estimator in the interior of the parameter space. Shrinking the variance in an unobserved direction increases the likelihood without bound. An algebraic expression such as A/n is not automatically an admissible parameter. When the estimator does exist, its divisor is n.

Bias and the sample covariance

Main article: Bias and the sample covariance

Whether \widehat\Sigma is unbiased is not settled by the likelihood equation. Expand the cross-product matrix about the true mean,

A=\sum_{i=1}^n(\mathbf{X}_i-\boldsymbol{\mu})(\mathbf{X}_i-\boldsymbol{\mu})^\top -n(\overline{\mathbf{X}}-\boldsymbol{\mu})(\overline{\mathbf{X}}-\boldsymbol{\mu})^\top.

If the observations are uncorrelated, with common mean \boldsymbol{\mu} and common covariance \Sigma, each outer product (\mathbf{X}_i-\boldsymbol{\mu})(\mathbf{X}_i-\boldsymbol{\mu})^\top has expectation \Sigma, so the first sum has expectation n\Sigma. The same uncorrelatedness gives \operatorname{Var}(\overline{\mathbf{X}})=\Sigma/n, and the subtracted term then has expectation \Sigma. Therefore

\mathbb{E}(A)=(n-1)\Sigma, \qquad \mathbb{E}(\widehat\Sigma)=\frac{n-1}{n}\Sigma.

Dividing A by n estimates \Sigma by (n-1)/n times \Sigma on average, so the maximum likelihood estimator is biased low. The unbiased matrix in common use divides by the degrees of freedom that remain after estimating the mean. It is the sample covariance matrix

S=\frac{1}{n-1}A.

The two estimates use the same measured scatter and differ only in their divisor. Dividing by n solves the likelihood problem. Dividing by n-1 corrects the expectation. Hotelling’s statistic below uses S=A/(n-1) because the stated F calibration is derived for that convention. Substituting a different divisor without adjusting the calibration changes the test.

A single mean is the wrong centre when the observations fall into G groups that may have different means and are assumed to share one covariance \Sigma. Let group g have n_g observations \mathbf{x}_{ig}, i=1,\ldots,n_g, with n=\sum_g n_g. Independence across groups lets the one-sample argument run inside each group, so the maximum likelihood estimator of \boldsymbol{\mu}_g is the group mean \bar{\mathbf{x}}_g. The unbiased estimator of the common \Sigma pools the corrected cross-product matrices. It is the pooled sample covariance

S=\frac{1}{n-G}\sum_{g=1}^G\sum_{i=1}^{n_g}(\mathbf{x}_{ig}-\bar{\mathbf{x}}_g)(\mathbf{x}_{ig}-\bar{\mathbf{x}}_g)^\top =\frac{1}{n-G}\sum_{g=1}^G A_g.

Each estimated group mean consumes one degree of freedom, which is why the divisor is n-G rather than n-1. The one-sample matrix is the case G=1.

Likelihood of a multivariate normal sample

The mean

The covariance and the invariance lemma

Bias and the sample covariance

Used in

See also