A random vector is a column of \(p\) random variables measured on the same unit:
\[ \mathbf{X} = \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_p \end{pmatrix}. \]Its mean vector is taken componentwise, and its dispersion matrix (also called the variance–covariance matrix) collects every variance and every covariance:
\[ \boldsymbol\mu = E(\mathbf{X}) = \begin{pmatrix} E(X_1) \\ \vdots \\ E(X_p)\end{pmatrix}, \qquad \boldsymbol\Sigma = E\!\left[(\mathbf{X} - \boldsymbol\mu)(\mathbf{X} - \boldsymbol\mu)'\right], \]so that \(\sigma_{ii} = \operatorname{Var}(X_i)\) and \(\sigma_{ij} = \operatorname{Cov}(X_i, X_j)\). Three properties do most of the work, and each is one line from the definition.
(i) Linear transformation. For a constant \(q \times p\) matrix \(\mathbf{A}\) and constant vector \(\mathbf{b}\), put \(\mathbf{Y} = \mathbf{AX} + \mathbf{b}\). Then
\[ E(\mathbf{Y}) = \mathbf{A}\boldsymbol\mu + \mathbf{b}, \qquad \operatorname{Cov}(\mathbf{Y}) = \mathbf{A}\boldsymbol\Sigma\mathbf{A}'. \]Proof of the second. By definition \(\operatorname{Cov}(\mathbf{Y}) = E[(\mathbf{Y} - E\mathbf{Y})(\mathbf{Y} - E\mathbf{Y})']\). Substituting \(\mathbf{Y} - E\mathbf{Y} = \mathbf{A}(\mathbf{X} - \boldsymbol\mu)\) gives \(E[\mathbf{A}(\mathbf{X}-\boldsymbol\mu)(\mathbf{X}-\boldsymbol\mu)'\mathbf{A}']\). The matrices \(\mathbf{A}\) and \(\mathbf{A}'\) are constants, so they come outside the expectation, leaving \(\mathbf{A}\,E[(\mathbf{X}-\boldsymbol\mu)(\mathbf{X}-\boldsymbol\mu)']\,\mathbf{A}' = \mathbf{A}\boldsymbol\Sigma\mathbf{A}'\). The additive constant \(\mathbf{b}\) cancels in the centring and therefore never appears. \(\blacksquare\)
(ii) Non-negative definiteness. For any constant vector \(\mathbf{a}\), rule (i) with \(\mathbf{A} = \mathbf{a}'\) gives \(\operatorname{Var}(\mathbf{a}'\mathbf{X}) = \mathbf{a}'\boldsymbol\Sigma\mathbf{a}\). A variance cannot be negative, so \(\mathbf{a}'\boldsymbol\Sigma\mathbf{a} \ge 0\) for every \(\mathbf{a}\): every dispersion matrix is non-negative definite. It is positive definite unless some linear combination \(\mathbf{a}'\mathbf{X}\) has zero variance — that is, unless the variables satisfy an exact linear relation.
(iii) Standardising. With \(\mathbf{D} = \operatorname{diag}(\sigma_{11}, \ldots, \sigma_{pp})\), the correlation matrix is \(\boldsymbol\rho = \mathbf{D}^{-1/2}\boldsymbol\Sigma\mathbf{D}^{-1/2}\), which is rule (i) applied to \(\mathbf{D}^{-1/2}\mathbf{X}\). So results proved for \(\boldsymbol\Sigma\) transfer to \(\boldsymbol\rho\) without repetition.
Repeat an experiment \(n\) times independently. Each trial falls into exactly one of \(k\) categories, with probabilities \(p_1, \ldots, p_k\) summing to 1. Let \(X_j\) count the trials landing in category \(j\). Then \((X_1, \ldots, X_k)\) has the multinomial distribution
\[ P(X_1 = x_1, \ldots, X_k = x_k) = \frac{n!}{x_1!\,x_2!\cdots x_k!}\; p_1^{x_1}p_2^{x_2}\cdots p_k^{x_k}, \qquad \sum_j x_j = n. \]It is the direct generalisation of the binomial, which is the case \(k = 2\). The multinomial coefficient counts the arrangements of a sequence of \(n\) outcomes containing \(x_1\) of the first kind, \(x_2\) of the second, and so on; each such arrangement has the same probability \(\prod_j p_j^{x_j}\) by independence.
Write \(X_j = \sum_{t=1}^{n} I_{jt}\), where \(I_{jt}\) is 1 if trial \(t\) falls in category \(j\) and 0 otherwise. Each \(I_{jt}\) is a Bernoulli variable with success probability \(p_j\), so
\[ E(X_j) = np_j, \qquad \operatorname{Var}(X_j) = np_j(1 - p_j). \]For the covariance, take \(i \ne j\). Within one trial \(I_{it}I_{jt} = 0\), because a trial cannot be in two categories at once, so
\[ E(I_{it}I_{jt}) = 0 \quad\Longrightarrow\quad \operatorname{Cov}(I_{it}, I_{jt}) = 0 - p_ip_j = -p_ip_j. \]Different trials are independent and contribute nothing, so summing over the \(n\) trials gives
\[ \operatorname{Cov}(X_i, X_j) = -np_ip_j \qquad (i \ne j). \]Read the sign. The covariance is negative for every pair, with no exceptions, and the reason is structural rather than statistical: the counts must add to \(n\), so a trial gained by one category is a trial lost to the others. The correlation follows at once:
\[ \rho_{ij} = \frac{-np_ip_j}{\sqrt{np_i(1-p_i)}\sqrt{np_j(1-p_j)}} = -\sqrt{\frac{p_ip_j}{(1-p_i)(1-p_j)}}. \]Notice that \(n\) has cancelled: the strength of the dependence does not change with sample size, only the counts do.
The distribution is singular. Because \(\sum_j X_j = n\) identically, the vector \(\mathbf{a} = (1, 1, \ldots, 1)'\) gives \(\operatorname{Var}(\mathbf{a}'\mathbf{X}) = \operatorname{Var}(n) = 0\), so by rule (ii) above \(\boldsymbol\Sigma\) is non-negative definite but not positive definite: its determinant is 0 and its rank is \(k - 1\). That is why the multinomial is normally written with \(k - 1\) free counts.
Given. \(n = 10\) independent trials with three outcomes of probabilities \(p_1 = 0.5\), \(p_2 = 0.3\), \(p_3 = 0.2\). Asked. The probability of the split \((5, 3, 2)\), and the covariance matrix of the counts.
Step 1 — the multinomial coefficient.
\[ \frac{10!}{5!\,3!\,2!} = \frac{3{,}628{,}800}{120 \times 6 \times 2} = \frac{3{,}628{,}800}{1440} = 2520. \]Step 2 — the three powers.
\[ (0.5)^{5} = 0.03125, \qquad (0.3)^{3} = 0.027000, \qquad (0.2)^{2} = 0.0400, \] \[ 0.03125 \times 0.027000 \times 0.0400 = 0.00003375. \]Step 3 — multiply.
\[ P(5, 3, 2) = 2520 \times 0.00003375 = 0.085050. \]Step 4 — the variances. \(np_j(1-p_j)\) for each \(j\):
\[ 10(0.5)(0.5) = 2.5, \qquad 10(0.3)(0.7) = 2.1, \qquad 10(0.2)(0.8) = 1.6. \]Step 5 — the covariances. \(-np_ip_j\):
\[ -10(0.5)(0.3) = -1.5, \qquad -10(0.5)(0.2) = -1.0, \qquad -10(0.3)(0.2) = -0.6. \]Step 6 — assemble and test.
\[ \boldsymbol\Sigma = \begin{pmatrix} 2.5 & -1.5 & -1.0 \\ -1.5 & 2.1 & -0.6 \\ -1.0 & -0.6 & 1.6 \end{pmatrix}. \]Every row sums to zero — \(2.5 - 1.5 - 1.0 = 0\), \(-1.5 + 2.1 - 0.6 = 0\), \(-1.0 - 0.6 + 1.6 = 0\) — which is the singularity of Step 6 made visible: \(\boldsymbol\Sigma\mathbf{1} = \mathbf{0}\), so \(\det\boldsymbol\Sigma = 0\) and the rank is 2, not 3.
Step 7 — a correlation, computed twice. From the matrix,
\[ \rho_{12} = \frac{-1.5}{\sqrt{2.5 \times 2.1}} = \frac{-1.5}{\sqrt{5.25}} = \frac{-1.5}{2.291288} = -0.654654. \]From the general formula, with no reference to \(n\),
\[ -\sqrt{\frac{0.5 \times 0.3}{0.5 \times 0.7}} = -\sqrt{\frac{0.15}{0.35}} = -\sqrt{0.428571} = -0.654654. \checkmark \]Interpretation. Categories 1 and 2 are strongly negatively correlated, and would be so at \(n = 10\) or \(n = 10{,}000\). This is the simplest case of a theme that runs through the whole course: a dispersion matrix can be singular, and when it is, the singularity is telling you about a constraint in the data rather than about a numerical accident.
Collapsing categories preserves the family. If the \(k\) categories are grouped into \(m\) groups, the group counts are multinomial with the group probabilities — because a trial falls in a group exactly when it falls in one of its categories, and those events are disjoint. In particular each single \(X_j\) is binomial \((n, p_j)\): group category \(j\) against all the rest.
Conditioning also preserves it. Given \(X_k = x_k\), the remaining \(n - x_k\) trials are distributed among the first \(k-1\) categories with the renormalised probabilities \(p_j/(1 - p_k)\), so
\[ (X_1, \ldots, X_{k-1}) \mid X_k = x_k \;\sim\; \text{Multinomial}\!\left(n - x_k,\ \frac{p_1}{1-p_k}, \ldots, \frac{p_{k-1}}{1-p_k}\right). \]This closure under both operations is exactly the behaviour the multivariate normal will show in the next two sections, and it is what makes a distribution usable in \(p\) dimensions.
\(\mathbf{X}\) has the \(p\)-variate normal distribution \(N_p(\boldsymbol\mu, \boldsymbol\Sigma)\), with \(\boldsymbol\Sigma\) positive definite, if
\[ f(\mathbf{x}) = \frac{1}{(2\pi)^{p/2}\,|\boldsymbol\Sigma|^{1/2}} \exp\!\left\{-\tfrac12 (\mathbf{x} - \boldsymbol\mu)'\boldsymbol\Sigma^{-1} (\mathbf{x} - \boldsymbol\mu)\right\}, \qquad \mathbf{x} \in \mathbb{R}^{p}. \]Read it one factor at a time, because every factor is the univariate density in disguise:
The moment generating function is the single most useful fact about the family:
\[ M_{\mathbf{X}}(\mathbf{t}) = E\!\left(e^{\mathbf{t}'\mathbf{X}}\right) = \exp\!\left(\mathbf{t}'\boldsymbol\mu + \tfrac12\,\mathbf{t}'\boldsymbol\Sigma\mathbf{t}\right). \](a) Every linear combination is normal. Let \(\mathbf{Y} = \mathbf{AX} + \mathbf{b}\) with \(\mathbf{A}\) a \(q \times p\) constant matrix. Then
\[ M_{\mathbf{Y}}(\mathbf{t}) = E\!\left(e^{\mathbf{t}'(\mathbf{AX}+\mathbf{b})}\right) = e^{\mathbf{t}'\mathbf{b}}\,M_{\mathbf{X}}(\mathbf{A}'\mathbf{t}) = \exp\!\left(\mathbf{t}'(\mathbf{A}\boldsymbol\mu + \mathbf{b}) + \tfrac12\,\mathbf{t}'\mathbf{A}\boldsymbol\Sigma\mathbf{A}'\mathbf{t}\right), \]which is the MGF of \(N_q(\mathbf{A}\boldsymbol\mu + \mathbf{b}, \mathbf{A}\boldsymbol\Sigma\mathbf{A}')\). By the uniqueness of moment generating functions, that is the distribution of \(\mathbf{Y}\). Taking \(\mathbf{A}\) to be a single row shows that \(\mathbf{a}'\mathbf{X}\) is univariate normal for every \(\mathbf{a}\) — a property which is often used as the definition, because it also covers the singular case.
(b) Every marginal is normal. Choose \(\mathbf{A}\) to select the rows you want. Picking out \(X_1\) alone gives \(X_1 \sim N(\mu_1, \sigma_{11})\) at once.
(c) For the normal, zero covariance is independence. Partition \(\mathbf{X}' = (\mathbf{X}_1', \mathbf{X}_2')\) and suppose \(\boldsymbol\Sigma_{12} = \mathbf{0}\). Then the quadratic form in the MGF splits,
\[ \mathbf{t}'\boldsymbol\Sigma\mathbf{t} = \mathbf{t}_1'\boldsymbol\Sigma_{11}\mathbf{t}_1 + \mathbf{t}_2'\boldsymbol\Sigma_{22}\mathbf{t}_2, \]so \(M_{\mathbf{X}}(\mathbf{t}) = M_{\mathbf{X}_1}(\mathbf{t}_1) M_{\mathbf{X}_2}(\mathbf{t}_2)\), which is exactly the factorisation that characterises independence. This is a property of the normal and of almost nothing else: in general, uncorrelated variables need not be independent, as Theory of Probability, Unit 3 shows with a counterexample.
The density is constant exactly where the quadratic form is constant, so the contours are the ellipsoids
\[ (\mathbf{x} - \boldsymbol\mu)'\boldsymbol\Sigma^{-1}(\mathbf{x} - \boldsymbol\mu) = c^{2}, \]centred at \(\boldsymbol\mu\). Their axes point along the eigenvectors of \(\boldsymbol\Sigma\) and have half-lengths \(c\sqrt{\lambda_i}\) — a fact that returns as principal component analysis in Unit 4, where it is the whole idea. Because the quadratic form has the \(\chi^{2}_{p}\) distribution (Result (a) with \(\boldsymbol\Sigma^{-1/2}\), then Distribution Theory, Unit 4), choosing \(c^{2}\) as the upper \(5\%\) point of \(\chi^{2}_{p}\) makes the ellipsoid contain exactly \(95\%\) of the probability.
Partition the vector, the mean and the dispersion matrix conformably, with \(\mathbf{X}_1\) of length \(q\) and \(\mathbf{X}_2\) of length \(p - q\):
\[ \mathbf{X} = \begin{pmatrix}\mathbf{X}_1\\ \mathbf{X}_2\end{pmatrix}, \qquad \boldsymbol\mu = \begin{pmatrix}\boldsymbol\mu_1\\ \boldsymbol\mu_2\end{pmatrix}, \qquad \boldsymbol\Sigma = \begin{pmatrix} \boldsymbol\Sigma_{11} & \boldsymbol\Sigma_{12}\\ \boldsymbol\Sigma_{21} & \boldsymbol\Sigma_{22}\end{pmatrix}. \]Then the marginal is \(\mathbf{X}_1 \sim N_q(\boldsymbol\mu_1, \boldsymbol\Sigma_{11})\), and the conditional is again normal:
\[ \mathbf{X}_1 \mid \mathbf{X}_2 = \mathbf{x}_2 \;\sim\; N_q\!\left(\boldsymbol\mu_{1\cdot 2},\ \boldsymbol\Sigma_{11\cdot 2}\right), \] \[ \boldsymbol\mu_{1\cdot 2} = \boldsymbol\mu_1 + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} (\mathbf{x}_2 - \boldsymbol\mu_2), \qquad \boldsymbol\Sigma_{11\cdot 2} = \boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{21}. \]Statement. The conditional distribution above holds whenever \(\boldsymbol\Sigma_{22}\) is non-singular.
Step 1 — build a part that is uncorrelated with \(\mathbf{X}_2\). Define
\[ \mathbf{Z} = \mathbf{X}_1 - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\mathbf{X}_2. \]This is a linear transformation of \(\mathbf{X}\), so by Result (a) of section 3 the pair \((\mathbf{Z}, \mathbf{X}_2)\) is jointly normal.
Step 2 — compute the covariance between them. Using bilinearity of covariance,
\[ \operatorname{Cov}(\mathbf{Z}, \mathbf{X}_2) = \operatorname{Cov}(\mathbf{X}_1, \mathbf{X}_2) - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} \operatorname{Cov}(\mathbf{X}_2, \mathbf{X}_2) = \boldsymbol\Sigma_{12} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{22} = \mathbf{0}. \]The choice of coefficient matrix in Step 1 was made precisely to produce this zero.
Step 3 — upgrade zero covariance to independence. \(\mathbf{Z}\) and \(\mathbf{X}_2\) are jointly normal and uncorrelated, so by Result (c) of section 3 they are independent. Therefore conditioning on \(\mathbf{X}_2\) does not change the distribution of \(\mathbf{Z}\).
Step 4 — read off the conditional mean. Since \(\mathbf{X}_1 = \mathbf{Z} + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\mathbf{X}_2\) and the second term is a constant once \(\mathbf{X}_2 = \mathbf{x}_2\) is fixed,
\[ E(\mathbf{X}_1 \mid \mathbf{x}_2) = E(\mathbf{Z}) + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\mathbf{x}_2 = \left(\boldsymbol\mu_1 - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} \boldsymbol\mu_2\right) + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\mathbf{x}_2, \]which rearranges to \(\boldsymbol\mu_1 + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} (\mathbf{x}_2 - \boldsymbol\mu_2)\).
Step 5 — read off the conditional variance. Fixing \(\mathbf{x}_2\) adds a constant, which does not affect a variance, so the conditional dispersion matrix is \(\operatorname{Cov}(\mathbf{Z})\). Expanding by rule (i) of section 1 with \(\mathbf{A} = (\mathbf{I}, -\boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1})\),
\[ \operatorname{Cov}(\mathbf{Z}) = \boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21} + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{22} \boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21}, \]and the last two terms cancel, leaving \(\boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{21}\). \(\blacksquare\)
Two things worth noticing. The conditional mean is linear in \(\mathbf{x}_2\) — which is why linear regression is the right model when the variables are jointly normal, and only then. And the conditional variance does not depend on \(\mathbf{x}_2\) at all: the multivariate normal is homoscedastic by construction.
Given.
\[ \boldsymbol\mu = \begin{pmatrix}10\\20\\30\end{pmatrix}, \qquad \boldsymbol\Sigma = \begin{pmatrix} 9 & 6 & 3\\ 6 & 16 & 4\\ 3 & 4 & 25\end{pmatrix}. \]Asked. The correlations; the distribution of \(X_1\) given \(X_2 = 24\) and \(X_3 = 35\); and the multiple correlation of \(X_1\) on the other two.
Step 1 — read the standard deviations and correlations. The diagonal gives \(\sigma_1 = 3\), \(\sigma_2 = 4\), \(\sigma_3 = 5\), so
\[ \rho_{12} = \frac{6}{3 \times 4} = 0.500000, \qquad \rho_{13} = \frac{3}{3 \times 5} = 0.200000, \qquad \rho_{23} = \frac{4}{4 \times 5} = 0.200000. \]Step 2 — check that \(\boldsymbol\Sigma\) is a legal dispersion matrix. Its determinant is
\[ |\boldsymbol\Sigma| = 9(16 \times 25 - 16) - 6(6 \times 25 - 12) + 3(24 - 48) = 9(384) - 6(138) + 3(-24) = 3456 - 828 - 72 = 2556, \]and the leading minors are \(9 > 0\), \(9(16) - 36 = 108 > 0\) and \(2556 > 0\), so by the leading-minor test of Linear Algebra, Unit 3 it is positive definite. The density exists.
Step 3 — partition. Take \(\mathbf{X}_1 = X_1\) and \(\mathbf{X}_2 = (X_2, X_3)'\):
\[ \boldsymbol\Sigma_{11} = 9, \qquad \boldsymbol\Sigma_{12} = \begin{pmatrix} 6 & 3\end{pmatrix}, \qquad \boldsymbol\Sigma_{22} = \begin{pmatrix}16 & 4\\ 4 & 25\end{pmatrix}. \]Step 4 — invert the \(2 \times 2\) block. \(|\boldsymbol\Sigma_{22}| = 400 - 16 = 384\), so
\[ \boldsymbol\Sigma_{22}^{-1} = \frac{1}{384}\begin{pmatrix}25 & -4\\ -4 & 16\end{pmatrix}. \]Step 5 — the regression coefficients.
\[ \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1} = \frac{1}{384}\begin{pmatrix}6 & 3\end{pmatrix} \begin{pmatrix}25 & -4\\ -4 & 16\end{pmatrix} = \frac{1}{384}\begin{pmatrix}150 - 12 & -24 + 48\end{pmatrix} = \frac{1}{384}\begin{pmatrix}138 & 24\end{pmatrix}, \] \[ = \begin{pmatrix}\tfrac{23}{64} & \tfrac{1}{16}\end{pmatrix} = \begin{pmatrix}0.359375 & 0.062500\end{pmatrix}. \]Step 6 — the conditional mean.
\[ E(X_1 \mid X_2 = 24,\, X_3 = 35) = 10 + 0.359375(24 - 20) + 0.062500(35 - 30) \] \[ = 10 + 1.437500 + 0.312500 = 11.750000. \]Step 7 — the conditional variance.
\[ \sigma_{11\cdot 23} = 9 - \begin{pmatrix}\tfrac{23}{64} & \tfrac{1}{16}\end{pmatrix} \begin{pmatrix}6\\3\end{pmatrix} = 9 - \left(\frac{138}{64} + \frac{3}{16}\right) = 9 - 2.343750 = 6.656250, \]with conditional standard deviation \(\sqrt{6.656250} = 2.579971\).
Step 8 — the multiple correlation. The variance explained is \(9 - 6.656250 = 2.343750\), so
\[ R^{2}_{1\cdot 23} = \frac{2.343750}{9} = 0.260417, \qquad R_{1\cdot 23} = \sqrt{0.260417} = 0.5103. \]Interpretation. Knowing \(X_2\) and \(X_3\) cuts the standard deviation of \(X_1\) from 3 to 2.58 — useful but not dramatic, and the number that says so is \(R = 0.51\). Note that the conditional variance 6.656250 is the same whatever values \(X_2\) and \(X_3\) take; only the centre moves. The same partition formula, with \(\boldsymbol\Sigma\) replaced by the correlation matrix, gives the multiple and partial correlations of Unit 2 and the Foundation treatment in Statistical Methods, Unit 3 — the difference is only that a matrix inverse replaces a three-variable formula.
Let \(\mathbf{x}_1, \ldots, \mathbf{x}_n\) be a random sample from \(N_p(\boldsymbol\mu, \boldsymbol\Sigma)\). The likelihood is
\[ L(\boldsymbol\mu, \boldsymbol\Sigma) = (2\pi)^{-np/2}|\boldsymbol\Sigma|^{-n/2} \exp\!\left\{-\tfrac12\sum_{i=1}^{n}(\mathbf{x}_i - \boldsymbol\mu)' \boldsymbol\Sigma^{-1}(\mathbf{x}_i - \boldsymbol\mu)\right\}. \]Step 1 — split the exponent around \(\bar{\mathbf{x}}\). Write \(\mathbf{x}_i - \boldsymbol\mu = (\mathbf{x}_i - \bar{\mathbf{x}}) + (\bar{\mathbf{x}} - \boldsymbol\mu)\) and expand. The cross term carries the factor \(\sum_i(\mathbf{x}_i - \bar{\mathbf{x}}) = \mathbf{0}\) and so vanishes, leaving
\[ \sum_{i}(\mathbf{x}_i - \boldsymbol\mu)'\boldsymbol\Sigma^{-1}(\mathbf{x}_i - \boldsymbol\mu) = \operatorname{tr}\!\left(\boldsymbol\Sigma^{-1}\mathbf{A}\right) + n(\bar{\mathbf{x}} - \boldsymbol\mu)'\boldsymbol\Sigma^{-1} (\bar{\mathbf{x}} - \boldsymbol\mu), \]where \(\mathbf{A} = \sum_i (\mathbf{x}_i - \bar{\mathbf{x}}) (\mathbf{x}_i - \bar{\mathbf{x}})'\) is the sum of squares and cross products matrix. The trace appears because a scalar equals its own trace and the trace is invariant under cyclic permutation.
Step 2 — maximise over \(\boldsymbol\mu\). Only the second term involves \(\boldsymbol\mu\), it is non-negative because \(\boldsymbol\Sigma^{-1}\) is positive definite, and it is zero exactly at \(\boldsymbol\mu = \bar{\mathbf{x}}\). So
\[ \hat{\boldsymbol\mu} = \bar{\mathbf{x}}, \qquad \text{whatever } \boldsymbol\Sigma. \]Step 3 — maximise over \(\boldsymbol\Sigma\). With \(\boldsymbol\mu\) at its maximum the log-likelihood is \(-\tfrac{n}{2}\ln|\boldsymbol\Sigma| - \tfrac12\operatorname{tr} (\boldsymbol\Sigma^{-1}\mathbf{A})\), which is maximised at
\[ \hat{\boldsymbol\Sigma} = \frac{\mathbf{A}}{n}. \]Note the divisor: the MLE uses \(n\), and is therefore biased. The unbiased estimator is \(\mathbf{S} = \mathbf{A}/(n-1)\), for the same reason as in one dimension — one degree of freedom has gone into estimating \(\boldsymbol\mu\).
(i) \(\bar{\mathbf{X}} \sim N_p\!\left(\boldsymbol\mu, \tfrac{1}{n}\boldsymbol\Sigma\right)\). This is immediate from Result (a) of section 3: \(\bar{\mathbf{X}}\) is a linear function of the stacked sample, its mean is \(\boldsymbol\mu\), and its dispersion is \(\tfrac{1}{n^{2}}\sum_i\boldsymbol\Sigma = \boldsymbol\Sigma/n\) by independence.
(ii) \(\bar{\mathbf{X}}\) and \(\mathbf{A}\) are independent, and \(\mathbf{A}\) has the Wishart distribution \(W_p(n-1, \boldsymbol\Sigma)\) of Unit 2.
Result (ii) is the exact multivariate analogue of the independence of \(\bar X\) and \(S^{2}\) proved by the Helmert transformation in Distribution Theory, Unit 3, and it is proved the same way: apply an orthogonal matrix whose first row is \((1/\sqrt n, \ldots, 1/\sqrt n)\) to the sample, so that the first transformed vector is \(\sqrt n\,\bar{\mathbf{X}}\) and \(\mathbf{A}\) is built from the remaining \(n-1\), which are independent of it. Without this result there is no Hotelling's \(T^{2}\), exactly as without its univariate version there is no \(t\) test.
Given. Four observations on two variables:
\[ (2, 5), \quad (4, 9), \quad (6, 7), \quad (8, 11). \]Asked. \(\hat{\boldsymbol\mu}\), \(\hat{\boldsymbol\Sigma}\), \(\mathbf{S}\), both generalized variances, and the correlation.
Step 1 — the mean vector.
\[ \bar x_1 = \frac{2+4+6+8}{4} = 5, \qquad \bar x_2 = \frac{5+9+7+11}{4} = 8. \]Step 2 — the deviations.
\[ (-3, -3), \quad (-1, 1), \quad (1, -1), \quad (3, 3). \]Step 3 — the sums of squares and cross products.
\[ a_{11} = 9 + 1 + 1 + 9 = 20, \qquad a_{22} = 9 + 1 + 1 + 9 = 20, \] \[ a_{12} = (-3)(-3) + (-1)(1) + (1)(-1) + (3)(3) = 9 - 1 - 1 + 9 = 16, \] \[ \mathbf{A} = \begin{pmatrix}20 & 16\\ 16 & 20\end{pmatrix}. \]Step 4 — the two estimators.
\[ \hat{\boldsymbol\Sigma} = \frac{\mathbf{A}}{4} = \begin{pmatrix}5 & 4\\ 4 & 5\end{pmatrix}, \qquad \mathbf{S} = \frac{\mathbf{A}}{3} = \begin{pmatrix}6.666667 & 5.333333\\ 5.333333 & 6.666667\end{pmatrix}. \]Step 5 — the generalized variances. Work from the exact fractions, not the decimals:
\[ |\hat{\boldsymbol\Sigma}| = 25 - 16 = 9, \qquad |\mathbf{S}| = \frac{|\mathbf{A}|}{3^{2}} = \frac{400 - 256}{9} = \frac{144}{9} = 16. \]A warning worth taking. Computing \(|\mathbf{S}|\) from the six-decimal entries printed above gives
\[ (6.666667)^{2} - (5.333333)^{2} = 44.444449 - 28.444441 = 16.000008, \]not \(16\). The error is tiny here, but a determinant multiplies its entries together, so rounding error multiplies too — and in a \(5 \times 5\) matrix it compounds five times over. Always take a determinant from the exact values.
Step 6 — the correlation. This is unaffected by the divisor:
\[ r = \frac{16}{\sqrt{20 \times 20}} = \frac{16}{20} = 0.8. \]Interpretation. \(|\hat{\boldsymbol\Sigma}| = 9\) is small relative to \(\hat\sigma_{11}\hat\sigma_{22} = 25\), and the ratio \(9/25 = 0.36 = 1 - r^{2}\) is not a coincidence: for \(p = 2\), \(|\boldsymbol\Sigma| = \sigma_{11}\sigma_{22}(1 - \rho^{2})\). A generalized variance collapses towards zero exactly as the variables approach an exact linear relation, which is what makes it the right one-number summary of multivariate spread — and what makes it the subject of the next unit.