Sampling Distributions and Bias of Gaussian MLEs

Why Sampling Distributions Matter

A maximum-likelihood estimate is a random quantity because it depends on the observed sample. Its sampling distribution describes how it varies across repeated samples from the same population. For a multivariate Gaussian sample, the sample mean and centered scatter matrix have especially useful exact distributions: the mean is Gaussian, the scatter is Wishart, and the two are independent. These facts also give a direct finite-sample calculation of the covariance MLE’s bias.

Let $X_1,\ldots,X_n$ be independent $p$-dimensional random vectors with $X_i\sim N_p(\mu,\Sigma)$, where $\Sigma$ is positive definite. Assume $n>p$, so the unrestricted covariance MLE exists and the centered scatter matrix is positive definite with probability one. Define the sample mean and centered scatter matrix by

\bar X=\frac{1}{n}\sum_{i=1}^n X_i, \qquad S=\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^T.

The covariance maximum-likelihood estimator is $\hat\Sigma_{\mathrm{ML}}=S/n$. The denominator is $n$, not $n-1$; that distinction is central to its finite-sample bias.

The Wishart Distribution and Gaussian Sampling Theorem

Definition: If $Z_1,\ldots,Z_\nu$ are independent $N_p(0,\Sigma)$ vectors, then the matrix $\sum_{j=1}^{\nu}Z_jZ_j^T$ has a Wishart distribution with scale matrix $\Sigma$ and $\nu$ degrees of freedom, written $W\sim\mathcal W_p(\Sigma,\nu)$. For $\nu\ge p$, its density on positive definite matrices is $$ f_W(W)=\frac{|W|^{(\nu-p-1)/2}\exp\{-\tfrac12\operatorname{tr}(\Sigma^{-1}W)\}} {2^{\nu p/2}|\Sigma|^{\nu/2}\Gamma_p(\nu/2)}, $$ where $\Gamma_p(a)=\pi^{p(p-1)/4}\prod_{j=1}^p\Gamma(a+(1-j)/2)$.

Theorem: For independent $X_i\sim N_p(\mu,\Sigma)$, the sample mean and centered scatter satisfy $$ \bar X\sim N_p(\mu,\Sigma/n),\qquad S\sim\mathcal W_p(\Sigma,n-1), \qquad \bar X\ \text{is independent of }S. $$

The degrees of freedom are $n-1$ because centering removes the one common direction corresponding to the sample mean. Here is a proof that establishes the distributions and independence together.

Rotate the observation index.
Arrange the observations as the rows of an $n\times p$ matrix $X$. Choose an orthogonal $n\times n$ matrix $Q$ whose first column is $n^{-1/2}\mathbf 1$, and set $Y=Q^TX$. Orthogonal transformation preserves the joint Gaussian law and has absolute Jacobian one.

Identify the transformed rows.
The first row of $Y$ is $\sqrt n\,\bar X^T$. Its distribution is $N_p(\sqrt n\,\mu,\Sigma)$. Each of the remaining $n-1$ rows has distribution $N_p(0,\Sigma)$.

Use Gaussian independence.
The transformed rows are jointly Gaussian. Their cross-covariances are zero because $QQ^T=I_n$, so the rows are independent. In particular, the first row is independent of all remaining rows.

Express the scatter using the residual rows.
The first transformed row represents the common mean direction; the other rows span its orthogonal complement. Consequently, $S$ is the sum of the outer products of those $n-1$ residual rows. It therefore has distribution $\mathcal W_p(\Sigma,n-1)$, and is independent of $\bar X$.

For completeness, the scatter identity in the last step follows by writing $H=I_n-n^{-1}\mathbf 1\mathbf 1^T$, the centering projection. Then $S=X^THX$. The chosen $Q$ gives $Q^THQ=\operatorname{diag}(0,1,\ldots,1)$, so $S=\sum_{j=2}^n Y_{j\cdot}^TY_{j\cdot}$, where $Y_{j\cdot}$ is row $j$ of $Y$. This also makes clear why the first transformed row, which determines $\bar X$, does not enter $S$. $\square$

