Multivariate Analysis
This page pairs well with Statistics, Probability and Linear Algebra. It condenses my notes for UNSW MATH5855 (Multivariate Analysis).
Multivariate analysis studies \(p\) variables measured on the same \(n\) units at once. Almost every classical method in it is one of two things: a likelihood-ratio test under the multivariate normal, or an eigenproblem for a symmetric matrix built from \(\mathbf{S}\). So the page opens with the matrix facts that make the eigenproblems tractable, builds the normal and its sampling theory, and then runs through the methods: inference on means, correlation structure, dimension reduction (PCA, CCA, factor analysis), linear models (MANOVA, structural equations), and finally classification, clustering and copulas, where normality is relaxed or dropped.
Notation throughout: vectors are bold columns, \(\mathcal{M}_{n,p}\) is the space of \(n\times p\) matrices, \(\mathbf{1}_n\) is the vector of ones, and \(\lambda_1\ge\dots\ge\lambda_p\) are eigenvalues in decreasing order.
Matrix toolkit
Symmetric matrices diagonalise in an orthonormal basis, and that one fact generates square roots, inverses, extremal problems and projections. Everything later leans on this section; Linear Algebra has the general theory.
Let \(A\in\mathcal{M}_{p,p}\) be symmetric. Its eigenvalues are real, and there is an orthogonal \(\Gamma=[\boldsymbol\gamma_1,\dots,\boldsymbol\gamma_p]\) (so \(\Gamma^\top\Gamma=\Gamma\Gamma^\top=I_p\)) with \[ A=\Gamma\Lambda\Gamma^\top=\sum_{i=1}^p\lambda_i\boldsymbol\gamma_i\boldsymbol\gamma_i^\top,\qquad \Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_p), \] where \(A\boldsymbol\gamma_i=\lambda_i\boldsymbol\gamma_i\). Consequently \(\operatorname{tr}A=\sum_i\lambda_i\) and \(|A|=\prod_i\lambda_i\).
If \(A\mathbf{x}=\lambda\mathbf{x}\) with \(\mathbf{x}\ne\mathbf{0}\) complex, then \(\bar{\mathbf{x}}^\top A\mathbf{x}=\lambda\,\bar{\mathbf{x}}^\top\mathbf{x}\), and the left side equals its own conjugate because \(A\) is real symmetric; since \(\bar{\mathbf{x}}^\top\mathbf{x}>0\), \(\lambda\) is real. Eigenvectors for distinct eigenvalues are orthogonal: \(\lambda_i\boldsymbol\gamma_j^\top\boldsymbol\gamma_i=\boldsymbol\gamma_j^\top A\boldsymbol\gamma_i=(A\boldsymbol\gamma_j)^\top\boldsymbol\gamma_i=\lambda_j\boldsymbol\gamma_j^\top\boldsymbol\gamma_i\). For the full decomposition, take one unit eigenvector \(\boldsymbol\gamma_1\), note that \(A\) maps \(\boldsymbol\gamma_1^\perp\) into itself (if \(\mathbf{x}\perp\boldsymbol\gamma_1\) then \(\boldsymbol\gamma_1^\top A\mathbf{x}=\lambda_1\boldsymbol\gamma_1^\top\mathbf{x}=0\)), and induct on the dimension. The trace and determinant follow from \(\operatorname{tr}(\Gamma\Lambda\Gamma^\top)=\operatorname{tr}(\Lambda\Gamma^\top\Gamma)\) and \(|\Gamma|^2=1\).
A symmetric \(A\) is positive definite if \(\mathbf{x}^\top A\mathbf{x}>0\) for every \(\mathbf{x}\ne\mathbf{0}\); writing \(\mathbf{x}=\Gamma\mathbf{y}\) gives \(\mathbf{x}^\top A\mathbf{x}=\sum_i\lambda_iy_i^2\), so this holds exactly when every \(\lambda_i>0\). For such \(A\) define \[ A^{-1}=\Gamma\Lambda^{-1}\Gamma^\top,\qquad A^{1/2}=\Gamma\Lambda^{1/2}\Gamma^\top,\qquad A^{-1/2}=\Gamma\Lambda^{-1/2}\Gamma^\top. \] \(A^{1/2}\) is symmetric positive definite, \(A^{1/2}A^{1/2}=A\), and \(A^{-1/2}AA^{-1/2}=I_p\). Every covariance matrix is positive semi-definite, since \(\mathbf{a}^\top\Sigma\mathbf{a}=\operatorname{Var}(\mathbf{a}^\top\mathbf{X})\ge0\).
A symmetric positive definite \(A\) factors uniquely as \(A=LL^\top\) with \(L\) lower triangular with positive diagonal. Software uses it to solve \(A\mathbf{x}=\mathbf{b}\), to compute \(|A|=\prod_iL_{ii}^2\), and to simulate: if \(\mathbf{Z}\sim N_p(\mathbf{0},I)\) then \(\boldsymbol\mu+L\mathbf{Z}\sim N_p(\boldsymbol\mu,A)\).
For symmetric \(A\) with spectral decomposition as above, \[ \max_{\mathbf{x}\ne\mathbf{0}}\frac{\mathbf{x}^\top A\mathbf{x}}{\mathbf{x}^\top\mathbf{x}}=\lambda_1,\qquad \min_{\mathbf{x}\ne\mathbf{0}}\frac{\mathbf{x}^\top A\mathbf{x}}{\mathbf{x}^\top\mathbf{x}}=\lambda_p, \] attained at \(\boldsymbol\gamma_1\) and \(\boldsymbol\gamma_p\). More generally the maximum over \(\mathbf{x}\perp\boldsymbol\gamma_1,\dots,\boldsymbol\gamma_{k-1}\) is \(\lambda_k\), attained at \(\boldsymbol\gamma_k\).
Put \(\mathbf{y}=\Gamma^\top\mathbf{x}\), so \(\mathbf{y}^\top\mathbf{y}=\mathbf{x}^\top\mathbf{x}\). The quotient becomes \(\sum_i\lambda_iy_i^2/\sum_iy_i^2\), a weighted average of the \(\lambda_i\) with weights \(y_i^2/\sum_jy_j^2\). A weighted average lies between the smallest and largest values and reaches them by putting all weight on one index. Orthogonality to \(\boldsymbol\gamma_1,\dots,\boldsymbol\gamma_{k-1}\) forces \(y_1=\dots=y_{k-1}=0\), leaving an average of \(\lambda_k,\dots,\lambda_p\).
\(P\in\mathcal{M}_{n,n}\) is an orthogonal projector onto a subspace \(\mathcal{V}\subseteq\mathbb{R}^n\) if \(P\mathbf{x}\in\mathcal{V}\) and \(\mathbf{x}-P\mathbf{x}\perp\mathcal{V}\) for every \(\mathbf{x}\).
\(P\) is an orthogonal projector if and only if it is symmetric and idempotent (\(P^\top=P=P^2\)). If the columns of \(X\in\mathcal{M}_{n,k}\) are a basis of \(\mathcal{V}\) then \[ P=X(X^\top X)^{-1}X^\top, \] \(I-P\) projects onto \(\mathcal{V}^\perp\), the eigenvalues of \(P\) are \(0\) and \(1\), and \(\operatorname{rank}P=\operatorname{tr}P=k\).
Write \(P\mathbf{x}=X\mathbf{b}\). The condition \(X^\top(\mathbf{x}-X\mathbf{b})=\mathbf{0}\) gives \(\mathbf{b}=(X^\top X)^{-1}X^\top\mathbf{x}\), which yields the formula; symmetry and idempotence are then direct. Conversely if \(P\) is symmetric idempotent, set \(\mathcal{V}=\operatorname{col}(P)\): \(P\mathbf{x}\in\mathcal{V}\), and for \(P\mathbf{y}\in\mathcal{V}\), \((P\mathbf{y})^\top(\mathbf{x}-P\mathbf{x})=\mathbf{y}^\top(P-P^2)\mathbf{x}=0\). From \(P\mathbf{v}=\lambda\mathbf{v}\) and \(P^2=P\) we get \(\lambda^2=\lambda\), and the trace is the number of unit eigenvalues.
\(H=I_n-\tfrac1n\mathbf{1}_n\mathbf{1}_n^\top\) is the projector onto \(\mathbf{1}_n^\perp\): it subtracts the mean from each coordinate. It has rank \(n-1\), which is where every “\(n-1\) degrees of freedom” in this page comes from.
Data and summaries
The data are an \(n\times p\) matrix whose rows are cases and columns are variables. Three summaries, all quadratic in the data, carry almost all the information the normal-theory methods use.
For a data matrix \(\mathbf{X}\in\mathcal{M}_{n,p}\) with rows \(\mathbf{x}_i^\top\), \[ \bar{\mathbf{x}}=\frac1n\sum_{i=1}^n\mathbf{x}_i=\frac1n\mathbf{X}^\top\mathbf{1}_n,\qquad \mathbf{S}=\frac1{n-1}\sum_{i=1}^n(\mathbf{x}_i-\bar{\mathbf{x}})(\mathbf{x}_i-\bar{\mathbf{x}})^\top=\frac1{n-1}\mathbf{X}^\top H\mathbf{X}, \] and \(\mathbf{R}=D^{-1/2}\mathbf{S}D^{-1/2}\) with \(D=\operatorname{diag}(s_{11},\dots,s_{pp})\), so \(r_{ik}=s_{ik}/\sqrt{s_{ii}s_{kk}}\). The sum-of-squares-and-products matrix is \(\mathbf{V}=\mathbf{X}^\top H\mathbf{X}=(n-1)\mathbf{S}\).
\(|r_{ik}|\le1\), with equality exactly when the two columns are affinely related, and \(r_{ik}\) is unchanged by changes of location and of positive scale in either variable. \(\mathbf{S}\) and \(\mathbf{R}\) are positive semi-definite, and positive definite when \(n>p\) and no column is an affine combination of the others.
\(r_{ik}\) is the cosine of the angle between the centred columns \(H\mathbf{x}_{(i)}\) and \(H\mathbf{x}_{(k)}\), so Cauchy-Schwarz gives the bound and its equality case. Scaling a column by \(c>0\) scales \(s_{ik}\) and \(\sqrt{s_{ii}}\) by \(c\), which cancels. For semi-definiteness, \(\mathbf{a}^\top\mathbf{V}\mathbf{a}=\lVert H\mathbf{X}\mathbf{a}\rVert^2\ge0\).
The generalised variance is \(|\mathbf{S}|\) and the total variance is \(\operatorname{tr}\mathbf{S}\). The Mahalanobis distance of \(\mathbf{x}\) from \(\boldsymbol\mu\) under covariance \(\Sigma\) is \[ d_\Sigma(\mathbf{x},\boldsymbol\mu)=\sqrt{(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)}. \] It is the Euclidean distance after the whitening map \(\mathbf{x}\mapsto\Sigma^{-1/2}\mathbf{x}\), so it discounts directions of large variance and accounts for correlation.
Before any model: look at the scatterplot matrix, the marginal histograms and \(\mathbf{R}\). Skewness should prompt a transformation (the microwave-oven example below takes fourth roots), and outliers show up as large Mahalanobis distances. Generalised variance hides structure: very different \(\mathbf{S}\) can share a determinant.
The multivariate normal
The multivariate normal is defined through its linear combinations, and the definition makes its closure properties almost free: linear maps, marginals and conditionals all stay normal, and zero correlation means independence. See Probability for characteristic functions in one dimension.
A random vector \(\mathbf{X}\in\mathbb{R}^p\) is multivariate normal, \(\mathbf{X}\sim N_p(\boldsymbol\mu,\Sigma)\), if \(\mathbf{a}^\top\mathbf{X}\) is univariate normal (possibly degenerate) for every \(\mathbf{a}\in\mathbb{R}^p\). Equivalently its characteristic function is \[ \varphi_{\mathbf{X}}(\mathbf{t})=\mathbb{E}\,e^{i\mathbf{t}^\top\mathbf{X}}=\exp\!\Big(i\mathbf{t}^\top\boldsymbol\mu-\tfrac12\mathbf{t}^\top\Sigma\mathbf{t}\Big). \] When \(\Sigma\) is positive definite it has density \[ f(\mathbf{x})=(2\pi)^{-p/2}|\Sigma|^{-1/2}\exp\!\Big(-\tfrac12(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)\Big), \] whose level sets are the ellipsoids of constant Mahalanobis distance, with axes along the eigenvectors of \(\Sigma\) and half-lengths proportional to \(\sqrt{\lambda_i}\).
(Of the equivalence.) If every \(\mathbf{t}^\top\mathbf{X}\) is normal, it has mean \(\mathbf{t}^\top\boldsymbol\mu\) and variance \(\mathbf{t}^\top\Sigma\mathbf{t}\), so \(\varphi_{\mathbf{X}}(\mathbf{t})=\varphi_{\mathbf{t}^\top\mathbf{X}}(1)\) is the displayed function. Conversely \(\varphi_{\mathbf{a}^\top\mathbf{X}}(s)=\varphi_{\mathbf{X}}(s\mathbf{a})\) is the characteristic function of \(N(\mathbf{a}^\top\boldsymbol\mu,\mathbf{a}^\top\Sigma\mathbf{a})\), and characteristic functions determine laws. For the density, \(\mathbf{X}=\boldsymbol\mu+\Sigma^{1/2}\mathbf{Z}\) with \(\mathbf{Z}\) a vector of independent \(N(0,1)\) variables has the right characteristic function; change variables with Jacobian \(|\Sigma|^{-1/2}\).
If \(\mathbf{X}\sim N_p(\boldsymbol\mu,\Sigma)\), \(A\in\mathcal{M}_{q,p}\) and \(\mathbf{b}\in\mathbb{R}^q\), then \(A\mathbf{X}+\mathbf{b}\sim N_q(A\boldsymbol\mu+\mathbf{b},A\Sigma A^\top)\). In particular every sub-vector is normal with the corresponding sub-blocks of \(\boldsymbol\mu\) and \(\Sigma\).
\(\varphi_{A\mathbf{X}+\mathbf{b}}(\mathbf{t})=e^{i\mathbf{t}^\top\mathbf{b}}\varphi_{\mathbf{X}}(A^\top\mathbf{t})=\exp\big(i\mathbf{t}^\top(A\boldsymbol\mu+\mathbf{b})-\tfrac12\mathbf{t}^\top A\Sigma A^\top\mathbf{t}\big)\). Marginals take \(A\) to be a selection matrix.
With \(\Sigma_{12}=0\) the quadratic form splits, \(\mathbf{t}^\top\Sigma\mathbf{t}=\mathbf{t}_1^\top\Sigma_{11}\mathbf{t}_1+\mathbf{t}_2^\top\Sigma_{22}\mathbf{t}_2\), so the joint characteristic function is the product of the marginal ones, which characterises independence. The converse holds for any distribution with finite second moments.
With the partition above and \(\Sigma_{22}\) invertible, \[ \mathbf{X}_1\mid\mathbf{X}_2=\mathbf{x}_2\ \sim\ N\Big(\boldsymbol\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf{x}_2-\boldsymbol\mu_2),\ \Sigma_{11\cdot2}\Big),\qquad \Sigma_{11\cdot2}=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}. \] The conditional mean is linear in \(\mathbf{x}_2\) and the conditional covariance does not depend on \(\mathbf{x}_2\) at all.
Let \(\mathbf{Z}=\mathbf{X}_1-\Sigma_{12}\Sigma_{22}^{-1}\mathbf{X}_2\). The pair \((\mathbf{Z},\mathbf{X}_2)\) is a linear map of \(\mathbf{X}\), hence jointly normal, and \[ \operatorname{Cov}(\mathbf{Z},\mathbf{X}_2)=\Sigma_{12}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{22}=0, \] so \(\mathbf{Z}\) is independent of \(\mathbf{X}_2\). It has mean \(\boldsymbol\mu_1-\Sigma_{12}\Sigma_{22}^{-1}\boldsymbol\mu_2\) and covariance \(\Sigma_{11}-2\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}+\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{22}\Sigma_{22}^{-1}\Sigma_{21}=\Sigma_{11\cdot2}\). Given \(\mathbf{X}_2=\mathbf{x}_2\), \(\mathbf{X}_1=\mathbf{Z}+\Sigma_{12}\Sigma_{22}^{-1}\mathbf{x}_2\) with \(\mathbf{Z}\) unaffected by the conditioning.
For \((X_1,X_2)\) bivariate normal with standard deviations \(\sigma_1,\sigma_2\) and correlation \(\rho\), \[ X_1\mid X_2=x_2\sim N\Big(\mu_1+\rho\frac{\sigma_1}{\sigma_2}(x_2-\mu_2),\ \sigma_1^2(1-\rho^2)\Big). \] The slope \(\rho\sigma_1/\sigma_2\) is the population regression coefficient, and knowing \(X_2\) removes the fraction \(\rho^2\) of the variance of \(X_1\).
For any \(Y\) with finite variance and any vector \(\mathbf{X}\), the function \(g\) minimising \(\mathbb{E}\{Y-g(\mathbf{X})\}^2\) is \(g^*(\mathbf{X})=\mathbb{E}(Y\mid\mathbf{X})\). If \((Y,\mathbf{X})\) is jointly normal with \(\operatorname{Var}\mathbf{X}=C\) and \(\operatorname{Cov}(\mathbf{X},Y)=\boldsymbol\sigma_0\), then \(g^*\) is linear: \[ g^*(\mathbf{x})=\mathbb{E}Y+\boldsymbol\sigma_0^\top C^{-1}(\mathbf{x}-\mathbb{E}\mathbf{X}). \]
\(\mathbb{E}\{Y-g\}^2=\mathbb{E}\{Y-g^*\}^2+\mathbb{E}\{g^*-g\}^2\), because the cross term \(\mathbb{E}[(Y-g^*)(g^*-g)]\) vanishes on conditioning on \(\mathbf{X}\). The normal case is the conditional mean of the theorem above with \(\mathbf{X}_1=Y\).
If \(\mathbf{X}\sim N_p(\boldsymbol\mu,\Sigma)\) with \(\Sigma\) positive definite then \((\mathbf{X}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{X}-\boldsymbol\mu)\sim\chi^2_p\). Hence the ellipsoid \(\{\mathbf{x}:(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)\le\chi^2_{p}(1-\alpha)\}\) carries probability \(1-\alpha\).
\(\mathbf{Z}=\Sigma^{-1/2}(\mathbf{X}-\boldsymbol\mu)\sim N_p(\mathbf{0},I_p)\) has independent standard normal components, and the quadratic form is \(\mathbf{Z}^\top\mathbf{Z}=\sum_jZ_j^2\).
Let \(g_{ij}=(\mathbf{x}_i-\bar{\mathbf{x}})^\top\mathbf{S}_n^{-1}(\mathbf{x}_j-\bar{\mathbf{x}})\) with \(\mathbf{S}_n=\mathbf{V}/n\). Mardia’s sample skewness and kurtosis are \[ \hat\beta_{1,p}=\frac1{n^2}\sum_{i,j}g_{ij}^3,\qquad \hat\beta_{2,p}=\frac1n\sum_ig_{ii}^2. \] Under multivariate normality \(\beta_{1,p}=0\) and \(\beta_{2,p}=p(p+2)\), and asymptotically \[ \frac{n\hat\beta_{1,p}}6\sim\chi^2_{p(p+1)(p+2)/6},\qquad \frac{\hat\beta_{2,p}-p(p+2)}{\sqrt{8p(p+2)/n}}\sim N(0,1). \]
Marginal Q-Q plots and univariate tests check only the margins; a vector can have normal margins without being jointly normal (the remark above is an example, and the copulas section manufactures many more). Mardia’s statistics, and a Q-Q plot of the \(g_{ii}\) against \(\chi^2_p\) quantiles, look at the joint law.
Estimation and the Wishart distribution
For a normal sample, \(\bar{\mathbf{X}}\) and \(\mathbf{V}\) are sufficient, maximum likelihood, and independent. The law of \(\mathbf{V}\) is the matrix generalisation of \(\sigma^2\chi^2\), and the scalar pivots extracted from it drive every test later on.
For \(\mathbf{x}_1,\dots,\mathbf{x}_n\) independent \(N_p(\boldsymbol\mu,\Sigma)\), \[ \ell(\boldsymbol\mu,\Sigma)=-\frac{np}2\log(2\pi)-\frac n2\log|\Sigma|-\frac12\operatorname{tr}\Big(\Sigma^{-1}\big[\mathbf{V}+n(\bar{\mathbf{x}}-\boldsymbol\mu)(\bar{\mathbf{x}}-\boldsymbol\mu)^\top\big]\Big). \] So \((\bar{\mathbf{x}},\mathbf{V})\) is sufficient for \((\boldsymbol\mu,\Sigma)\).
Each quadratic form is a scalar, so equals its trace, and \(\operatorname{tr}(\mathbf{a}^\top B\mathbf{a})=\operatorname{tr}(B\mathbf{a}\mathbf{a}^\top)\). Summing and adding and subtracting \(\bar{\mathbf{x}}\) inside each term gives \(\sum_i(\mathbf{x}_i-\boldsymbol\mu)(\mathbf{x}_i-\boldsymbol\mu)^\top=\mathbf{V}+n(\bar{\mathbf{x}}-\boldsymbol\mu)(\bar{\mathbf{x}}-\boldsymbol\mu)^\top\); the cross terms vanish because \(\sum_i(\mathbf{x}_i-\bar{\mathbf{x}})=\mathbf{0}\).
For symmetric positive definite \(A\in\mathcal{M}_{p,p}\), the function \(h(G)=-n\log|G|-\operatorname{tr}(G^{-1}A)\) on positive definite \(G\) is maximised at \(G=A/n\), with maximum \(np\log n-n\log|A|-np\).
Let \(\eta_1,\dots,\eta_p>0\) be the eigenvalues of \(A^{1/2}G^{-1}A^{1/2}\). Then \(\operatorname{tr}(G^{-1}A)=\sum_j\eta_j\) and \(|G|^{-1}=|A|^{-1}\prod_j\eta_j\), so \[ h(G)=-n\log|A|+\sum_{j=1}^p\big(n\log\eta_j-\eta_j\big). \] Each summand is maximised at \(\eta_j=n\). All eigenvalues equal to \(n\) means \(A^{1/2}G^{-1}A^{1/2}=nI\), that is \(G=A/n\).
If \(n>p\), the MLEs are \(\hat{\boldsymbol\mu}=\bar{\mathbf{x}}\) and \(\hat\Sigma=\mathbf{V}/n=\frac{n-1}n\mathbf{S}\). By invariance the MLE of the correlation matrix is \(\mathbf{R}\) (the factor \((n-1)/n\) cancels).
For any \(\Sigma\) the trace term involving \(\boldsymbol\mu\) is a non-negative Mahalanobis form, minimised to zero by \(\boldsymbol\mu=\bar{\mathbf{x}}\). What is left is \(\tfrac12h(\Sigma)\) with \(A=\mathbf{V}\), positive definite with probability one when \(n>p\); apply Anderson’s lemma.
If \(\mathbf{Z}_1,\dots,\mathbf{Z}_m\) are independent \(N_p(\mathbf{0},\Sigma)\), the law of \(W=\sum_{j=1}^m\mathbf{Z}_j\mathbf{Z}_j^\top\) is the Wishart distribution \(W_p(\Sigma,m)\) with \(m\) degrees of freedom. For \(p=1\) it is \(\sigma^2\chi^2_m\).
For a sample of size \(n\) from \(N_p(\boldsymbol\mu,\Sigma)\):
- \(\bar{\mathbf{X}}\sim N_p(\boldsymbol\mu,\Sigma/n)\);
- \(\mathbf{V}=n\hat\Sigma=(n-1)\mathbf{S}\sim W_p(\Sigma,n-1)\);
- \(\bar{\mathbf{X}}\) and \(\mathbf{V}\) are independent.
Put the observations as the columns of \(\mathbf{Y}=\mathbf{X}^\top\in\mathcal{M}_{p,n}\). Choose an orthogonal \(\Gamma\in\mathcal{M}_{n,n}\) whose first column is \(\mathbf{1}_n/\sqrt n\) (Gram-Schmidt from any basis starting with \(\mathbf{1}_n\)) and rotate: \(\mathbf{U}=\mathbf{Y}\Gamma\). Because the columns of \(\mathbf{Y}\) are independent with common covariance \(\Sigma\), the covariance of the stacked vector is \(I_n\otimes\Sigma\), and the rotation sends it to \((\Gamma^\top\Gamma)\otimes\Sigma=I_n\otimes\Sigma\): the columns \(\mathbf{u}_1,\dots,\mathbf{u}_n\) of \(\mathbf{U}\) are again independent \(N_p(\cdot,\Sigma)\). Their means are \(\mathbb{E}\mathbf{u}_1=\sqrt n\boldsymbol\mu\) and \(\mathbb{E}\mathbf{u}_j=\boldsymbol\mu\mathbf{1}_n^\top\boldsymbol\gamma_j=\mathbf{0}\) for \(j\ge2\) by orthogonality to \(\mathbf{1}_n\). Now \(\mathbf{u}_1=\sqrt n\,\bar{\mathbf{X}}\), and \(\mathbf{Y}\mathbf{Y}^\top=\mathbf{U}\mathbf{U}^\top\) gives \(\mathbf{V}=\mathbf{Y}\mathbf{Y}^\top-n\bar{\mathbf{X}}\bar{\mathbf{X}}^\top=\sum_{j=2}^n\mathbf{u}_j\mathbf{u}_j^\top\). That is a function of \(\mathbf{u}_2,\dots,\mathbf{u}_n\) only, hence independent of \(\mathbf{u}_1\), and is \(W_p(\Sigma,n-1)\) by definition.
Let \(W\sim W_p(\Sigma,m)\).
- For \(B\in\mathcal{M}_{k,p}\), \(BWB^\top\sim W_k(B\Sigma B^\top,m)\).
- For a fixed \(\mathbf{c}\ne\mathbf{0}\), \(\mathbf{c}^\top W\mathbf{c}/\mathbf{c}^\top\Sigma\mathbf{c}\sim\chi^2_m\).
- If \(m\ge p\), \(\mathbf{c}^\top\Sigma^{-1}\mathbf{c}/\mathbf{c}^\top W^{-1}\mathbf{c}\sim\chi^2_{m-p+1}\).
- Partitioning with an \(r\times r\) leading block, \(W_{11}\sim W_r(\Sigma_{11},m)\) and the Schur complement \(W_{11\cdot2}=W_{11}-W_{12}W_{22}^{-1}W_{21}\sim W_r(\Sigma_{11\cdot2},m-p+r)\).
(1) \(BWB^\top=\sum_j(B\mathbf{Z}_j)(B\mathbf{Z}_j)^\top\) with \(B\mathbf{Z}_j\sim N_k(\mathbf{0},B\Sigma B^\top)\). (2) is (1) with \(B=\mathbf{c}^\top\), then \(p=1\). (3) from (4): for \(\mathbf{c}=\mathbf{e}_1\) the partitioned inverse gives \((W^{-1})_{11}=1/W_{11\cdot2}\) and \((\Sigma^{-1})_{11}=1/\sigma_{11\cdot2}\), so the ratio is \(W_{11\cdot2}/\sigma_{11\cdot2}\sim\chi^2_{m-p+1}\) by (4) with \(r=1\) and (2). A general \(\mathbf{c}\) is moved to \(\mathbf{e}_1\) by a non-singular change of variables, which leaves both sides unchanged by (1). (4) is the conditional-distribution theorem applied column by column to the regression of the first block on the second, losing one degree of freedom per regressor; not proved here.
Inference on a mean vector
Hotelling’s \(T^2\) is the squared \(t\)-statistic maximised over all linear combinations of the variables. It is at once the union-intersection test, the likelihood-ratio test and the pivot behind simultaneous confidence intervals.
To test \(H_0:\boldsymbol\mu=\boldsymbol\mu_0\) from a sample of size \(n>p\), \[ T^2=n(\bar{\mathbf{x}}-\boldsymbol\mu_0)^\top\mathbf{S}^{-1}(\bar{\mathbf{x}}-\boldsymbol\mu_0). \] It is the sample Mahalanobis distance of \(\boldsymbol\mu_0\) from \(\bar{\mathbf{x}}\) in the metric of \(\mathbf{S}/n\), and it is invariant under \(\mathbf{x}\mapsto A\mathbf{x}+\mathbf{b}\) for non-singular \(A\).
Under \(H_0\), \[ \frac{n-p}{(n-1)p}\,T^2\sim F_{p,\,n-p}, \] and the level-\(\alpha\) test rejects when \(T^2>\frac{(n-1)p}{n-p}F_{p,n-p}(1-\alpha)\). As \(n\to\infty\), \(T^2\to\chi^2_p\) in distribution.
Write \(\mathbf{d}=\sqrt n(\bar{\mathbf{X}}-\boldsymbol\mu_0)\sim N_p(\mathbf{0},\Sigma)\), independent of \(\mathbf{V}=(n-1)\mathbf{S}\sim W_p(\Sigma,n-1)\). Then \[ T^2=\mathbf{d}^\top\mathbf{S}^{-1}\mathbf{d}=(n-1)\,\mathbf{d}^\top\mathbf{V}^{-1}\mathbf{d}=(n-1)\,\frac{\mathbf{d}^\top\Sigma^{-1}\mathbf{d}}{\mathbf{d}^\top\Sigma^{-1}\mathbf{d}/\mathbf{d}^\top\mathbf{V}^{-1}\mathbf{d}}. \] The numerator is \(\chi^2_p\). Conditionally on \(\mathbf{d}\), the denominator is \(\chi^2_{n-p}\) by Wishart property 3 with \(\mathbf{c}=\mathbf{d}\); that law does not depend on \(\mathbf{d}\), so the denominator is \(\chi^2_{n-p}\) and independent of \(\mathbf{d}\). Hence \(T^2/(n-1)=\chi^2_p/\chi^2_{n-p}\) with independent pieces, which is \(\frac p{n-p}F_{p,n-p}\).
For every \(\mathbf{a}\ne\mathbf{0}\) let \(t_{\mathbf{a}}^2=n(\mathbf{a}^\top\bar{\mathbf{x}}-\mathbf{a}^\top\boldsymbol\mu_0)^2/\mathbf{a}^\top\mathbf{S}\mathbf{a}\) be the squared \(t\)-statistic for the scalar hypothesis \(\mathbf{a}^\top\boldsymbol\mu=\mathbf{a}^\top\boldsymbol\mu_0\). Then \[ \max_{\mathbf{a}\ne\mathbf{0}}t_{\mathbf{a}}^2=T^2,\qquad\text{attained at }\mathbf{a}\propto\mathbf{S}^{-1}(\bar{\mathbf{x}}-\boldsymbol\mu_0). \]
Put \(\mathbf{b}=\mathbf{S}^{1/2}\mathbf{a}\) and \(\mathbf{d}=\mathbf{S}^{-1/2}(\bar{\mathbf{x}}-\boldsymbol\mu_0)\). Then \(t_{\mathbf{a}}^2=n(\mathbf{b}^\top\mathbf{d})^2/\mathbf{b}^\top\mathbf{b}\le n\,\mathbf{d}^\top\mathbf{d}=T^2\) by Cauchy-Schwarz, with equality when \(\mathbf{b}\propto\mathbf{d}\).
With \(\hat\Sigma=\mathbf{V}/n\) and \(\hat\Sigma_0=\frac1n\sum_i(\mathbf{x}_i-\boldsymbol\mu_0)(\mathbf{x}_i-\boldsymbol\mu_0)^\top\), the likelihood ratio is \(\Lambda=(|\hat\Sigma|/|\hat\Sigma_0|)^{n/2}\), and Wilks’ lambda satisfies \[ \Lambda^{2/n}=\frac{|\hat\Sigma|}{|\hat\Sigma_0|}=\Big(1+\frac{T^2}{n-1}\Big)^{-1}. \] So the LRT rejects for large \(T^2\): the two tests are the same test.
The ratio of maxima follows from Anderson’s lemma applied with \(A=n\hat\Sigma_0\) and \(A=n\hat\Sigma\). For the identity, \(n\hat\Sigma_0=\mathbf{V}+\mathbf{u}\mathbf{u}^\top\) with \(\mathbf{u}=\sqrt n(\bar{\mathbf{x}}-\boldsymbol\mu_0)\), and the matrix determinant lemma gives \(|\mathbf{V}+\mathbf{u}\mathbf{u}^\top|=|\mathbf{V}|(1+\mathbf{u}^\top\mathbf{V}^{-1}\mathbf{u})=|\mathbf{V}|(1+T^2/(n-1))\).
A \(100(1-\alpha)\%\) confidence region for \(\boldsymbol\mu\) is the ellipsoid \[ \Big\{\boldsymbol\mu:\ n(\bar{\mathbf{x}}-\boldsymbol\mu)^\top\mathbf{S}^{-1}(\bar{\mathbf{x}}-\boldsymbol\mu)\le c^2\Big\},\qquad c^2=\frac{(n-1)p}{n-p}F_{p,n-p}(1-\alpha), \] centred at \(\bar{\mathbf{x}}\) with axes along the eigenvectors of \(\mathbf{S}\). Simultaneously for all \(\mathbf{a}\in\mathbb{R}^p\), with probability at least \(1-\alpha\), \[ \mathbf{a}^\top\boldsymbol\mu\in\Big[\mathbf{a}^\top\bar{\mathbf{x}}\pm c\sqrt{\mathbf{a}^\top\mathbf{S}\mathbf{a}/n}\Big]. \] For a handful \(m\) of pre-chosen combinations, Bonferroni intervals \(\mathbf{a}_k^\top\bar{\mathbf{x}}\pm t_{n-1}(1-\tfrac{\alpha}{2m})\sqrt{\mathbf{a}_k^\top\mathbf{S}\mathbf{a}_k/n}\) are usually shorter.
The region is the inversion of the \(T^2\) test. The intervals: \(\boldsymbol\mu\) is in the ellipsoid if and only if \(\max_{\mathbf{a}}t_{\mathbf{a}}^2\le c^2\) (union of \(t\)-tests, with \(\boldsymbol\mu\) in place of \(\boldsymbol\mu_0\)), that is if and only if every interval covers. Bonferroni follows from \(P(\bigcup_kE_k)\le\sum_kP(E_k)\).
Radiation was measured on \(n=42\) ovens with the door closed and open. The raw data are right-skewed, so take fourth roots; the transformed sample has \(\bar{\mathbf{x}}=(0.5643,\,0.6030)^\top\) and \[ \mathbf{S}=\begin{pmatrix}0.014350&0.011715\\0.011715&0.014545\end{pmatrix}. \] To test \(H_0:\boldsymbol\mu=(0.6,0.6)^\top\), compute \(T^2=12.46\). The 5% critical value is \(\frac{41\cdot2}{40}F_{2,40}(0.95)=2.05\times3.23=6.62\), so \(H_0\) is rejected. The determinant route gives the same number without inverting \(\mathbf{S}\): \(|\hat\Sigma_0|/|\hat\Sigma|=1.3039\), so \(T^2=41(1.3039-1)=12.46\) and Wilks’ lambda is \(0.767\). All conclusions are about the mean of the fourth root of radiation, not of radiation.
The same pivot covers paired comparisons (apply \(T^2\) to differences), repeated measures (test \(C\boldsymbol\mu=\mathbf{0}\) for a contrast matrix \(C\) with \(p-1\) rows, giving \(F_{p-1,n-p+1}\)), and two independent samples with a common covariance (replace \(\mathbf{S}\) by the pooled \(\mathbf{S}_{\mathrm{p}}\) and \(n\) by \(n_1n_2/(n_1+n_2)\), with \(n_1+n_2-2\) error degrees of freedom).
Partial and multiple correlation
Correlation between two variables can be produced entirely by a third. Partial correlation removes the third by conditioning; multiple correlation asks how well one variable can be predicted linearly from many. Both are read off \(\Sigma\) or its inverse.
Partition \(\mathbf{X}=(\mathbf{X}_1^\top,\mathbf{X}_2^\top)^\top\) and let \(\Sigma_{11\cdot2}=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\), the conditional covariance of \(\mathbf{X}_1\) given \(\mathbf{X}_2\) under normality. The partial correlation of \(X_i\) and \(X_j\) (both in \(\mathbf{X}_1\)) given \(\mathbf{X}_2\) is \[ \rho_{ij\cdot2}=\frac{(\Sigma_{11\cdot2})_{ij}}{\sqrt{(\Sigma_{11\cdot2})_{ii}(\Sigma_{11\cdot2})_{jj}}}. \] Its MLE replaces \(\Sigma\) by \(\hat\Sigma\).
With a single conditioning variable, \[ \rho_{12\cdot3}=\frac{\rho_{12}-\rho_{13}\rho_{23}}{\sqrt{(1-\rho_{13}^2)(1-\rho_{23}^2)}}. \] More generally, if \(P=\rho^{-1}\) is the inverse correlation matrix, the partial correlation of \(X_i\) and \(X_j\) given all the other variables is \(-P_{ij}/\sqrt{P_{ii}P_{jj}}\). A zero in \(\Sigma^{-1}\) therefore means conditional independence of that pair under normality: the starting point of Gaussian graphical models.
For schoolchildren with \(X_1\) intelligence, \(X_2\) weight, \(X_3\) age, take \(\rho_{12}=0.6162\), \(\rho_{13}=0.8267\), \(\rho_{23}=0.7321\). Then \[ \rho_{12\cdot3}=\frac{0.6162-0.8267\times0.7321}{\sqrt{(1-0.8267^2)(1-0.7321^2)}}=\frac{0.0110}{0.3833}=0.029. \] The marginal correlation of \(0.62\) was produced by both variables growing with age; holding age fixed leaves essentially none.
Let \((Y,\mathbf{X}^\top)\) have covariance \(\begin{pmatrix}\sigma_Y^2&\boldsymbol\sigma_0^\top\\boldsymbol\sigma_0&C\end{pmatrix}\). The multiple correlation coefficient is \[ R=\max_{\boldsymbol\alpha}\operatorname{cor}(Y,\boldsymbol\alpha^\top\mathbf{X})=\sqrt{\frac{\boldsymbol\sigma_0^\top C^{-1}\boldsymbol\sigma_0}{\sigma_Y^2}}, \] and \(R^2\) is the population coefficient of determination.
The maximising \(\boldsymbol\alpha\) is proportional to \(C^{-1}\boldsymbol\sigma_0\), the vector of best linear predictor coefficients, and \[ 1-R^2=\frac{\sigma_Y^2-\boldsymbol\sigma_0^\top C^{-1}\boldsymbol\sigma_0}{\sigma_Y^2}=\frac{\operatorname{Var}(Y\mid\mathbf{X})}{\operatorname{Var}Y}=\frac1{\rho^{YY}}, \] where \(\rho^{YY}\) is the leading entry of the inverse correlation matrix of \((Y,\mathbf{X})\).
\(\operatorname{cor}(Y,\boldsymbol\alpha^\top\mathbf{X})^2=(\boldsymbol\alpha^\top\boldsymbol\sigma_0)^2/(\sigma_Y^2\,\boldsymbol\alpha^\top C\boldsymbol\alpha)\). Cauchy-Schwarz in the inner product \(\langle\mathbf{u},\mathbf{v}\rangle=\mathbf{u}^\top C\mathbf{v}\), applied to \(\boldsymbol\alpha\) and \(C^{-1}\boldsymbol\sigma_0\), gives \((\boldsymbol\alpha^\top\boldsymbol\sigma_0)^2\le(\boldsymbol\alpha^\top C\boldsymbol\alpha)(\boldsymbol\sigma_0^\top C^{-1}\boldsymbol\sigma_0)\), with equality at \(\boldsymbol\alpha\propto C^{-1}\boldsymbol\sigma_0\). The last identity is the partitioned-inverse formula: the \((1,1)\) entry of the inverse is the reciprocal of the Schur complement, here computed on the correlation scale.
Principal component analysis
PCA rotates the coordinate system so the first axis carries as much variance as possible, the second as much of the rest, and so on. It is the spectral decomposition of \(\Sigma\) read statistically. The machine-learning view is on the PCA page.
Let \(\operatorname{Var}\mathbf{X}=\Sigma=\Gamma\Lambda\Gamma^\top\). The \(k\)-th principal component is \(Y_k=\boldsymbol\gamma_k^\top(\mathbf{X}-\boldsymbol\mu)\), the coefficients \(\boldsymbol\gamma_k\) are the loadings, and \(\mathbf{Y}=\Gamma^\top(\mathbf{X}-\boldsymbol\mu)\) is the vector of all components.
\(Y_1\) maximises \(\operatorname{Var}(\mathbf{a}^\top\mathbf{X})\) over unit vectors \(\mathbf{a}\); \(Y_k\) maximises it over unit vectors with \(\operatorname{Cov}(\mathbf{a}^\top\mathbf{X},Y_j)=0\) for \(j<k\). Moreover \(\operatorname{Var}Y_k=\lambda_k\), the components are uncorrelated, and \[ \sum_{k=1}^p\operatorname{Var}Y_k=\operatorname{tr}\Sigma=\sum_{k=1}^p\operatorname{Var}X_k, \] so the \(k\)-th component explains the proportion \(\lambda_k/\sum_j\lambda_j\) of the total variance.
\(\operatorname{Var}(\mathbf{a}^\top\mathbf{X})=\mathbf{a}^\top\Sigma\mathbf{a}\), and \(\operatorname{Cov}(\mathbf{a}^\top\mathbf{X},Y_j)=\mathbf{a}^\top\Sigma\boldsymbol\gamma_j=\lambda_j\mathbf{a}^\top\boldsymbol\gamma_j\), so (for \(\lambda_j>0\)) the constraint is \(\mathbf{a}\perp\boldsymbol\gamma_j\). Both statements are then the Rayleigh-quotient theorem. \(\operatorname{Var}\mathbf{Y}=\Gamma^\top\Sigma\Gamma=\Lambda\) gives the rest, and the trace is invariant under the rotation.
Equivalently, projecting onto the span of \(\boldsymbol\gamma_1,\dots,\boldsymbol\gamma_k\) minimises the expected squared reconstruction error \(\mathbb{E}\lVert(\mathbf{X}-\boldsymbol\mu)-P(\mathbf{X}-\boldsymbol\mu)\rVert^2\) over all rank-\(k\) orthogonal projectors \(P\); the minimum is \(\lambda_{k+1}+\dots+\lambda_p\). On data this is the truncated singular value decomposition of the centred data matrix.
Sample principal components use the eigen-decomposition of \(\mathbf{S}\) (or \(\hat\Sigma\); the eigenvectors are the same). PCA is not scale-invariant: a variable measured in millimetres rather than metres dominates the first component. When units are not comparable, standardise first, which amounts to PCA of the correlation matrix \(\mathbf{R}\); then the total variance is \(p\).
Common rules: keep enough components to explain a target proportion (say 80 to 90%) of total variance; with \(\mathbf{R}\), keep those with \(\lambda_k>1\) (Kaiser: better than an average single variable); or look for the elbow in the scree plot of \(\lambda_k\) against \(k\). Under normality one can test whether the trailing eigenvalues are equal: with \(A_0\) and \(G_0\) the arithmetic and geometric means of the last \(p-k\) sample eigenvalues of \(\mathbf{S}\), \[ n(p-k)\log\frac{A_0}{G_0}\ \dot\sim\ \chi^2_{(p-k+2)(p-k-1)/2}. \] (Testing that they are zero makes no sense: a singular \(\Sigma\) gives zero sample eigenvalues almost surely.)
For \(\mathbf{R}=\begin{pmatrix}1&r\\r&1\end{pmatrix}\) with \(r>0\), the eigenvalues are \(1\pm r\) with eigenvectors \((1,1)^\top/\sqrt2\) and \((1,-1)^\top/\sqrt2\). The first component is the standardised sum, explaining the fraction \((1+r)/2\); the second is the difference. The loadings do not depend on \(r\) at all, only the variance split does, which is a reminder that loadings alone say little without their eigenvalues.
The same algebra appears in finance. For asset returns with mean \(\boldsymbol\mu\) and covariance \(\Sigma\), the minimum-variance portfolio with weights summing to one is \(\mathbf{w}=\Sigma^{-1}\mathbf{1}/\mathbf{1}^\top\Sigma^{-1}\mathbf{1}\), obtained by a Lagrange multiplier on \(\min\mathbf{w}^\top\Sigma\mathbf{w}\) subject to \(\mathbf{w}^\top\mathbf{1}=1\). Adding a target return \(\mathbf{w}^\top\boldsymbol\mu=m\) traces out the efficient frontier. The leading principal component of a covariance matrix of stock returns is typically close to an equally weighted “market” portfolio.
Canonical correlation analysis
CCA generalises multiple correlation to two sets of variables: find the linear combination of each set that are most correlated with each other, then the next pair uncorrelated with the first, and so on. It reduces to an eigenproblem for a matrix built from both blocks.
Let \(\mathbf{X}^{(1)}\in\mathbb{R}^r\) and \(\mathbf{X}^{(2)}\in\mathbb{R}^q\) have joint covariance \(\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\Sigma_{21}&\Sigma_{22}\end{pmatrix}\). The first canonical variables are \(U_1=\mathbf{a}_1^\top\mathbf{X}^{(1)}\), \(V_1=\mathbf{b}_1^\top\mathbf{X}^{(2)}\) maximising \(\operatorname{cor}(U,V)\); the maximum \(\mu_1\) is the first canonical correlation. The \(k\)-th pair maximises the same correlation subject to being uncorrelated with all earlier \(U_j\) and \(V_j\).
Let \(\mu_1^2\ge\dots\ge\mu_q^2\) be the eigenvalues of \[ M=\Sigma_{22}^{-1/2}\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}\Sigma_{22}^{-1/2} \] with orthonormal eigenvectors \(\hat{\mathbf{b}}_k\). Then the \(k\)-th canonical correlation is \(\mu_k\), with \(\mathbf{b}_k=\Sigma_{22}^{-1/2}\hat{\mathbf{b}}_k\) and \(\mathbf{a}_k=\mu_k^{-1}\Sigma_{11}^{-1}\Sigma_{12}\mathbf{b}_k\). There are at most \(\min(r,q)\) non-zero canonical correlations.
Correlation is scale-free, so normalise \(\operatorname{Var}U=\operatorname{Var}V=1\) and substitute \(\mathbf{c}=\Sigma_{11}^{1/2}\mathbf{a}\), \(\mathbf{d}=\Sigma_{22}^{1/2}\mathbf{b}\), both unit vectors. Then \(\operatorname{cor}(U,V)=\mathbf{c}^\top K\mathbf{d}\) with \(K=\Sigma_{11}^{-1/2}\Sigma_{12}\Sigma_{22}^{-1/2}\). For fixed \(\mathbf{d}\), Cauchy-Schwarz gives \(\max_{\mathbf{c}}\mathbf{c}^\top K\mathbf{d}=\lVert K\mathbf{d}\rVert\), and \(\lVert K\mathbf{d}\rVert^2=\mathbf{d}^\top K^\top K\mathbf{d}\) with \(K^\top K=M\), so the Rayleigh-quotient theorem finishes it; the later pairs follow from the constrained version. Equivalently, the canonical correlations are the singular values of \(K\), so their number is at most \(\operatorname{rank}\Sigma_{12}\le\min(r,q)\).
With \(q=1\) there is one canonical correlation and it is the multiple correlation \(R\). Under joint normality, \(\mu_1=0\) if and only if \(\Sigma_{12}=0\), if and only if the two sets are independent; the union-intersection test of independence rejects for a large first sample canonical correlation \(r_1^2\), the largest eigenvalue of \(\mathbf{S}_{22}^{-1/2}\mathbf{S}_{21}\mathbf{S}_{11}^{-1}\mathbf{S}_{12}\mathbf{S}_{22}^{-1/2}\). In MANOVA, the eigenvalues \(\lambda_i\) of \(\mathbf{E}^{-1}\mathbf{H}\) are \(r_i^2/(1-r_i^2)\) for the canonical correlations between responses and the design.
The multivariate linear model and MANOVA
Stacking \(p\) regressions that share the same design gives a model whose estimator is just ordinary least squares column by column. The novelty is in testing: the error and hypothesis sums of squares become Wishart matrices, and there is no single best way to compare two matrices.
\[ \mathbf{Y}=\mathbf{X}B+\mathbf{E},\qquad \mathbf{Y}\in\mathcal{M}_{n,p},\ \mathbf{X}\in\mathcal{M}_{n,k}\text{ of full rank},\ B\in\mathcal{M}_{k,p}, \] with the rows of \(\mathbf{E}\) independent \(N_p(\mathbf{0},\Sigma)\). One-way MANOVA is the case where \(\mathbf{X}\) codes \(g\) groups.
The MLE (and least-squares estimator) is \(\hat B=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{Y}\), the univariate OLS fit applied to each response separately. With \(P=\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\), the error matrix \[ \mathbf{E}_{\mathrm{SS}}=\mathbf{Y}^\top(I-P)\mathbf{Y}\sim W_p(\Sigma,n-k), \] and \(\hat\Sigma=\mathbf{E}_{\mathrm{SS}}/n\). Correlation between responses changes the inference, never the point estimate.
\(\operatorname{tr}\big[\Sigma^{-1}(\mathbf{Y}-\mathbf{X}B)^\top(\mathbf{Y}-\mathbf{X}B)\big]\) is minimised over \(B\) by making \((\mathbf{Y}-\mathbf{X}B)^\top(\mathbf{Y}-\mathbf{X}B)\) as small as possible in the positive semi-definite order, and \((\mathbf{Y}-\mathbf{X}B)^\top(\mathbf{Y}-\mathbf{X}B)=\mathbf{Y}^\top(I-P)\mathbf{Y}+(\hat B-B)^\top\mathbf{X}^\top\mathbf{X}(\hat B-B)\). For the Wishart law, diagonalise the rank-\((n-k)\) projector \(I-P=\Gamma\operatorname{diag}(I_{n-k},0)\Gamma^\top\) and rotate the rows exactly as in the proof that \(\mathbf{V}\sim W_p(\Sigma,n-1)\). The regression analogue is on the Regression page.
A linear hypothesis \(H_0:CB=0\) on \(\nu_h\) degrees of freedom has a hypothesis matrix \(\mathbf{H}\), the increase in the residual matrix when \(H_0\) is imposed, with \(\mathbf{H}\sim W_p(\Sigma,\nu_h)\) under \(H_0\), independent of the error matrix \(\mathbf{E}\sim W_p(\Sigma,\nu_e)\). In one-way MANOVA \(\mathbf{H}=\sum_jn_j(\bar{\mathbf{y}}_j-\bar{\mathbf{y}})(\bar{\mathbf{y}}_j-\bar{\mathbf{y}})^\top\) (between groups) and \(\mathbf{E}=\sum_j\sum_i(\mathbf{y}_{ij}-\bar{\mathbf{y}}_j)(\mathbf{y}_{ij}-\bar{\mathbf{y}}_j)^\top\) (within groups), with \(\nu_h=g-1\), \(\nu_e=n-g\).
Let \(\lambda_1\ge\dots\ge\lambda_s>0\) be the non-zero eigenvalues of \(\mathbf{E}^{-1}\mathbf{H}\), \(s=\min(p,\nu_h)\). Every invariant test is a function of them; the standard four are \[
\begin{aligned} \text{Wilks:}&\quad \Lambda=\frac{|\mathbf{E}|}{|\mathbf{E}+\mathbf{H}|}=\prod_i\frac1{1+\lambda_i}\quad(\text{small rejects}),\\ \text{Pillai:}&\quad V=\operatorname{tr}\big[\mathbf{H}(\mathbf{H}+\mathbf{E})^{-1}\big]=\sum_i\frac{\lambda_i}{1+\lambda_i},\\ \text{Lawley-Hotelling:}&\quad U=\operatorname{tr}(\mathbf{E}^{-1}\mathbf{H})=\sum_i\lambda_i,\\ \text{Roy:}&\quad \lambda_1. \end{aligned}
\]
In one dimension \(\mathbf{H}\) and \(\mathbf{E}\) are numbers and the \(F\)-ratio is the obvious comparison. With matrices there is no total order, and each statistic summarises the \(s\) roots differently. Roy’s root is most powerful when the group means differ along a single direction; Pillai’s trace is the most robust to non-normality and unequal covariances; Wilks’s lambda is the likelihood-ratio test. They agree when \(s=1\).
Bartlett’s approximation: under \(H_0\), \[ -\Big[\nu_e-\tfrac12(p-\nu_h+1)\Big]\log\Lambda\ \dot\sim\ \chi^2_{p\nu_h}. \] When \(\nu_h=1\) (two groups, or one predictor) all four statistics are monotone functions of \(\lambda_1\), and \[ \frac{\nu_e-p+1}{p}\,\lambda_1\sim F_{p,\,\nu_e-p+1}\quad\text{exactly}, \] which is the two-sample Hotelling \(T^2\) test again. Otherwise software reports \(F\) approximations to each of the four.
Roy’s root has the same union-intersection meaning as \(T^2\): \(\lambda_1=\max_{\mathbf{a}}\mathbf{a}^\top\mathbf{H}\mathbf{a}/\mathbf{a}^\top\mathbf{E}\mathbf{a}\), the largest univariate \(F\)-ratio (up to degrees of freedom) over all linear combinations of the responses. The maximising \(\mathbf{a}\) is the first canonical discriminant direction, the projection that separates the groups best.
Tests on covariance matrices
Likelihood-ratio tests for structure in \(\Sigma\) all come from Anderson’s lemma: the maximised likelihood is a power of a determinant, so each statistic compares determinants of fitted and unrestricted covariance estimates.
For \(H_0:\Sigma=\Sigma_0\), first transform \(\mathbf{x}\mapsto\Sigma_0^{-1/2}\mathbf{x}\), which reduces to \(H_0:\Sigma=I_p\). Then \[ -2\log\Lambda=\operatorname{tr}\mathbf{V}-n\log|\mathbf{V}|+np\log n-np=n\sum_{j=1}^p(\ell_j-\log\ell_j-1), \] where \(\ell_j\) are the eigenvalues of \(\hat\Sigma=\mathbf{V}/n\). Under \(H_0\) it is asymptotically \(\chi^2_{p(p+1)/2}\).
The numerator maximises over \(\boldsymbol\mu\) only, giving \(e^{-\operatorname{tr}\mathbf{V}/2}\) up to constants; the denominator is the unrestricted maximum from Anderson’s lemma, \(|\mathbf{V}/n|^{-n/2}e^{-np/2}\). Take logs. Since \(x-\log x-1\ge0\) with equality only at \(x=1\), the statistic is zero exactly when \(\hat\Sigma=I\). The degrees of freedom count the free parameters of a symmetric \(p\times p\) matrix.
For \(H_0:\Sigma=\sigma^2I_p\) (spherical contours), \(\hat\sigma^2=\operatorname{tr}\mathbf{V}/(np)\) and \[ -2\log\Lambda=np\log\frac{\bar\ell}{\tilde\ell}, \] where \(\bar\ell\) and \(\tilde\ell\) are the arithmetic and geometric means of the eigenvalues of \(\hat\Sigma\); asymptotically \(\chi^2_{p(p+1)/2-1}\). The AM-GM inequality makes it non-negative. Sphericity is not scale-invariant, so it only makes sense when all variables share units.
For \(k\) groups of sizes \(n_i\), \(n=\sum n_i\), and \(H_0:\Sigma_1=\dots=\Sigma_k\), \[ -2\log\Lambda=n\log|\hat\Sigma_{\mathrm{p}}|-\sum_{i=1}^kn_i\log|\hat\Sigma_i|\ \dot\sim\ \chi^2_{(k-1)p(p+1)/2}. \] Box’s \(M\) corrects this for bias and for small samples: with \(N_i=n_i-1\), \(N=n-k\), unbiased \(\mathbf{S}_i\) and pooled \(\mathbf{S}_{\mathrm{p}}\), \[ M=N\log|\mathbf{S}_{\mathrm{p}}|-\sum_iN_i\log|\mathbf{S}_i|,\qquad \rho=1-\Big(\sum_i\frac1{N_i}-\frac1N\Big)\frac{2p^2+3p-1}{6(p+1)(k-1)}, \] and \(\rho M\dot\sim\chi^2_{(k-1)p(p+1)/2}\).
Take \(k=2\), \(p=2\), \(n_1=n_2=11\) and \(\mathbf{S}_1=\operatorname{diag}(4,1)\), \(\mathbf{S}_2=\operatorname{diag}(1,4)\). Both have determinant \(4\), but \(\mathbf{S}_{\mathrm{p}}=\operatorname{diag}(2.5,2.5)\) has determinant \(6.25\), so \[ M=20\log6.25-10\log4-10\log4=20\log\frac{6.25}{4}=8.93. \] Here \(\rho=1-(\tfrac1{10}+\tfrac1{10}-\tfrac1{20})\frac{8+6-1}{18}=1-0.15\times0.722=0.892\), so \(\rho M=7.96\) on \(3\) degrees of freedom, against a 5% critical value of \(7.81\): reject. Equal generalised variances do not make \(M\) small; the pooled matrix is “rounder” than either group’s, and \(M\) measures that via the determinant inequality \(|\tfrac12(A+B)|\ge|A|^{1/2}|B|^{1/2}\).
Box’s \(M\) is very sensitive to non-normality, particularly heavy tails, so a rejection may reflect kurtosis rather than unequal covariances. With large samples it rejects for differences that do not matter in practice. Use it as a diagnostic alongside plots, not as a gate.
Factor analysis
Factor analysis explains the covariance among \(p\) observed variables by a smaller number of unobserved factors. Unlike PCA it is a statistical model with an error term, so it can be fitted by maximum likelihood and tested, at the price of an unidentifiable rotation.
\[ \mathbf{X}=\boldsymbol\mu+\Lambda\mathbf{f}+\mathbf{e},\qquad \Lambda\in\mathcal{M}_{p,k},\ k<p, \] with common factors \(\mathbf{f}\sim N_k(\mathbf{0},I_k)\), specific errors \(\mathbf{e}\sim N_p(\mathbf{0},\Psi)\) with \(\Psi=\operatorname{diag}(\psi_1,\dots,\psi_p)\), and \(\mathbf{f}\) independent of \(\mathbf{e}\). The entries \(\lambda_{rj}\) are the loadings of variable \(r\) on factor \(j\).
The model implies \[ \Sigma=\Lambda\Lambda^\top+\Psi,\qquad \sigma_{rr}=\underbrace{\textstyle\sum_j\lambda_{rj}^2}_{\text{communality}}+\underbrace{\psi_r}_{\text{uniqueness}},\qquad \operatorname{Cov}(\mathbf{X},\mathbf{f})=\Lambda. \] For any orthogonal \(P\in\mathcal{M}_{k,k}\), the loadings \(\Lambda P\) with factors \(P^\top\mathbf{f}\) give the same \(\Sigma\), so \(\Lambda\) is identified only up to rotation.
\(\operatorname{Var}(\Lambda\mathbf{f}+\mathbf{e})=\Lambda I\Lambda^\top+\Psi\) because \(\mathbf{f}\) and \(\mathbf{e}\) are independent; and \((\Lambda P)(\Lambda P)^\top=\Lambda PP^\top\Lambda^\top=\Lambda\Lambda^\top\), with \(P^\top\mathbf{f}\sim N_k(\mathbf{0},I)\) still.
\(\Sigma\) has \(p(p+1)/2\) free entries; the model has \(pk+p\) parameters less \(k(k-1)/2\) removed by the rotation. The model is a genuine restriction only when \[ d=\tfrac12\big[(p-k)^2-(p+k)\big]>0. \] For \(p=6\) this allows at most \(k=3\) factors.
Fix \(\Psi\) and impose the identification condition that \(\Lambda^\top\Psi^{-1}\Lambda\) is diagonal. Let \(\omega_1\ge\dots\ge\omega_p\) be the eigenvalues of \(\Psi^{-1/2}\mathbf{S}\Psi^{-1/2}\) with unit eigenvectors \(\mathbf{u}_j\). The likelihood equations are solved by \[ \Lambda_{\cdot j}=\sqrt{\omega_j-1}\;\Psi^{1/2}\mathbf{u}_j,\qquad j=1,\dots,k, \] using the \(k\) largest eigenvalues (assumed greater than one). Alternating this step with an update of \(\Psi=\operatorname{diag}(\mathbf{S}-\Lambda\Lambda^\top)\) is Lawley’s algorithm.
The iteration can drive some \(\psi_r\) to zero or below (a Heywood or ultra-Heywood case): a variable claims all, or more than all, of its variance as common. It signals too many factors, too small a sample, or a variable that is nearly a linear combination of the others.
For \(H_0\): \(k\) factors suffice, against an unrestricted \(\Sigma\), \[ -2\log\Lambda=n\Big[\log|\hat\Sigma|-\log|\mathbf{S}_n|+\operatorname{tr}(\hat\Sigma^{-1}\mathbf{S}_n)-p\Big]\ \dot\sim\ \chi^2_d,\qquad \hat\Sigma=\hat\Lambda\hat\Lambda^\top+\hat\Psi, \] with \(d\) as above (Bartlett’s multiplier \(n-1-(2p+4k+5)/6\) in place of \(n\) improves the approximation). Fit \(k=1,2,\dots\) and stop at the first \(k\) not rejected.
Rotation is used to make the loadings interpretable. Varimax chooses the orthogonal \(P\) maximising the variance of squared loadings within each factor, \[ \sum_{j=1}^k\Big[\sum_{r=1}^p\lambda_{rj}^4-\frac1p\Big(\sum_{r=1}^p\lambda_{rj}^2\Big)^2\Big], \] pushing each loading towards either zero or large, so each factor is defined by a few variables.
The principal component solution takes \(\Lambda=[\sqrt{\lambda_1}\boldsymbol\gamma_1,\dots,\sqrt{\lambda_k}\boldsymbol\gamma_k]\) from \(\mathbf{S}\) or \(\mathbf{R}\) and sets \(\Psi\) from the diagonal remainder; the principal factor solution repeats this on \(\mathbf{R}-\hat\Psi\) with communalities on the diagonal (initialised at \(1-1/r^{rr}\), the squared multiple correlation of each variable on the rest). PCA explains total variance and has no model; factor analysis explains only the covariances, through a model that can be rejected. With small uniquenesses the two give similar loadings.
Structural equation models
Structural equation modelling turns factor analysis inside out: write down a model for how latent variables cause each other and how observed variables measure them, derive the covariance matrix it implies, and fit by matching it to \(\mathbf{S}\). Regression, factor analysis and path analysis are special cases.
A structural equation among latent variables and two measurement equations: \[ \boldsymbol\eta=B\boldsymbol\eta+\Gamma\boldsymbol\xi+\boldsymbol\zeta,\qquad \mathbf{Y}=\Lambda_Y\boldsymbol\eta+\boldsymbol\epsilon,\qquad \mathbf{X}=\Lambda_X\boldsymbol\xi+\boldsymbol\delta, \] where \(\boldsymbol\eta\) are endogenous and \(\boldsymbol\xi\) exogenous latent variables, \(I-B\) is non-singular, and the errors \(\boldsymbol\zeta,\boldsymbol\epsilon,\boldsymbol\delta\) are mean zero, mutually uncorrelated and uncorrelated with \(\boldsymbol\xi\). Write \(\Phi=\operatorname{Var}\boldsymbol\xi\), \(\Psi=\operatorname{Var}\boldsymbol\zeta\), \(\Theta_\epsilon\), \(\Theta_\delta\) for the error covariances, and \(\boldsymbol\theta\) for all free parameters.
With \(A=(I-B)^{-1}\), \[ \Sigma(\boldsymbol\theta)=\begin{pmatrix}\Lambda_YA(\Gamma\Phi\Gamma^\top+\Psi)A^\top\Lambda_Y^\top+\Theta_\epsilon & \Lambda_YA\Gamma\Phi\Lambda_X^\top\\ \Lambda_X\Phi\Gamma^\top A^\top\Lambda_Y^\top & \Lambda_X\Phi\Lambda_X^\top+\Theta_\delta\end{pmatrix}. \]
Solve the structural equation: \(\boldsymbol\eta=A(\Gamma\boldsymbol\xi+\boldsymbol\zeta)\), so \(\operatorname{Var}\boldsymbol\eta=A(\Gamma\Phi\Gamma^\top+\Psi)A^\top\) and \(\operatorname{Cov}(\boldsymbol\eta,\boldsymbol\xi)=A\Gamma\Phi\). Substitute into the measurement equations, using that the measurement errors are uncorrelated with everything else.
Under normality, \((n-1)\mathbf{S}\sim W(\Sigma,n-1)\), and maximising the Wishart likelihood is equivalent to minimising \[ F_{\mathrm{ML}}(\boldsymbol\theta)=\log|\Sigma(\boldsymbol\theta)|+\operatorname{tr}\big[\mathbf{S}\Sigma(\boldsymbol\theta)^{-1}\big]-\log|\mathbf{S}|-(p+q)\ \ge0, \] with equality exactly when \(\Sigma(\boldsymbol\theta)=\mathbf{S}\). The model is tested by \[ T=(n-1)F_{\mathrm{ML}}(\hat{\boldsymbol\theta})\ \dot\sim\ \chi^2_{\mathrm{df}},\qquad \mathrm{df}=\tfrac12(p+q)(p+q+1)-\dim\boldsymbol\theta. \]
(Non-negativity.) With \(\eta_j\) the eigenvalues of \(\Sigma^{-1/2}\mathbf{S}\Sigma^{-1/2}\), \(F_{\mathrm{ML}}=\sum_j(\eta_j-\log\eta_j-1)\), and each term is non-negative, zero only at \(\eta_j=1\). This is Anderson’s lemma once more.
Dropping the structural part (\(B=0\), \(\Gamma=0\)) leaves confirmatory factor analysis: unlike the exploratory version above, specific loadings are fixed at zero in advance, which removes the rotational indeterminacy and makes the hypothesis about which variable measures which factor testable. With no latent variables and \(\Lambda=I\), \(\Theta=0\) it is a system of simultaneous regressions (path analysis). In practice the \(\chi^2\) test rejects almost any model at large \(n\), so fit indices such as RMSEA and CFI are reported alongside it, and identification must be checked: a positive df is necessary but not sufficient.
Classification and discriminant analysis
Given labelled training data from \(g\) populations, allocate a new observation to one of them. The optimal rule minimises expected cost; under normality it is linear in \(\mathbf{x}\) when covariances are equal and quadratic when they are not.
Populations \(\pi_1,\dots,\pi_g\) have densities \(f_i\) and prior probabilities \(p_i\). A rule partitions \(\mathbb{R}^p\) into regions \(R_1,\dots,R_g\). With \(P(j\mid i)=\int_{R_j}f_i\) the probability of allocating a member of \(\pi_i\) to \(\pi_j\) and \(c(j\mid i)\) its cost, the expected cost of misclassification is \[ \mathrm{ECM}=\sum_{i=1}^gp_i\sum_{j\ne i}c(j\mid i)P(j\mid i). \] With unit costs it is the total probability of misclassification (TPM).
For \(g=2\) the ECM is minimised by \[ R_1=\Big\{\mathbf{x}:\ \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})}\ge\frac{c(1\mid2)}{c(2\mid1)}\cdot\frac{p_2}{p_1}\Big\}. \] With equal costs this is “allocate to the population with the larger posterior probability \(p_if_i(\mathbf{x})/\sum_jp_jf_j(\mathbf{x})\)”, and the same holds for \(g\) groups.
Since \(P(2\mid1)=1-\int_{R_1}f_1\), \[ \mathrm{ECM}=c(2\mid1)p_1+\int_{R_1}\big[c(1\mid2)p_2f_2(\mathbf{x})-c(2\mid1)p_1f_1(\mathbf{x})\big]\,d\mathbf{x}, \] which is smallest when \(R_1\) contains exactly the points where the integrand is non-positive.
If \(\pi_i=N_p(\boldsymbol\mu_i,\Sigma)\), the minimum-ECM rule allocates \(\mathbf{x}\) to \(\pi_1\) when \[ (\boldsymbol\mu_1-\boldsymbol\mu_2)^\top\Sigma^{-1}\mathbf{x}-\tfrac12(\boldsymbol\mu_1-\boldsymbol\mu_2)^\top\Sigma^{-1}(\boldsymbol\mu_1+\boldsymbol\mu_2)\ \ge\ \log\Big[\frac{c(1\mid2)}{c(2\mid1)}\cdot\frac{p_2}{p_1}\Big]. \] The boundary is a hyperplane. With unequal covariances the \(\mathbf{x}^\top\Sigma_i^{-1}\mathbf{x}\) terms no longer cancel, the left side gains \(-\tfrac12\mathbf{x}^\top(\Sigma_1^{-1}-\Sigma_2^{-1})\mathbf{x}\) and a \(\log(|\Sigma_1|/|\Sigma_2|)\) constant, and the boundary is a quadric (quadratic discriminant analysis).
Take logs of the density ratio: \(\log f_1-\log f_2=-\tfrac12(\mathbf{x}-\boldsymbol\mu_1)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu_1)+\tfrac12(\mathbf{x}-\boldsymbol\mu_2)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu_2)\). The \(\mathbf{x}^\top\Sigma^{-1}\mathbf{x}\) terms cancel and expanding the rest gives the left side.
Replacing parameters by \(\bar{\mathbf{x}}_1,\bar{\mathbf{x}}_2\) and the pooled \(\mathbf{S}_{\mathrm{p}}=\frac{(n_1-1)\mathbf{S}_1+(n_2-1)\mathbf{S}_2}{n_1+n_2-2}\) (unbiased for \(\Sigma\)) gives the sample rule based on \[ \hat{\mathbf{a}}^\top\mathbf{x}-\tfrac12\hat{\mathbf{a}}^\top(\bar{\mathbf{x}}_1+\bar{\mathbf{x}}_2),\qquad \hat{\mathbf{a}}=\mathbf{S}_{\mathrm{p}}^{-1}(\bar{\mathbf{x}}_1-\bar{\mathbf{x}}_2). \] Fisher derived the same direction without normality, as the \(\mathbf{a}\) maximising the separation \((\mathbf{a}^\top\bar{\mathbf{x}}_1-\mathbf{a}^\top\bar{\mathbf{x}}_2)^2/\mathbf{a}^\top\mathbf{S}_{\mathrm{p}}\mathbf{a}\) (Cauchy-Schwarz again).
With equal priors and costs and common \(\Sigma\), the minimum TPM is \[ \mathrm{OER}=\Phi\Big(-\frac\Delta2\Big),\qquad \Delta^2=(\boldsymbol\mu_1-\boldsymbol\mu_2)^\top\Sigma^{-1}(\boldsymbol\mu_1-\boldsymbol\mu_2), \] the squared Mahalanobis distance between the populations.
Let \(Y=\mathbf{l}^\top\mathbf{X}\) with \(\mathbf{l}=\Sigma^{-1}(\boldsymbol\mu_1-\boldsymbol\mu_2)\). Under \(\pi_i\), \(Y\) is normal with mean \(\mathbf{l}^\top\boldsymbol\mu_i\) and variance \(\mathbf{l}^\top\Sigma\mathbf{l}=\Delta^2\). The rule allocates to \(\pi_1\) when \(Y\ge m=\tfrac12\mathbf{l}^\top(\boldsymbol\mu_1+\boldsymbol\mu_2)\), the midpoint of the two means, which are \(\Delta^2\) apart. So \(P(2\mid1)=P(Y<m)=\Phi\big((-\Delta^2/2)/\Delta\big)=\Phi(-\Delta/2)\), and by symmetry \(P(1\mid2)\) is the same.
The apparent error rate, the misclassification rate on the training data, is optimistic. Use a test set or leave-one-out cross-validation (Lachenbruch’s holdout). For two classes the confusion matrix gives sensitivity and specificity; sweeping the threshold on the right of the rule traces the ROC curve.
Support vector machines
SVMs drop the distributional model entirely and choose the separating hyperplane with the widest margin. They come from statistical learning theory, which bounds the true risk by the empirical risk plus a complexity term. The ML view is on the SVM page.
For a classifier \(f\) with zero-one loss, the risk is \(R(f)=P(f(\mathbf{X})\ne Y)\) and the empirical risk is the training error \(\hat R(f)\). If the class \(\mathcal{F}\) has VC dimension \(h\), then with probability \(1-\eta\), simultaneously for all \(f\in\mathcal{F}\), \[ R(f)\le\hat R(f)+\phi\Big(\frac hn,\frac{\log\eta}n\Big), \] with \(\phi\) increasing in \(h\). Linear rules in \(\mathbb{R}^p\) have \(h=p+1\). Minimising the bound, not \(\hat R\) alone, is structural risk minimisation; the margin is a way to control complexity independently of \(p\).
Training data \((\mathbf{x}_i,y_i)\) with \(y_i\in\{-1,+1\}\) are linearly separable if some \((\mathbf{w},b)\) has \(y_i(\mathbf{w}^\top\mathbf{x}_i+b)>0\) for all \(i\). Fixing the scale by \(\min_iy_i(\mathbf{w}^\top\mathbf{x}_i+b)=1\), the distance from the hyperplane to the nearest point is \(1/\lVert\mathbf{w}\rVert\), and the SVM solves \[ \min_{\mathbf{w},b}\ \tfrac12\lVert\mathbf{w}\rVert^2\quad\text{subject to}\quad y_i(\mathbf{w}^\top\mathbf{x}_i+b)\ge1,\ i=1,\dots,n. \]
The KKT conditions give \(\mathbf{w}=\sum_i\alpha_iy_i\mathbf{x}_i\) and \(\sum_i\alpha_iy_i=0\), with \(\alpha_i\big[y_i(\mathbf{w}^\top\mathbf{x}_i+b)-1\big]=0\). The multipliers solve the quadratic programme \[ \max_{\boldsymbol\alpha}\ \sum_i\alpha_i-\tfrac12\sum_{i,j}\alpha_i\alpha_jy_iy_j\,\mathbf{x}_i^\top\mathbf{x}_j\quad\text{subject to}\quad\alpha_i\ge0,\ \sum_i\alpha_iy_i=0. \] Only points on the margin, the support vectors, have \(\alpha_i>0\). The soft-margin version, which adds slacks \(\xi_i\ge0\) and the penalty \(C\sum_i\xi_i\), has the same dual with the box constraint \(0\le\alpha_i\le C\).
The Lagrangian is \(L=\tfrac12\lVert\mathbf{w}\rVert^2-\sum_i\alpha_i[y_i(\mathbf{w}^\top\mathbf{x}_i+b)-1]\). Setting \(\partial L/\partial\mathbf{w}=0\) and \(\partial L/\partial b=0\) gives the two stationarity conditions; substituting back, \(\tfrac12\lVert\mathbf{w}\rVert^2-\mathbf{w}^\top\mathbf{w}=-\tfrac12\lVert\mathbf{w}\rVert^2\), which is the quadratic term. Complementary slackness forces \(\alpha_i=0\) unless the constraint is active. For the soft margin, stationarity in \(\xi_i\) gives \(C-\alpha_i-\mu_i=0\) with \(\mu_i\ge0\), hence \(\alpha_i\le C\), and the \(\xi_i\) terms drop out.
The data enter the dual and the predictor \(\hat f(\mathbf{x})=\operatorname{sign}\big(\sum_i\alpha_iy_i\mathbf{x}_i^\top\mathbf{x}+b\big)\) only through inner products. Replacing \(\mathbf{x}_i^\top\mathbf{x}_j\) by a kernel \(K(\mathbf{x}_i,\mathbf{x}_j)\) fits a linear separator in a feature space without computing the features, giving non-linear boundaries in the original space.
Cluster analysis
Clustering is classification without labels: find groups from the data alone. There is no likelihood to anchor it, so the answer depends on the distance, the algorithm, and the number of clusters, all of which the analyst chooses. See also the k-means and Gaussian mixtures pages.
For a partition into \(S_1,\dots,S_K\) with centroids \(\bar{\mathbf{x}}_{S_k}\), the within-cluster sum of squares is \[ W=\sum_{k=1}^K\sum_{i\in S_k}\lVert\mathbf{x}_i-\bar{\mathbf{x}}_{S_k}\rVert^2=\sum_{k=1}^K\frac1{2|S_k|}\sum_{i,j\in S_k}\lVert\mathbf{x}_i-\mathbf{x}_j\rVert^2. \] \(K\)-means alternates: assign each point to its nearest centroid, then recompute centroids. Each step cannot increase \(W\), so it converges, but only to a local minimum; restart from several initialisations. \(K\)-medoids (PAM) uses actual data points as centres and any dissimilarity, which makes it robust to outliers.
(Pairwise form.) For any finite set with mean \(\bar{\mathbf{x}}\), \(\sum_{i,j}\lVert\mathbf{x}_i-\mathbf{x}_j\rVert^2=\sum_{i,j}\lVert(\mathbf{x}_i-\bar{\mathbf{x}})-(\mathbf{x}_j-\bar{\mathbf{x}})\rVert^2=2m\sum_i\lVert\mathbf{x}_i-\bar{\mathbf{x}}\rVert^2\), with \(m\) the size of the set, because the cross terms sum to zero.
Agglomerative clustering starts with \(n\) singletons and repeatedly merges the two closest clusters, recording the merge heights in a dendrogram; cutting it at a height gives a partition. The linkage defines the distance between clusters: single (closest pair; chains), complete (farthest pair; compact), average (mean pair), or Ward (increase in \(W\)). All fit the Lance-Williams recursion \[ d(S_1\cup S_2,S_3)=\alpha_1d(S_1,S_3)+\alpha_2d(S_2,S_3)+\beta d(S_1,S_2)+\gamma|d(S_1,S_3)-d(S_2,S_3)|, \] so the distance matrix is updated after each merge without revisiting the data. Single linkage has \(\alpha_1=\alpha_2=\tfrac12\), \(\beta=0\), \(\gamma=-\tfrac12\), which is the minimum.
Merging \(S_1\) and \(S_2\), of sizes \(n_1,n_2\), increases \(W\) by \[ \Delta W=\frac{n_1n_2}{n_1+n_2}\lVert\bar{\mathbf{x}}_{S_1}-\bar{\mathbf{x}}_{S_2}\rVert^2. \]
The merged centroid is \(\bar{\mathbf{x}}=(n_1\bar{\mathbf{x}}_1+n_2\bar{\mathbf{x}}_2)/(n_1+n_2)\). By the parallel-axis decomposition, \(\sum_{i\in S_1}\lVert\mathbf{x}_i-\bar{\mathbf{x}}\rVert^2=\sum_{i\in S_1}\lVert\mathbf{x}_i-\bar{\mathbf{x}}_1\rVert^2+n_1\lVert\bar{\mathbf{x}}_1-\bar{\mathbf{x}}\rVert^2\), and likewise for \(S_2\). So \(\Delta W=n_1\lVert\bar{\mathbf{x}}_1-\bar{\mathbf{x}}\rVert^2+n_2\lVert\bar{\mathbf{x}}_2-\bar{\mathbf{x}}\rVert^2\). Since \(\bar{\mathbf{x}}_1-\bar{\mathbf{x}}=\frac{n_2}{n_1+n_2}(\bar{\mathbf{x}}_1-\bar{\mathbf{x}}_2)\) and \(\bar{\mathbf{x}}_2-\bar{\mathbf{x}}=-\frac{n_1}{n_1+n_2}(\bar{\mathbf{x}}_1-\bar{\mathbf{x}}_2)\), this is \(\frac{n_1n_2^2+n_2n_1^2}{(n_1+n_2)^2}\lVert\bar{\mathbf{x}}_1-\bar{\mathbf{x}}_2\rVert^2\), as claimed.
For point \(i\) in cluster \(G_i\), let \(a(i)\) be its average distance to the other members of \(G_i\) and \(b(i)\) its smallest average distance to the members of another cluster. Its silhouette is \[ s(i)=\frac{b(i)-a(i)}{\max\{a(i),b(i)\}}\in[-1,1] \] (zero for a singleton). Values near \(1\) mean well clustered, negative values mean the point sits nearer another cluster. The average silhouette width is used to choose \(K\).
Hide the species of Fisher’s 150 iris flowers (four measurements, three species of 50) and cluster. Ward’s method, \(K\)-means and \(K\)-medoids all give nearly the same three clusters: setosa recovered exactly, and about 15 virginica placed with versicolor, for 16 misclustered flowers (10.7%). Average silhouette width at \(K=3\) is \(0.55\) for all three. Yet the silhouette prefers \(K=2\) (width about \(0.69\)): setosa is far from the other two species, and versicolor and virginica overlap where the density is high, so a distance-based method has no reason to cut there. The lesson is that the “right” \(K\) is a scientific question the data may not answer.
Copulas
Every joint distribution is its margins glued together by a copula, a joint distribution on the unit cube with uniform margins. Separating the two lets one model heavy-tailed margins and non-normal dependence, including dependence in the tails, which the normal cannot capture.
A \(p\)-dimensional copula is a joint distribution function \(C\) on \([0,1]^p\) whose margins are uniform on \([0,1]\).
If \(X\) has continuous distribution function \(F\), then \(F(X)\sim U(0,1)\). Conversely, if \(U\sim U(0,1)\) then \(F^{-1}(U)\) has distribution function \(F\), where \(F^{-1}\) is the generalised inverse.
For continuous strictly increasing \(F\), \(P(F(X)\le u)=P(X\le F^{-1}(u))=F(F^{-1}(u))=u\); the general continuous case uses the generalised inverse, and the converse follows from \(F^{-1}(u)\le x\iff u\le F(x)\).
For any joint distribution function \(F\) with margins \(F_1,\dots,F_p\) there is a copula \(C\) with \[ F(x_1,\dots,x_p)=C\big(F_1(x_1),\dots,F_p(x_p)\big)\quad\text{for all }\mathbf{x}, \] unique when the margins are continuous. Conversely, for any copula \(C\) and distribution functions \(F_j\), the right side is a joint distribution function with margins \(F_j\). If densities exist, \(f(\mathbf{x})=c\big(F_1(x_1),\dots,F_p(x_p)\big)\prod_jf_j(x_j)\) with \(c=\partial^pC/\partial u_1\cdots\partial u_p\).
(Continuous margins.) Let \(U_j=F_j(X_j)\), uniform by the lemma, and let \(C\) be the joint distribution function of \((U_1,\dots,U_p)\), a copula. Since \(F_j\) is non-decreasing and \(X_j=F_j^{-1}(U_j)\) almost surely, \(\{X_j\le x_j\}=\{U_j\le F_j(x_j)\}\) up to null sets, so \(F(\mathbf{x})=C(F_1(x_1),\dots,F_p(x_p))\). Uniqueness: \(C\) is pinned down on \(\prod_j\operatorname{range}(F_j)=[0,1]^p\). For the converse, take \(\mathbf{U}\sim C\) and set \(X_j=F_j^{-1}(U_j)\). The density formula is the chain rule.
- Independence: \(C(\mathbf{u})=\prod_ju_j\).
- Gaussian: \(C_R(\mathbf{u})=\Phi_R\big(\Phi^{-1}(u_1),\dots,\Phi^{-1}(u_p)\big)\) for a correlation matrix \(R\): the dependence of a multivariate normal, without its margins. It has no tail dependence.
- \(t\): the same construction from the multivariate \(t_\nu\); it has symmetric tail dependence that grows as \(\nu\) falls.
- Archimedean: \(C(\mathbf{u})=\phi^{-1}\big(\phi(u_1)+\dots+\phi(u_p)\big)\) for a continuous, strictly decreasing, convex generator \(\phi:[0,1]\to[0,\infty]\) with \(\phi(1)=0\) (and \(p\)-monotone inverse when \(p>2\)). The Gumbel-Hougaard copula, \(\phi(t)=(-\log t)^\theta\), \(\theta\ge1\), gives
\[ C_\theta(\mathbf{u})=\exp\Big\{-\Big[\sum_j(-\log u_j)^\theta\Big]^{1/\theta}\Big\}, \] independence at \(\theta=1\) and upper-tail dependence increasing with \(\theta\); Clayton is the lower-tail counterpart.
Fit in two stages: estimate each margin (parametrically, or by the empirical distribution function rescaled by \(n/(n+1)\), giving pseudo-observations \(\hat u_{ij}\)), then maximise the copula likelihood \(\sum_i\log c(\hat{\mathbf{u}}_i;\boldsymbol\theta)\). To simulate, draw \(\mathbf{U}\sim C\) and transform each coordinate by \(F_j^{-1}\). Rank correlations (Kendall’s \(\tau\), Spearman’s \(\rho\)) depend only on the copula, which is why they, not Pearson’s \(r\), are the natural dependence measures here.
Where it leads
- High dimensions. When \(p\) is comparable to \(n\), \(\mathbf{S}\) is singular or badly conditioned, its eigenvalues spread out (Marchenko-Pastur), and \(T^2\) and LDA break down. Shrinkage estimators of \(\Sigma\), sparse PCA and regularised discriminant analysis are the responses.
- Graphical models. Zeros in \(\Sigma^{-1}\) are conditional independences under normality; estimating a sparse precision matrix (the graphical lasso) is partial correlation at scale.
- Bayesian multivariate analysis. The Wishart is the conjugate prior for the precision matrix of a normal, and the normal-inverse-Wishart family makes the MLE results above the large-sample limit of the posterior.
- Mixtures and model-based clustering. Replacing the \(K\)-means criterion by a likelihood for a mixture of normals with structured covariances gives clustering with a principled choice of \(K\) (by BIC); see Gaussian mixtures.
- Time series and spatial data. Vector autoregressions and Gaussian processes are multivariate normals with structured covariances indexed by time or space.
Backlinks (2)
1. Bayesian Inference /wiki/mathematics/statistics/bayesian/
2. Statistics /wiki/mathematics/statistics/
This page pairs well with Probability.
Two child pages go further:
- Bayesian Inference: priors and posteriors, decisions, Monte Carlo and MCMC, Bayes factors, hierarchical and mixture models.
- Multivariate Analysis: the multivariate normal and Wishart distributions, Hotelling’s \(T^2\), tests on mean vectors and covariances, PCA, MANOVA, factor analysis, classification and copulas.
Table of Distributions
| Distribution | mass/density function | \[S_X\] | \[\mathbb{E}(X)\] | \[\mathrm{Var}(X)\] | \[\phi_X(s)\] |
|---|---|---|---|---|---|
| Bernoulli \[\mathrm{Bern}(\pi)\] | \[P(X=1)=\pi\\ P(X=0)=1-\pi\] | \[\{0,1\}\] | \[\pi\] | \[\pi(1-\pi)\] | \[(1-\pi)+\pi e^{s}\] |
| Binomial \[\mathrm{Bin}(n,\pi)\] | \[p_X(x)=\binom{n}{x}\pi^{x}(1-\pi)^{n-x}\] | \[\{0,1,\dots,n\}\] | \[n\pi\] | \[n\pi(1-\pi)\] | \[(1-\pi+\pi e^{s})^{n}\] |
| Geometric \[\mathrm{Geo}(\pi)\] | \[p_X(x)=\pi(1-\pi)^{x-1}\] | \[\{1,2,\dots\}\] | \[\pi^{-1}\] | \[(1-\pi)\pi^{-2}\] | \[\frac{\pi}{e^{-s}-1+\pi}\] |
| Poisson \[\mathcal{P}(\lambda)\] | \[p_X(x)=e^{-\lambda}\lambda^{x}/x!\] | \[\{0,1,\dots\}\] | \[\lambda\] | \[\lambda\] | \[\exp\{\lambda(e^{s}-1)\}\] |
| Uniform \[U[\alpha,\beta]\] | \[f_X(x)=(\beta-\alpha)^{-1}\] | \[[\alpha,\beta]\] | \[\frac{1}{2}(\alpha+\beta)\] | \[\frac{1}{12}(\beta-\alpha)^2\] | \[\frac{e^{\beta s}-e^{\alpha s}}{s(\beta-\alpha)}\] |
| Exponential \[\mathrm{Exp}(\lambda)\] | \[f_X(x)=\lambda e^{-\lambda x}\] | \[[0,\infty)\] | \[\lambda^{-1}\] | \[\lambda^{-2}\] | \[\frac{\lambda}{\lambda-s}\] |
| Gaussian \[\mathcal{N}(\mu,\sigma^{2})\] | \[f_X(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^{2}}\right\}\] | \[\mathbb{R}\] | \[\mu\] | \[\sigma^{2}\] | \[e^{\mu s+\frac{1}{2}\sigma^{2}s^{2}}\] |
| Gamma \[\Gamma(\alpha,\lambda)\] | \[f_X(x)=\frac{1}{\Gamma(\alpha)}\lambda^{\alpha}x^{\alpha-1}e^{-\lambda x}\] | \[[0,\infty)\] | \[\alpha\lambda^{-1}\] | \[\alpha\lambda^{-2}\] | \[\left(\frac{\lambda}{\lambda-s}\right)^{\alpha}\] |
Statistical Inference
Let \(X=(X_1,\ldots,X_n)\) be i.i.d. from a parametric family \(\{F_\theta:\theta\in\Theta\subset\mathbb{R}^p\}\). The parameter \(\theta\) is unknown; inference uses the randomness of \(X\) to learn about \(\theta\).
Comments