Profile Likelihood for an Unknown-Covariance Mean

Why Profile the Covariance?

In a multivariate normal model, evidence against a specified mean depends on how far the observed sample mean lies from that value relative to the observed covariance. When the covariance matrix is unknown, a likelihood-ratio test must account for its estimation rather than treating it as fixed. Profiling provides a direct way to do this: for each candidate mean, maximize the likelihood over every positive-definite covariance matrix, then compare the resulting maximized likelihoods.

The calculation connects the likelihood-ratio test to Hotelling’s T-squared. It also makes clear that the likelihood ratio is not itself T-squared: it is a decreasing function of it, with a finite-sample form that differs from the usual large-sample chi-square approximation.

Model, Notation, and Profile Likelihood

Let $X_1,\ldots,X_n$ be independent observations from $N_p(\mu,\Sigma)$, where $\mu\in\mathbb{R}^p$ and $\Sigma$ is positive definite. We test $H_0:\mu=\mu_0$ against the unrestricted alternative. Assume $n>p$, so the sample covariance is positive definite with probability one under this model. Write

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

For observed vectors $x_1,\ldots,x_n$, define $A=\sum_{i=1}^n(x_i-\bar x)(x_i-\bar x)'=(n-1)S$. For any candidate mean $\mu$, let $B(\mu)=\sum_{i=1}^n(x_i-\mu)(x_i-\mu)'$. The likelihood, as a function of both parameters, is

L(\mu,\Sigma) =(2\pi)^{-np/2}|\Sigma|^{-n/2} \exp\left\{-\frac12\operatorname{tr}\bigl(\Sigma^{-1}B(\mu)\bigr)\right\}.

Definition: The profile likelihood for $\mu$ is $L_p(\mu)=\sup_{\Sigma\succ0}L(\mu,\Sigma)$, where the supremum is over all positive-definite covariance matrices. The likelihood ratio for $H_0:\mu=\mu_0$ is $\Lambda=L_p(\mu_0)/\sup_{\mu}L_p(\mu)$. Thus $0<\Lambda\leq1$, and smaller values are stronger evidence against the null.

Profiling the Covariance Matrix

For a fixed $\mu$, the Gaussian log-likelihood, up to a constant independent of the parameters, is

\ell(\mu,\Sigma) =-\frac n2\log|\Sigma| -\frac12\operatorname{tr}\bigl(\Sigma^{-1}B(\mu)\bigr)+\text{constant}.

When $B(\mu)$ is positive definite, the maximizing covariance is $B(\mu)/n$. To verify this, set $B=B(\mu)$ and $C=B^{1/2}\Sigma^{-1}B^{1/2}$. Then $C$ ranges over all positive-definite matrices as $\Sigma$ does, and

-\frac n2\log|\Sigma|-\frac12\operatorname{tr}(\Sigma^{-1}B) =-\frac n2\log|B|+\frac n2\log|C|-\frac12\operatorname{tr}(C).

If $\lambda_1,\ldots,\lambda_p$ are the positive eigenvalues of $C$, the terms depending on $C$ are $\sum_{j=1}^p\{(n/2)\log\lambda_j-\lambda_j/2\}$. Each summand has derivative $n/(2\lambda_j)-1/2$, which is positive for $\lambda_j<n$, zero at $n$, and negative for $\lambda_j>n$. Its unique maximum is therefore at $\lambda_j=n$. The maximum occurs at $C=nI$, equivalent to $\hat\Sigma_\mu=B(\mu)/n$. Substitution gives

\ell_p(\mu) =-\frac{np}{2}\{\log(2\pi)+1\} -\frac n2\log|B(\mu)|+\frac{np}{2}\log n.

The terms not involving $\mu$ will cancel in the likelihood ratio. The unrestricted maximum over $\mu$ occurs at $\bar x$: the decomposition

B(\mu)=A+n(\bar x-\mu)(\bar x-\mu)'

