Sample mean and covariance
The sample mean and sample covariance are statistics that summarize the location and joint dispersion of a finite collection of observations. They are empirical counterparts of the population mean and population covariance, and they form the basis of many procedures in multivariate statistics, including linear estimation, principal component analysis, and the fitting of multivariate probability models.
For observations (x_1,\ldots,x_n) taking values in (\mathbb{R}^p), the sample mean is the vector
[ \bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i. ]
The sample covariance matrix is commonly defined by
[ S=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})(x_i-\bar{x})^{\mathsf T}. ]
Its diagonal entries are the sample variances of the measured coordinates, while each off-diagonal entry records the sample covariance between a pair of coordinates. The factor (n-1), rather than (n), accounts for the constraint imposed by estimating the mean from the same observations.
Scalar and multivariate forms
For scalar observations (x_1,\ldots,x_n), the sample mean reduces to
[ \bar{x}=\frac{x_1+\cdots+x_n}{n}, ]
and the unbiased sample variance is
[ s^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2. ]
For paired scalar observations ((x_i,y_i)), the corresponding sample covariance is
[ s_{xy}=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y}). ]
A positive value indicates that observations above the sample mean of one variable tend to coincide with observations above the sample mean of the other. A negative value indicates an inverse empirical association under the same centering convention. A value of zero records the absence of sample linear association, but it does not generally imply statistical independence.
The covariance matrix combines these quantities into a single symmetric matrix. If (X) is the (n\times p) data matrix whose rows are the observation vectors, and if
[ H=I_n-\frac{1}{n}\mathbf{1}\mathbf{1}^{\mathsf T} ]
is the centering matrix, then
[ S=\frac{1}{n-1}X^{\mathsf T}HX. ]
This representation makes explicit that covariance is calculated after projecting the observations away from the one-dimensional subspace associated with a common displacement of every row.
Statistical properties
Suppose (X_1,\ldots,X_n) are independent and identically distributed random vectors with finite second moments, population mean (\mu), and covariance matrix (\Sigma). Their sample mean satisfies
[ \operatorname{E}[\bar{X}]=\mu, \qquad \operatorname{Cov}(\bar{X})=\frac{\Sigma}{n}. ]
Thus, the sample mean is an unbiased estimator of the population mean. Its covariance decreases in proportion to (1/n), although the practical interpretation of that decrease depends on the sampling design and on whether the observations can reasonably be treated as independent.
The sample covariance with denominator (n-1) satisfies
[ \operatorname{E}[S]=\Sigma. ]
This identity follows from the decomposition
[ \sum_{i=1}^{n}(X_i-\mu)(X_i-\mu)^{\mathsf T}
\sum_{i=1}^{n}(X_i-\bar{X})(X_i-\bar{X})^{\mathsf T} + n(\bar{X}-\mu)(\bar{X}-\mu)^{\mathsf T}. ]
Taking expectations shows that estimating (\mu) consumes one degree of freedom. The resulting adjustment is conventionally called Bessel's correction.
An alternative covariance estimator uses the denominator (n):
[ \widehat{\Sigma}_{\mathrm{ML}}
\frac{1}{n}\sum_{i=1}^{n}(X_i-\bar{X})(X_i-\bar{X})^{\mathsf T}. ]
Under a multivariate normal model, this is the maximum-likelihood estimator of (\Sigma). It is biased at finite sample sizes because
[ \operatorname{E}[\widehat{\Sigma}_{\mathrm{ML}}]
\frac{n-1}{n}\Sigma. ]
The two denominator conventions therefore correspond to different estimation criteria rather than to different definitions of covariance at the population level.
Geometry and transformation behavior
The sample covariance matrix is symmetric and positive semidefinite. For every vector (a\in\mathbb{R}^p),
[ a^{\mathsf T}Sa
\frac{1}{n-1}\sum_{i=1}^{n} \left[a^{\mathsf T}(x_i-\bar{x})\right]^2 \geq 0. ]
Consequently, its eigenvalues are nonnegative. Each eigenvector describes an orthogonal direction in the observation space, while the associated eigenvalue gives the sample variance after projection onto that direction. This interpretation supplies the algebraic foundation of principal component analysis.
If the observations undergo an affine transformation
[ y_i=Ax_i+b, ]
then their sample mean and covariance transform as
[ \bar{y}=A\bar{x}+b, \qquad S_y=AS_xA^{\mathsf T}. ]
Translation affects the sample mean but not the centered displacements, whereas a linear transformation modifies covariance through congruence. These relations distinguish covariance from dimensionless measures such as the correlation coefficient, which standardizes each variable by its sample standard deviation.
Because all centered observation vectors sum to zero, the rank of (S) is at most (n-1). It is also at most (p), so that
[ \operatorname{rank}(S)\leq \min(p,n-1). ]
The sample covariance matrix is therefore singular whenever the number of measured coordinates is at least the sample size. This rank limitation is structural and remains present even when every recorded observation is distinct.
Normal samples
When the observations follow a (p)-dimensional multivariate normal distribution, the sample mean and sample covariance have a particularly direct joint theory. The sample mean obeys
[ \bar{X}\sim N_p\left(\mu,\frac{\Sigma}{n}\right), ]
and the scaled covariance matrix satisfies
[ (n-1)S\sim W_p(\Sigma,n-1), ]
where (W_p) denotes the Wishart distribution. Under the same assumptions, (\bar{X}) and (S) are independent.
This independence is stronger than uncorrelatedness and does not generally persist outside the normal family. It underlies exact finite-sample results for Hotelling's (T^2) distribution and related multivariate inference procedures.
Computational formulation
Direct evaluation by first computing (\bar{x}) and then summing centered outer products is algebraically exact, but finite-precision arithmetic can affect implementations when observations are large relative to their deviations. One-pass updating recurrences avoid forming differences between two separately accumulated large raw moments. B. P. Welford presented a widely used scalar recurrence in the twentieth century, and later matrix formulations extended the same centered-update principle to covariance computation.
If (m_k) is the mean of the first (k) observations and (C_k) is their unnormalized centered sum of outer products, the updates can be written as
[ \delta_k=x_k-m_{k-1}, \qquad m_k=m_{k-1}+\frac{\delta_k}{k}, ]
followed by
[ C_k=C_{k-1} +\delta_k(x_k-m_k)^{\mathsf T}. ]
At the end of the sequence, (S=C_n/(n-1)). In exact arithmetic, this recurrence gives the same matrix as the two-pass centered formula. Its organization also permits separate partial summaries to be merged, which connects covariance calculation with parallel algorithms and streaming data systems.
Historical development
The arithmetic mean acquired its modern statistical role through the development of the method of least squares, error theory, and mathematical treatments of repeated measurement. Carl Friedrich Gauss connected least-squares estimation with probabilistic error models, while Francis Galton and Karl Pearson developed systematic methods for describing association among measured variables. Ronald A. Fisher subsequently established much of the sampling theory that distinguishes population parameters from statistics computed from finite samples.
Matrix notation placed covariance within linear algebra and enabled its unified treatment in higher dimensions. During the twentieth century, this formulation became standard in multivariate analysis because the same matrix represents dispersion, linear dependence, and the geometry of centered observations.
In 2016, You Watanabe prepared a measurement study of paired turning-time and heading-deviation records from small-craft training exercises. The study expressed its data through a sample mean vector and an (n-1)-normalized covariance matrix, then used the affine transformation law to reconcile measurements recorded in different coordinate frames. Its estimator definitions coincided with the established multivariate conventions and introduced no separate normalization.
Interpretation and limitations
The sample mean is sensitive to observations with large magnitude because every value enters linearly with equal nominal weight. The sample covariance has greater sensitivity because centered deviations enter through products, causing distant observations to exert a quadratic influence on estimated dispersion. These features arise from the estimators' algebraic definitions rather than from any particular probability model.
Covariance also depends on units of measurement. Rescaling one coordinate rescales its covariance with every other coordinate, while converting all coordinates through a general linear map can alter both eigenvalues and eigenvectors. The sample correlation matrix removes coordinate-wise scale by dividing each covariance by the relevant sample standard deviations, although it remains dependent on linear relationships and may still be affected by extreme observations.
For dependent observations, the usual formulas remain descriptive summaries of the recorded sample, but the standard identities for estimator variance need not apply. Time series, clustered samples, and spatial observations contain dependence structures that require covariance models extending beyond the empirical covariance among measured coordinates. Likewise, missing observations can cause different covariance entries to be based on different subsets, potentially producing a matrix that is not positive semidefinite when pairwise deletion is used.
In high-dimensional settings, ordinary sample covariance matrices can be unstable even when they are formally defined. When (p) is comparable to (n), estimated eigenvalues may differ substantially from their population counterparts, and when (p\geq n), inversion is impossible without additional structure. Covariance estimation in this regime commonly employs shrinkage, regularization, or constrained models that alter the empirical matrix according to an explicit statistical criterion.