Tutorials › Multivariate Statistics (Graduate) › Deriving the Unrestricted Gaussian Maximum-Likelihood Estimates

Home  /  Tutorials  /  Gaussian Likelihood and Invariance

Gaussian Likelihood and Invariance · Tutorial 164 of 975

Deriving the Unrestricted Gaussian Maximum-Likelihood Estimates

Differentiate and profile the multivariate Gaussian log-likelihood to derive the sample mean and the covariance estimate with divisor n, then verify the maximum.

Intermediate 10 min read

What You'll Learn

  • Write the unrestricted multivariate Gaussian log-likelihood for independent observations.
  • Differentiate with respect to the mean vector and symmetric covariance matrix.
  • Derive the mean estimate using a residual sum-of-squares decomposition.
  • Derive the covariance estimate using the covariance score equation.
  • Verify the joint maximum when the sample covariance matrix is positive definite.
  • Distinguish the Gaussian maximum-likelihood covariance estimate from the unbiased sample covariance.

From the Gaussian Likelihood to Its Estimates

For independent Gaussian observations, estimating the mean and covariance together is a joint optimization problem: changing the mean changes the residuals that enter the covariance estimate. The likelihood has a simple solution, but deriving it requires keeping track of both parameters and checking that the stationary point is actually a maximum.

Let \(x_1,\ldots,x_n\) be observed vectors in \(\mathbb{R}^p\), from independent observations with distribution \(N_p(\mu,\Sigma)\). The parameter space is \(\mu\in\mathbb{R}^p\) and \(\Sigma\) symmetric positive definite. Define the sample mean and residual sum-of-products matrix by

$$ \bar{x}=\frac{1}{n}\sum_{i=1}^n x_i, \qquad Q(\mu)=\sum_{i=1}^n (x_i-\mu)(x_i-\mu)^\top. $$

The sample covariance matrix with divisor \(n\) is \(S_n=Q(\bar{x})/n\). We will first derive the estimates and then verify global maximization when \(S_n\) is positive definite. That condition ensures the estimate lies in the stated parameter space.

Definition: The unrestricted Gaussian maximum-likelihood estimates are the values of \(\mu\) and positive-definite \(\Sigma\) that maximize the likelihood over all allowed mean vectors and covariance matrices. When the maximum is attained at an interior parameter value, its log-likelihood derivatives with respect to both parameters must vanish.

The Log-Likelihood and Its Mean Derivative

The joint density of the independent observations gives the log-likelihood

$$ \ell(\mu,\Sigma) =-\frac{np}{2}\log(2\pi) -\frac{n}{2}\log|\Sigma| -\frac{1}{2}\operatorname{tr}\!\left(\Sigma^{-1}Q(\mu)\right). $$

Here \(|\Sigma|\) denotes the determinant, and the trace expression equals the sum of the squared residuals measured using \(\Sigma^{-1}\). To differentiate with respect to \(\mu\), use the residual decomposition

$$ Q(\mu)=Q(\bar{x})+n(\bar{x}-\mu)(\bar{x}-\mu)^\top. $$

To verify it, write \(x_i-\mu=(x_i-\bar{x})+(\bar{x}-\mu)\), expand the outer products, and sum. The cross terms vanish because \(\sum_i(x_i-\bar{x})=0\). This identity also shows directly how the mean affects the likelihood.

Result: For every fixed positive-definite \(\Sigma\), the log-likelihood is uniquely maximized over \(\mu\) at \(\bar{x}\). Indeed, the decomposition gives
$$ \ell(\mu,\Sigma) = \ell(\bar{x},\Sigma) -\frac{n}{2}(\mu-\bar{x})^\top\Sigma^{-1}(\mu-\bar{x}). $$

Since \(\Sigma^{-1}\) is positive definite, the quadratic form is nonnegative and equals zero only when \(\mu=\bar{x}\). Equivalently, differentiating gives the mean score

$$ \nabla_\mu\ell(\mu,\Sigma) =\Sigma^{-1}\sum_{i=1}^n(x_i-\mu) =n\Sigma^{-1}(\bar{x}-\mu). $$

Thus the mean score vanishes exactly at \(\hat{\mu}=\bar{x}\), for every fixed \(\Sigma\). The decomposition is stronger than the score equation alone: it proves that this choice maximizes the likelihood in the mean direction.

Worked Example: The Mean Score in Two Dimensions

Consider the four observations \(x_1=(1,0)^\top\), \(x_2=(3,0)^\top\), \(x_3=(1,2)^\top\), and \(x_4=(3,2)^\top\). Their sample mean is