shows, by the determinant lemma, that $|B(\mu)|\geq |A|$, with equality at $\mu=\bar x$. The unrestricted estimates are consequently $\hat\mu=\bar x$ and $\hat\Sigma=A/n=(n-1)S/n$. Notice that this covariance maximum-likelihood estimate uses divisor $n$ relative to the centered scatter matrix, whereas $S$ uses divisor $n-1$.

The Likelihood Ratio in Terms of Hotelling’s T-Squared

Theorem: Under the model and nonsingularity assumptions above, define $T^2=n(\bar x-\mu_0)'S^{-1}(\bar x-\mu_0)$. The likelihood ratio for testing $H_0:\mu=\mu_0$ with unrestricted $\Sigma$ is $$ \Lambda=\left(1+\frac{T^2}{n-1}\right)^{-n/2}, \qquad -2\log\Lambda=n\log\left(1+\frac{T^2}{n-1}\right). $$

Proof. Put $d=\bar x-\mu_0$. Under the null, $B(\mu_0)=A+ndd'$, while at the unrestricted maximum $B(\bar x)=A$. Since $A=(n-1)S$, the matrix determinant lemma, $|A+uv'|=|A|(1+v'A^{-1}u)$, gives

\frac{|B(\mu_0)|}{|A|} =1+nd'A^{-1}d =1+\frac{n}{n-1}d'S^{-1}d =1+\frac{T^2}{n-1}.

The profile log-likelihood differs between the null and unrestricted fits only through the log determinants. Therefore

\Lambda =\exp\{\ell_p(\mu_0)-\ell_p(\bar x)\} =\left(\frac{|A|}{|B(\mu_0)|}\right)^{n/2} =\left(1+\frac{T^2}{n-1}\right)^{-n/2}.

Taking $-2$ times the logarithm yields the stated statistic. $\square$

This result shows that the likelihood-ratio statistic is a strictly increasing function of $T^2$: as $T^2$ increases, $-2\log\Lambda$ increases and $\Lambda$ decreases. The likelihood-ratio test and a test that rejects for large Hotelling’s T-squared therefore rank the data in exactly the same way.

Worked Example: Computing the Likelihood Ratio from Summary Statistics

Suppose $n=8$, $p=2$, the null mean is $\mu_0=(0,0)'$, the observed mean is $\bar x=(1,2)'$, and $S=\begin{pmatrix}2&0\\0&8\end{pmatrix}$. Compute $T^2$, $\Lambda$, and $-2\log\Lambda$.

Here $d=(1,2)'$ and $S^{-1}=\begin{pmatrix}1/2&0\\0&1/8\end{pmatrix}$. Hence $d'S^{-1}d=1/2+4/8=1$, so $T^2=8(1)=8$. The determinant factor is $1+T^2/(n-1)=1+8/7=15/7$. It follows that

\Lambda=\left(\frac{15}{7}\right)^{-4} =\left(\frac{7}{15}\right)^4 =\frac{2401}{50625}\approx0.04743, \qquad -2\log\Lambda=8\log(15/7)\approx6.097.

The likelihood under the null is about 4.7% of the unrestricted maximized likelihood. This ratio alone is not a p-value; a reference distribution is needed to calibrate a test.

Worked Example: Checking the Determinant Calculation with Explicit Data

Consider $n=5$, $p=2$, and the five observations $(1+2/\sqrt3,\sqrt{2/3})'$, $(1-2/\sqrt3,\sqrt{2/3})'$, $(1,-2\sqrt{2/3})'$, $(1,0)'$, and $(1,0)'$. Test $H_0:\mu=(0,0)'$.

The first coordinates sum to $5$, and the second coordinates sum to zero, so $\bar x=(1,0)'$. The first-coordinate residuals are $2/\sqrt3,-2/\sqrt3,0,0,0$; their squared sum is $8/3$. The second-coordinate residuals are $\sqrt{2/3},\sqrt{2/3},-2\sqrt{2/3},0,0$; their squared sum is $4$. The cross-product sum is $(2/\sqrt3)\sqrt{2/3}+(-2/\sqrt3)\sqrt{2/3}=0$. Thus $A=\operatorname{diag}(8/3,4)$ and $S=A/4=\operatorname{diag}(2/3,1)$.

With $d=(1,0)'$, $d'S^{-1}d=3/2$, giving $T^2=5(3/2)=15/2$. Therefore the determinant ratio is $1+(15/2)/4=23/8$, and

\Lambda=\left(\frac{23}{8}\right)^{-5/2} \approx0.07136, \qquad -2\log\Lambda=5\log(23/8)\approx5.280.

As a direct check, $|A|=(8/3)(4)=32/3$. The null scatter matrix is $B(\mu_0)=A+5dd'=\operatorname{diag}(23/3,4)$, so $|B(\mu_0)|=92/3$. Their ratio is $(92/3)/(32/3)=23/8$, as required. The reciprocal likelihood ratio is about $14.014$, so its reciprocal is about $0.07136$.

Finite-Sample Calibration and Interpretation

Under $H_0$ and the multivariate normal model, Hotelling’s exact one-sample result is $(n-p)T^2/\{p(n-1)\}\sim F_{p,n-p}$. Since $-2\log\Lambda$ is an increasing function of $T^2$, an exact finite-sample likelihood-ratio test can use this F calibration: reject for sufficiently large $T^2$, equivalently for sufficiently small $\Lambda$. The likelihood-ratio calculation does not require a new reference distribution; it expresses the same ordering of outcomes in likelihood terms.

Wilks’ theorem gives a different, asymptotic perspective. With fixed $p$, regularity conditions satisfied, and the null model correctly specified, $-2\log\Lambda$ converges in distribution to $\chi^2_p$ as $n$ grows. The exact F calibration is available here because of the Gaussian model; the chi-square approximation is not the exact finite-sample law of this statistic.

Worked Example: Converting a Larger T-Squared Value

Suppose $n=10$, $p=3$, and the observed value is $T^2=12$. The likelihood-ratio statistic is

-2\log\Lambda =10\log\left(1+\frac{12}{9}\right) =10\log(7/3)\approx8.473, \qquad \Lambda=(7/3)^{-5}=(3/7)^5\approx0.01446.

For the exact F calibration, the corresponding statistic is $(10-3)12/\{3(10-1)\}=84/27=28/9$. Whether this is significant depends on the chosen test level and the $F_{3,7}$ reference distribution. The likelihood ratio $0.01446$ is not itself a 1.446% p-value.

Common Pitfall: Mixing Up the Ratio, Its Log, and the F Statistic

A frequent error is to identify $-2\log\Lambda$ with $T^2$, or to treat $\Lambda$ as a p-value. Neither is correct. The exact relationship is $-2\log\Lambda=n\log\{1+T^2/(n-1)\}$, and a p-value is a tail probability under a specified null distribution. Another error is to forget that the exponent in the likelihood ratio is $n/2$, not $1/2$: the sample contributes $n$ Gaussian density factors. Finally, replacing $S$ by the maximum-likelihood covariance estimate without tracking its divisor changes the determinant factor. Here $S$ uses $n-1$, and the unrestricted maximum-likelihood estimate is $(n-1)S/n$.

Key takeaway: With an unrestricted unknown covariance matrix, profiling the Gaussian likelihood yields $\Lambda=\{1+T^2/(n-1)\}^{-n/2}$. This is a monotone transformation of Hotelling’s T-squared, so both tests give the same rejection ordering when calibrated consistently.

Check Your Understanding

Use the assumptions and notation defined above.

For fixed $\mu$, what covariance matrix maximizes the Gaussian likelihood, expressed using $B(\mu)$?
Why does the unrestricted fit use $\hat\mu=\bar x$, and what is the corresponding covariance maximum-likelihood estimate in terms of $S$?
If $T^2=0$, what are $\Lambda$ and $-2\log\Lambda$? Explain what this says about the fitted mean.
For fixed $n$, is $\Lambda$ increasing or decreasing in $T^2$? What does that imply about rejection regions?
Why is a likelihood ratio such as $0.05$ not automatically a p-value?