- Estimation of covariance matrices
statistics, sometimes the covariance matrixof a multivariate random variableis not known but has to be estimated. Estimation of covariance matrices then deals with the question of how to approximate the actual covariance matrix on the basis of a sample from the multivariate distribution. Simple cases, where observations are complete, can be dealt with by using the sample covariance matrix. The sample covariance matrix is an unbiased estimatorof the covariance matrix. In addition, if the random variable has normal distribution, the sample covariance matrix has Wishart distributionand a slightly differently scaled version of it is the maximum likelihood estimate. Cases involving missing datarequire deeper considerations. Statistical analyses of multivariate data often involve exploratory studies of the way in which the variables change in relation to one another and this may be followed up by explicit statistical models involving the covariance matrix of the variables. Thus the estimation of covariance matrices directly from observational data plays two roles::* to provide initial estimates that can be used to study the inter-relationships;:* to provide sample estimates that can be used for model checking.
Estimates of covariance matrices are required at the intial stages of
principal component analysisand factor analysis, and are also involved in versions of regression analysisthat treat the dependent variables in a data-set, jointly with the independent variableas the outcome of a random sample.
Estimation in a general context
is the sample covariance matrix
is the sample mean.This is true regardless of the distribution of the random vector "X", provided of course that the theoretical means and covariances exist. The reason for the factor "n" − 1 rather than "n" is essentially the same as the reason for the same factor appearing in unbiased estimates of sample variances and covariances, which relates to the fact that the mean is not known and is replaced by the sample mean.
In cases where the distribution of the
random variable"X" is known to be within a certain family of distributions, other estimate may be derived on the basis of that assumption. A well-known instance is when the random variable"X" is normally distributed: in this case the maximum likelihood estimatorof the covariance matrix is slightly different from the unbiased estimate, and is given by
A derivation of this result is given below. Clearly, the difference between the unbiased estimator and the maximum likelihood estimator diminishes for large "n".
In the general case, the unbiased estimate of the covariance matrix provides an acceptable estimate when the data vectors in the observed data set are all complete: that is they contain no missing elements. One approach to estimating the covariance matrix is to treat the estimation of each variance or pairwise covariance separately, and to use all the observations for which both variables have valid values. Assuming the missing data are
missing at randomthis results in an estimate for the covariance matrix which is unbiased. However, for many applications this estimate may not be acceptable because the estimated covariance matrix is not guaranteed to be positive semi-definite. This could lead to estimated correlations having absolute values which are greater than one.
Maximum-likelihood estimation for the multivariate normal distribution
A random vector "X" ∈ R"p"×1 (a "p"×1 "column vector") has a multivariate normal distribution with a nonsingular covariance matrix "Σ" precisely if "Σ" ∈ R"p" × "p" is a
positive-definite matrixand the probability density functionof "X" is
where μ ∈ R"p"×1 is the
expected value. The matrix "Σ" is the higher-dimensional analog of what in one dimension would be the variance.
Suppose now that "X"1, ..., "X""n" are independent and identically distributed with the distribution above. Based on the
observed values "x"1, ..., "x""n" of this sample, we wish to estimate "Σ". (Note that we adhere to the convention of writing random variables as capital letters and data as lower-case letters).
It is fairly readily shown that the maximum-likelihood estimate of the mean vector μ is the "
sample mean" vector:
See the section on estimation in the article on the normal distribution for details; the process here is similar.
Since the estimate of μ does not depend on "Σ", we can just substitute it for μ in the
and then seek the value of "Σ" that maximizes this (in practice it is easier to work with ).
The trace of a 1 × 1 matrix
Now we come to the first surprising step.
scalaras the trace of a 1×1 matrix!
This makes it possible to use the identity tr("AB") = tr("BA") whenever "A" and "B" are matrices so shaped that both products exist.We get
(so now we are taking the trace of a "p"×"p" matrix!)
: is sometimes called the
Using the spectral theorem
It follows from the
spectral theoremof linear algebrathat a positive-definite symmetric matrix "S" has a unique positive-definite symmetric square root "S"1/2. We can again use the "cyclic property" of the trace to write
Let "B" = "S"1/2 "Σ"−1 "S"1/2. Then the expression above becomes
The positive-definite matrix "B" can be diagonalized, and then the problem of finding the value of "B" that maximizes
reduces to the problem of finding the values of the diagonal entries λ1, ..., λ"p" that maximize
This is just a calculus problem and we get λ"i" = "n", so that "B" = "n Ip", i.e., "n" times the "p"×"p" identity matrix.
Finally we get
i.e., the "p"×"p" "sample covariance matrix"
is the maximum-likelihood estimator of the "population covariance matrix" "Σ". At this point we are using a capital "X" rather than a lower-case "x" because we are thinking of it "as an estimator rather than as an estimate", i.e., as something random whose probability distribution we could profit by knowing. The random matrix "S" can be shown to have a
Wishart distributionwith "n" − 1 degrees of freedom. [
K.V. Mardia, J.T. Kent, and J.M. Bibby (1979) "
Multivariate Analysis", Academic Press.] That is:
An alternative derivation of the maximum likelihood estimator can be performed via
matrix calculusformulae (see also differential of a determinant and differential of the inverse matrix). It also verifies the aforementioned fact about the maximum likelihood estimate of the mean. Re-write the likelihood in the log form using the trace trick:
The differential of this log-likelihood is
It naturally breaks down into the part related to the estimation of the mean, and to the part related to the estimation of the variance. The first order condition for maximum, , is satisfied when the terms multiplying and are identically zero. Assuming (the maximum likelihood estimate of) is non-singular, the first order condition for the estimate of the mean vector is
which leads to the maximum likelihood estimator
This lets us simplify as defined above. Then the terms involving in can be combined as
The first order condition will hold when the term in the square bracket is (matrix-valued) zero. Pre-multiplying the latter by and dividing by gives
which of course coincides with the canonical derivation given earlier.
Dwyer cite journal| title=Some applications of matrix derivatives in multivariate analysis| journal=Journal of the American Statistical Association| month=June| year=1967| volume=62| number=318 |pages=607–625| doi=10.2307/2283988| author=Dwyer, Paul S.] points out that decomposition into two terms such as appears above is "unnecessary" and derives the estimator in two lines of working.
If the sample size "n" is small and the number of considered variables "p" is large, the above empirical estimators of covariance and correlation are very unstable. Specifically, it is possible to furnish estimators that improve considerably upon the maximum likelihood estimate in terms of mean squared error. Moreover, for "n" < "p", the empirical estimate of the covariance matrix becomes singular, i.e. it cannot be inverted to compute the
As an alternative, many methods have been suggested to improve the estimation of the covariance matrix. All of these approaches rely on the concept of shrinkage. This is implicit in
Bayesianmethods, in penalized maximum likelihoodmethods, and explicit in the Stein-type shrinkage approach.
A simple version of a shrinkage estimator of the covariance matrix is constructed as follows. One considers a
convex combinationof the empirical estimator with some suitable chosen target, e.g., the diagonal matrix. Subsequently, the mixing parameteris selected to maximize the expected accuracy of the shrunken estimator. This can be done by cross-validation, or by using an analytic estimate of the shrinkage intensity. The resulting regularized estimator can be shown to outperform the maximum likelihood estimator for small samples. For large samples, the shrinkage intensity will reduce to zero, hence inthis case the shrinkage estimator will be identical to the empirical estimator. Apart from increased efficiency the shrinkage estimate has the additional advantage that it is always positive definite and well conditioned.
A review on this topic is given, e.g., in Schäfer and Strimmer 2005. [J. Schäfer and K. Strimmer (2005) " [http://www.bepress.com/sagmb/vol4/iss1/art32 A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics] ", Statistical Applications in Genetics and Molecular Biology: Vol. 4: No. 1, Article 32.]
A covariance shrinkage estimator is implemented in the R package [http://cran.r-project.org/web/packages/corpcor/index.html "corpcor"] .
Sample mean and sample covariance
Wikimedia Foundation. 2010.