Expected Scatter and the Bias of the Covariance MLE

Theorem: Under the Gaussian sampling assumptions above, $$ E[S]=(n-1)\Sigma,\qquad E[\hat\Sigma_{\mathrm{ML}}]=\frac{n-1}{n}\Sigma, \qquad \operatorname{Bias}(\hat\Sigma_{\mathrm{ML}}) =E[\hat\Sigma_{\mathrm{ML}}]-\Sigma=-\frac{\Sigma}{n}. $$ Thus $S/(n-1)$ is unbiased for $\Sigma$, while the covariance MLE is biased downward in expectation.

To prove the expectation formula, use the Wishart representation $S=\sum_{j=1}^{n-1}Z_jZ_j^T$, where the $Z_j$ are independent $N_p(0,\Sigma)$. For each $j$, the $(a,b)$ entry of $Z_jZ_j^T$ is $Z_{ja}Z_{jb}$, whose expectation is $\Sigma_{ab}$. Therefore $E[S_{ab}]=\sum_{j=1}^{n-1}\Sigma_{ab}=(n-1)\Sigma_{ab}$ for every pair of indices, proving $E[S]=(n-1)\Sigma$. Dividing by $n$ gives $E[\hat\Sigma_{\mathrm{ML}}]=((n-1)/n)\Sigma$; subtracting $\Sigma$ gives the stated bias. $\square$

This is a matrix-valued bias statement: every entry has expectation $((n-1)/n)$ times the corresponding entry of $\Sigma$. The correction factor $n/(n-1)$ turns the MLE into the unbiased estimator $S/(n-1)$, but it does not make that corrected estimator the maximum-likelihood estimator.

Worked Example: Exact Distributions and Bias for a Two-Variable Sample

Suppose $p=2$, $n=4$, and the population parameters are $$ \mu=\begin{pmatrix}1\\-1\end{pmatrix}, \qquad \Sigma=\begin{pmatrix}2/3&1/6\\1/6&1/3\end{pmatrix}. $$ The determinant of $\Sigma$ is $(2/3)(1/3)-(1/6)^2=2/9-1/36=7/36>0$, so this is positive definite.

The sample mean has distribution $\bar X\sim N_2(\mu,\Sigma/4)$, where $$ \frac{\Sigma}{4} =\begin{pmatrix}1/6&1/24\\1/24&1/12\end{pmatrix}. $$ For example, the upper-left entry is $(2/3)/4=2/12=1/6$, and the off-diagonal entry is $(1/6)/4=1/24$.

The centered scatter has distribution $S\sim\mathcal W_2(\Sigma,3)$, independently of $\bar X$. Its expectation is $$ E[S]=3\Sigma =\begin{pmatrix}2&1/2\\1/2&1\end{pmatrix}. $$ Thus the covariance MLE has expectation $$ E[\hat\Sigma_{\mathrm{ML}}]=\frac{3}{4}\Sigma =\begin{pmatrix}1/2&1/8\\1/8&1/4\end{pmatrix}, $$ and bias $$ -\frac{\Sigma}{4} =\begin{pmatrix}-1/6&-1/24\\-1/24&-1/12\end{pmatrix}. $$ Each entry follows by dividing the corresponding entry of $\Sigma$ by four; in particular, $(2/3)/4=1/6$, not $2/3$.

Scalar Consequences and Further Calculations

The matrix result can be checked or interpreted one direction at a time. For any fixed nonzero vector $a\in\mathbb R^p$, the scalar projection $a^TX_i$ is normal with variance $a^T\Sigma a$. The corresponding projected scatter satisfies $$ \frac{a^TSa}{a^T\Sigma a}\sim\chi^2_{n-1}, $$ because it is the sum of squares of $n-1$ independent standard normal residual projections. This gives a familiar chi-square sampling statement without requiring the full matrix calculation.

Worked Example: Projecting the Scatter Onto a Contrast