$$ \bar{x} =\frac{1}{4} \begin{pmatrix}8\\4\end{pmatrix} = \begin{pmatrix}2\\1\end{pmatrix}. $$

For any fixed positive-definite \(\Sigma\), the mean score is \(4\Sigma^{-1}(\bar{x}-\mu)\), so it is zero at \(\mu=(2,1)^\top\). If instead \(\mu=(1,1)^\top\), the log-likelihood difference from using the sample mean is

$$ \ell((1,1)^\top,\Sigma)-\ell((2,1)^\top,\Sigma) =-2 \begin{pmatrix}-1&0\end{pmatrix} \Sigma^{-1} \begin{pmatrix}-1\\0\end{pmatrix}<0. $$

The strict inequality follows from positive definiteness. Thus the sample mean is not merely a solution to the score equation; it gives a strictly higher likelihood than this alternative mean.

The Covariance Derivative

Now hold \(\mu\) fixed and differentiate with respect to the symmetric matrix \(\Sigma\). For a positive-definite matrix, the differential identities are \(d\log|\Sigma|=\operatorname{tr}(\Sigma^{-1}d\Sigma)\) and \(d\Sigma^{-1}=-\Sigma^{-1}(d\Sigma)\Sigma^{-1}\). Applying them to the log-likelihood yields

$$ d\ell = \frac{1}{2}\operatorname{tr}\!\left[ \left\{ \Sigma^{-1}Q(\mu)\Sigma^{-1}-n\Sigma^{-1} \right\}d\Sigma \right]. $$

The matrix in braces is symmetric. Since \(d\Sigma\) may be any symmetric perturbation at an interior maximum, the covariance score vanishes precisely when

$$ \Sigma^{-1}Q(\mu)\Sigma^{-1}-n\Sigma^{-1}=0. $$

Multiplying on the left and right by \(\Sigma\) gives \(Q(\mu)=n\Sigma\). Combining this equation with the mean estimate gives the joint stationary point

$$ \hat{\mu}=\bar{x}, \qquad \hat{\Sigma}=\frac{1}{n}Q(\bar{x})=S_n. $$

A stationary point need not be a maximum. The following argument verifies the global maximum when \(S_n\) is positive definite.

Theorem: If \(S_n\) is positive definite, the unrestricted Gaussian likelihood over \(\mu\in\mathbb{R}^p\) and positive-definite \(\Sigma\) has a unique global maximum at \(\hat{\mu}=\bar{x}\) and \(\hat{\Sigma}=S_n\).

Proof. For each \(\Sigma\), the mean decomposition proves \(\ell(\mu,\Sigma)\leq\ell(\bar{x},\Sigma)\), with equality only for \(\mu=\bar{x}\). It remains to maximize the log-likelihood with \(\mu=\bar{x}\). Write \(S=S_n\), which is positive definite by hypothesis, and let \(C=S^{-1/2}\Sigma S^{-1/2}\). This is positive definite; denote its eigenvalues by \(\lambda_1,\ldots,\lambda_p>0\). Since \(Q(\bar{x})=nS\), the part of the negative twice log-likelihood that depends on \(\Sigma\), after removing constants, is

$$ n\log|\Sigma|+n\operatorname{tr}(\Sigma^{-1}S) = n\log|S|+n\sum_{j=1}^p\left(\log\lambda_j+\frac{1}{\lambda_j}\right). $$

