Benutzer:JonskiC/Schätzung der Kovarianzmatrix

Dieser Artikel ist im Entstehen und noch nicht Bestandteil der freien Enzyklopädie Wikipedia.

Solltest du über eine Suchmaschine darauf gestoVorlage:SSen sein, bedenke, dass der Text noch unvollständig sein und Fehler oder ungeprüfte Aussagen enthalten kann. Wenn du Fragen zum Thema hast, nimm Kontakt mit dem Autor JonskiC auf.

In statistics, sometimes the covariance matrix of a multivariate random variable is 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 (SCM) is an unbiased and efficient estimator of the covariance matrix if the space of covariance matrices is viewed as an extrinsic convex cone in R^p×p; however, measured using the intrinsic geometry of positive-definite matrices, the SCM is a biased and inefficient estimator. In addition, if the random variable has normal distribution, the sample covariance matrix has Wishart distribution and a slightly differently scaled version of it is the maximum likelihood estimate. Cases involving missing data require deeper considerations. Another issue is the robustness to outliers, to which sample covariance matrices are highly sensitive.

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 initial stages of principal component analysis and factor analysis, and are also involved in versions of regression analysis that treat the dependent variables in a data-set, jointly with the independent variable as the outcome of a random sample.

Estimation in a general context

Given a sample consisting of n independent observations x₁,..., x_n of a p-dimensional random vector x ∈ R^p×1 (a p×1 column-vector), an unbiased estimator of the (p×p) covariance matrix

$\operatorname {\Sigma } =\operatorname {E} \left[\left(X-\operatorname {E} [X]\right)\left(X-\operatorname {E} [X]\right)^{\mathrm {T} }\right]$

is the sample covariance matrix

$\mathbf {Q} ={1 \over {n-1}}\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} },$

where $x_{i}$ is the i-th observation of the p-dimensional random vector, and the vector

${\overline {x}}={1 \over {n}}\sum _{i=1}^{n}x_{i}$

is the sample mean. This is true regardless of the distribution of the random variable 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 sample 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 estimates 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 estimator of the covariance matrix is slightly different from the unbiased estimate, and is given by

$\mathbf {Q_{n}} ={1 \over n}\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} }.$

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 random this 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, and/or a non-invertible covariance matrix.

When estimating the cross-covariance of a pair of signals that are wide-sense stationary, missing samples do not need be random (e.g., sub-sampling by an arbitrary factor is valid).Vorlage:Citation needed

Maximum-likelihood estimation for the multivariate normal distribution

A random vector X ∈ R^p (a p×1 "column vector") has a multivariate normal distribution with a nonsingular covariance matrix Σ precisely if Σ ∈ R^{p × p} is a positive-definite matrix and the probability density function of X is

$f(x)=(2\pi )^{-{\frac {p}{2}}}\,\det(\Sigma )^{-{\frac {1}{2}}}\exp \left(-{1 \over 2}(x-\mu )^{\mathrm {T} }\Sigma ^{-1}(x-\mu )\right)$

where μ ∈ R^p×1 is the expected value of X. The covariance matrix Σ is the multidimensional analog of what in one dimension would be the variance, and

$(2\pi )^{-{\frac {p}{2}}}\det(\Sigma )^{-{\frac {1}{2}}}$

normalizes the density $f(x)$ so that it integrates to 1.

Suppose now that X₁, ..., X_n are independent and identically distributed samples from the distribution above. Based on the observed values x₁, ..., x_n of this sample, we wish to estimate Σ.

First steps

The likelihood function is:

${\mathcal {L}}(\mu ,\Sigma )=(2\pi )^{-{\frac {np}{2}}}\,\prod _{i=1}^{n}\det(\Sigma )^{-{\frac {1}{2}}}\exp \left(-{\frac {1}{2}}(x_{i}-\mu )^{\mathrm {T} }\Sigma ^{-1}(x_{i}-\mu )\right)$

It is fairly readily shown that the maximum-likelihood estimate of the mean vector μ is the "sample mean" vector:

${\overline {x}}={\frac {x_{1}+\cdots +x_{n}}{n}}.$

See the section on estimation in the article on the normal distribution for details; the process here is similar.

Since the estimate ${\bar {x}}$ does not depend on Σ, we can just substitute it for μ in the likelihood function, getting

${\mathcal {L}}({\overline {x}},\Sigma )\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}(x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})\right),$

and then seek the value of Σ that maximizes the likelihood of the data (in practice it is easier to work with log ${\mathcal {L}}$).

