Solution (source code)

= Solution

With the divisor-$n$ <sample covariance matrix> defined above, the <maximum-likelihood estimators> are
$$
\boxed{\widehat\mu=\bar x,\qquad \widehat V=S.}
$$
Their joint <sampling distribution> is specified by
$$
\boxed{\bar X\sim N_p(\mu,V/n),\qquad nS\sim W_p(n-1,V),
\qquad\bar X\ \text{and}\ S\ \text{are independent}.}
$$
Here $W_p(\nu,V)$ uses the <Wishart distribution> convention $\mathbb E W=\nu V$. Thus the covariance <maximum-likelihood estimator> has <expected value> $(n-1)V/n$; the <unbiased estimator> of the <covariance matrix> is $nS/(n-1)$. These distributions follow from <orthogonal projections> of the jointly <multivariate normal> observations onto the constant-observation direction and its <orthogonal complement>, though no proof is required here. The covariance maximum in the positive-definite parameter space requires $S$ nonsingular, which holds almost surely when $n>p$ and $V$ is positive definite. Otherwise the <rank of a centered sample covariance matrix> shows why the usual unrestricted finite maximum does not exist.