For \(h(\lambda)=\log\lambda+1/\lambda\), the derivative is \(h'(\lambda)=(\lambda-1)/\lambda^2\). It is negative for \(0<\lambda<1\) and positive for \(\lambda>1\), so \(h\) has its unique global minimum at \(\lambda=1\). Each term in the sum is therefore minimized exactly when \(\lambda_j=1\). A positive-definite matrix with all eigenvalues equal to one is the identity, so \(C=I\), which is equivalent to \(\Sigma=S\). This gives a unique global maximum in the covariance parameter, and the mean inequality then establishes the unique joint maximum. \(\square\)

Worked Example: A Covariance Estimate With Nonzero Off-Diagonal Entries

Take \(n=3\) observations in two dimensions: \(x_1=(1,0)^\top\), \(x_2=(0,1)^\top\), and \(x_3=(-1,-1)^\top\). Their mean is zero, so

$$ Q(\bar{x}) = \begin{pmatrix}1&0\\0&0\end{pmatrix} + \begin{pmatrix}0&0\\0&1\end{pmatrix} + \begin{pmatrix}1&1\\1&1\end{pmatrix} = \begin{pmatrix}2&1\\1&2\end{pmatrix}. $$

The unrestricted covariance estimate is therefore

$$ \hat{\Sigma} =\frac{1}{3} \begin{pmatrix}2&1\\1&2\end{pmatrix}. $$

Its determinant is \((4-1)/9=1/3>0\), and its diagonal entries are positive, so it is positive definite. The covariance score equation checks directly: \(Q(\bar{x})=3\hat{\Sigma}\). The profiled log-likelihood at this estimate is

$$ \ell(\bar{x},\hat{\Sigma}) =-3\log(2\pi)-\frac{3}{2}\log(1/3)-3 =-3\log(2\pi)+\frac{3}{2}\log 3-3. $$

The off-diagonal entries are not set to zero: the unrestricted model estimates the covariance between the two coordinates as \(1/3\).

Checking the Joint Estimate in a Second Data Set

Worked Example: Verifying Both Estimates and the Likelihood Value

Let the observations be \(x_1=(0,0)^\top\), \(x_2=(2,0)^\top\), and \(x_3=(1,3)^\top\). Their sample mean is \(\bar{x}=(1,1)^\top\). The residual vectors are \((-1,-1)^\top\), \((1,-1)^\top\), and \((0,2)^\top\), giving

$$ Q(\bar{x}) = \begin{pmatrix}2&0\\0&6\end{pmatrix}, \qquad S_n=\hat{\Sigma} = \begin{pmatrix}2/3&0\\0&2\end{pmatrix}. $$

The estimate is positive definite because both diagonal entries are positive. At this estimate, \(\hat{\Sigma}^{-1}=\operatorname{diag}(3/2,1/2)\), and \(\operatorname{tr}(\hat{\Sigma}^{-1}Q(\bar{x}))=3+3=6=np\), as required by the covariance score equation. If we keep \(\Sigma=\hat{\Sigma}\) but replace the mean by \(\mu=(0,0)^\top\), then \(\bar{x}-\mu=(1,1)^\top\) and

$$ \ell(\mu,\hat{\Sigma})-\ell(\bar{x},\hat{\Sigma}) =-\frac{3}{2} \left(\frac{3}{2}+\frac{1}{2}\right) =-3. $$

Thus this alternative mean has a lower log-likelihood by 3. At the joint estimate, the likelihood is

$$ \ell(\bar{x},\hat{\Sigma}) =-3\log(2\pi)-\frac{3}{2}\log(4/3)-3, $$

because \(|\hat{\Sigma}|=4/3\). This agrees with the general expression obtained by substituting \(Q(\bar{x})=n\hat{\Sigma}\) into the log-likelihood.

Why the Divisor Is \(n\), Not \(n-1\)

The likelihood estimate is \(Q(\bar{x})/n\). The divisor follows from the covariance score equation \(Q(\bar{x})=n\Sigma\); it is not chosen by an unbiasedness requirement. Under the usual independent Gaussian sampling model, the sample covariance with divisor \(n-1\) is unbiased for \(\Sigma\), but that is a different property and a different estimator. The distinction matters whenever a likelihood calculation, likelihood comparison, or fitted Gaussian model is required.

The global-maximization proof above assumes \(S_n\) is positive definite. If it is not, the displayed stationary covariance is not an allowed positive-definite parameter value, so the theorem does not establish an attained maximum. The conditions and behavior in that case are a separate question; the important point here is not to treat a formal score solution as an MLE without checking the parameter space.

Common pitfall: Dividing the residual sum-of-products matrix by \(n-1\) because it is an unbiased covariance estimate. That argument answers an unbiasedness question, not a likelihood maximization question. For unrestricted Gaussian maximum likelihood, the covariance score gives \(Q(\bar{x})=n\hat{\Sigma}\), hence the divisor is \(n\). A second pitfall is to stop after solving the score equations: a stationary point must still be shown to lie in the parameter space and to maximize the likelihood.

Check Your Understanding

Use the likelihood and verification arguments above to answer these questions.

  1. For fixed positive-definite \(\Sigma\), why is \(\bar{x}\) the unique maximizing value of \(\mu\)?
  2. Write the covariance score equation and show how it leads to \(Q(\mu)=n\Sigma\).
  3. Why does the unrestricted Gaussian covariance MLE use divisor \(n\), rather than \(n-1\)?
  4. In the global-maximization proof, what property of \(h(\lambda)=\log\lambda+1/\lambda\) identifies the maximizing covariance?
  5. What condition on \(S_n\) allows the theorem to conclude that the displayed covariance estimate is an allowed parameter value and a global maximum?

Previous

Next

Back to Tutorials
Clinical Biostats