The trace of a 1 × 1 matrix

Now we come to the first surprising step: regard the scalar $(x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})$ as 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

${\begin{aligned}{\mathcal {L}}({\overline {x}},\Sigma )&\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\operatorname {tr} \left(\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\left(x_{i}-{\overline {x}}\right)\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\operatorname {tr} \left(\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S\Sigma ^{-1}\right)\right)\end{aligned}}$

where

$S=\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} }\in \mathbf {R} ^{p\times p}.$

$S$ is sometimes called the scatter matrix, and is positive definite if there exists a subset of the data consisting of $p$ affinely independent observations (which we will assume).

Using the spectral theorem

It follows from the spectral theorem of linear algebra that 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

$\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S^{\frac {1}{2}}\Sigma ^{-1}S^{\frac {1}{2}}\right)\right).$

Let B = S^1/2 Σ ⁻¹ S^1/2. Then the expression above becomes

$\det(S)^{-{\frac {n}{2}}}\det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right).$

The positive-definite matrix B can be diagonalized, and then the problem of finding the value of B that maximizes

$\det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right)$

Since the trace of a square matrix equals the sum of eigen-values ("trace and eigenvalues"), the equation reduces to the problem of finding the eigen values λ₁, ..., λ_p that maximize

$\lambda _{i}^{\frac {n}{2}}\exp \left(-{\frac {\lambda _{i}}{2}}\right).$

This is just a calculus problem and we get λ_i = n for all i. Thus, assume Q is the matrix of eigen vectors, then

$B=Q(nI_{p})Q^{-1}=nI_{p}$

i.e., n times the p×p identity matrix.

Concluding steps

Finally we get

$\Sigma =S^{\frac {1}{2}}B^{-1}S^{\frac {1}{2}}=S^{\frac {1}{2}}\left({\frac {1}{n}}I_{p}\right)S^{\frac {1}{2}}={\frac {S}{n}},$

i.e., the p×p "sample covariance matrix"

${S \over n}={1 \over n}\sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }$

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 distribution with n − 1 degrees of freedom. That is:

$\sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }\sim W_{p}(\Sigma ,n-1).$

Alternative derivation

An alternative derivation of the maximum likelihood estimator can be performed via matrix calculus formulae (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:

$\ln {\mathcal {L}}(\mu ,\Sigma )=\operatorname {const} -{n \over 2}\ln \det(\Sigma )-{1 \over 2}\operatorname {tr} \left[\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }\right].$

The differential of this log-likelihood is

$d\ln {\mathcal {L}}(\mu ,\Sigma )=-{\frac {n}{2}}\operatorname {tr} \left[\Sigma ^{-1}\left\{d\Sigma \right\}\right]-{1 \over 2}\operatorname {tr} \left[-\Sigma ^{-1}\{d\Sigma \}\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }-2\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )\{d\mu \}^{\mathrm {T} }\right].$

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, $d\ln {\mathcal {L}}(\mu ,\Sigma )=0$, is satisfied when the terms multiplying $d\mu$ and $d\Sigma$ are identically zero. Assuming (the maximum likelihood estimate of) $\Sigma$ is non-singular, the first order condition for the estimate of the mean vector is

$\sum _{i=1}^{n}(x_{i}-\mu )=0,$

which leads to the maximum likelihood estimator

${\widehat {\mu }}={\bar {X}}={1 \over n}\sum _{i=1}^{n}X_{i}.$

This lets us simplify

$\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }=\sum _{i=1}^{n}(x_{i}-{\bar {x}})(x_{i}-{\bar {x}})^{\mathrm {T} }=S$

as defined above. Then the terms involving $d\Sigma$ in $d\ln L$ can be combined as

$-{1 \over 2}\operatorname {tr} \left(\Sigma ^{-1}\left\{d\Sigma \right\}\left[nI_{p}-\Sigma ^{-1}S\right]\right).$

The first order condition $d\ln {\mathcal {L}}(\mu ,\Sigma )=0$ will hold when the term in the square bracket is (matrix-valued) zero. Pre-multiplying the latter by $\Sigma$ and dividing by $n$ gives

${\widehat {\Sigma }}={1 \over n}S,$

which of course coincides with the canonical derivation given earlier.

Dwyer points out that decomposition into two terms such as appears above is "unnecessary" and derives the estimator in two lines of working. Note that it may be not trivial to show that such derived estimator is the unique global maximizer for likelihood function.