Use the covariance matrix in the previous example and take $a=(1,-1)^T$. The population variance of this contrast is $$ a^T\Sigma a=\frac23+\frac13-2\left(\frac16\right)=1-\frac13=\frac23. $$ With $n=4$, it follows that $a^TSa/(2/3)\sim\chi^2_3$. Since a chi-square variable with three degrees of freedom has expectation three, $E[a^TSa]=3(2/3)=2$.

The MLE of the variance of this contrast is $a^T\hat\Sigma_{\mathrm{ML}}a=a^TSa/4$, so its expectation is $2/4=1/2$. The true projected variance is $2/3$, giving bias $1/2-2/3=-1/6$. This agrees with the matrix bias calculation: $$ a^T\left(-\frac{\Sigma}{4}\right)a=-\frac14\left(\frac23\right)=-\frac16. $$

Worked Example: Comparing MLE and Unbiased Scatter in Observed Data

Consider four observed vectors in two dimensions: $$ x_1=(1,0)^T,\quad x_2=(-1,0)^T,\quad x_3=(0,1)^T,\quad x_4=(0,-1)^T. $$ Their sample mean is $\bar x=(0,0)^T$. Since the mean is zero, the centered outer products are the ordinary outer products. Adding them gives $$ S=\sum_{i=1}^4 x_ix_i^T =\begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}0&0\\0&1\end{pmatrix} +\begin{pmatrix}0&0\\0&1\end{pmatrix} =\begin{pmatrix}2&0\\0&2\end{pmatrix}. $$

The covariance MLE divides by $n=4$, giving $\hat\Sigma_{\mathrm{ML}}=S/4=\operatorname{diag}(1/2,1/2)$. The unbiased sample covariance divides by $n-1=3$, giving $S/3=\operatorname{diag}(2/3,2/3)$. These are two different estimates computed from the same data. Unbiasedness describes the average over repeated samples, not whether a particular realized estimate lies above or below the unknown population covariance.

Common Pitfalls and Interpretation

Using the wrong denominator when identifying the MLE. The unrestricted Gaussian covariance MLE is $S/n$. The familiar correction $S/(n-1)$ gives an unbiased estimator, not the MLE.
Assuming independence without Gaussian sampling. For independent Gaussian observations, $\bar X$ and $S$ are independent. Their zero covariance alone would not establish independence for general, non-Gaussian observations.
Confusing the scatter matrix with a covariance estimate. $S$ has expectation $(n-1)\Sigma$; it is not itself centered to estimate $\Sigma$. The division by $n$ or $n-1$ determines which estimator is being considered.
Interpreting downward bias as a guarantee for each sample. The equality $E[\hat\Sigma_{\mathrm{ML}}]=(n-1)\Sigma/n$ is about repeated-sampling expectation. It does not say that every observed MLE is smaller than $\Sigma$ in every matrix direction.

The exact distributions are useful because they separate two sources of sampling variation: the location estimate is governed by a Gaussian distribution, while the residual variation is governed by a Wishart distribution. Their independence is also what makes many later Gaussian calculations tractable. The finite-sample bias is modest when $n$ is large, but the factor $(n-1)/n$ can matter when the sample size is limited.

Key takeaway: For an independent Gaussian sample, $\bar X$ and $S$ are independent, $\bar X\sim N_p(\mu,\Sigma/n)$, and $S\sim\mathcal W_p(\Sigma,n-1)$. Therefore the covariance MLE $S/n$ has bias $-\Sigma/n$, whereas $S/(n-1)$ is unbiased.

Check Your Understanding

Use the Gaussian sampling assumptions stated above unless a question says otherwise.

For a sample of size $n$, state the distributions of $\bar X$ and $S$, and say whether they are independent.
Why does the Wishart distribution of the centered scatter have $n-1$, rather than $n$, degrees of freedom?
If $E[S]=(n-1)\Sigma$, derive the bias of the covariance MLE $S/n$.
For a fixed vector $a$, what is the distribution of $a^TSa/(a^T\Sigma a)$?
Explain why $S/(n-1)$ can be unbiased without being the maximum-likelihood estimator.