Estimación de matrices de covarianza

Concepto de estadística

En estadística , a veces no se conoce la matriz de covarianza de una variable aleatoria multivariante, pero debe estimarse . La estimación de matrices de covarianza se ocupa entonces de la cuestión de cómo aproximar la matriz de covarianza real sobre la base de una muestra de la distribución multivariante . Los casos simples, donde las observaciones están completas, pueden abordarse utilizando la matriz de covarianza de muestra . La matriz de covarianza de muestra (SCM) es un estimador imparcial y eficiente de la matriz de covarianza si el espacio de matrices de covarianza se considera como un cono convexo extrínseco en R ^{p × p} ; sin embargo, medido utilizando la geometría intrínseca de matrices definidas positivas , el SCM es un estimador sesgado e ineficiente. ^[1] Además, si la variable aleatoria tiene una distribución normal , la matriz de covarianza de muestra tiene una distribución Wishart y una versión escalada ligeramente diferente de esta es la estimación de máxima verosimilitud . Los casos que involucran datos faltantes , heterocedasticidad o residuos autocorrelacionados requieren consideraciones más profundas. Otro problema es la robustez ante valores atípicos , a los que las matrices de covarianza de muestras son altamente sensibles. ^{[2] [3] [4]}

Los análisis estadísticos de datos multivariados suelen implicar estudios exploratorios de la forma en que las variables cambian entre sí, y esto puede ir seguido de modelos estadísticos explícitos que involucran la matriz de covarianza de las variables. Por lo tanto, la estimación de matrices de covarianza directamente a partir de datos observacionales cumple dos funciones:

• proporcionar estimaciones iniciales que puedan utilizarse para estudiar las interrelaciones;
• Proporcionar estimaciones de muestra que puedan utilizarse para la verificación del modelo.

Las estimaciones de matrices de covarianza son necesarias en las etapas iniciales del análisis de componentes principales y del análisis factorial , y también se utilizan en versiones del análisis de regresión que tratan las variables dependientes en un conjunto de datos, junto con la variable independiente , como el resultado de una muestra aleatoria.

Estimación en un contexto general

Dada una muestra que consta de n observaciones independientes x ₁ ,..., x _n de un vector aleatorio p -dimensional X ∈ R ^{p ×1} (un vector columna p ×1), un estimador insesgado de la matriz de covarianza ( p × p )

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

es la matriz de covarianza de muestra

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

donde es la i -ésima observación del vector aleatorio p -dimensional, y el vector ${\displaystyle x_{i}}$

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

es la media de la muestra . Esto es cierto independientemente de la distribución de la variable aleatoria X , siempre que existan las medias y covarianzas teóricas. La razón para el factor n  − 1 en lugar de n es esencialmente la misma que la razón para que el mismo factor aparezca en estimaciones no sesgadas de varianzas y covarianzas de la muestra , que se relaciona con el hecho de que la media no se conoce y se reemplaza por la media de la muestra (ver la corrección de Bessel ).

En los casos en que se sabe que la distribución de la variable aleatoria X se encuentra dentro de una determinada familia de distribuciones, se pueden derivar otras estimaciones sobre la base de ese supuesto. Un caso bien conocido es cuando la variable aleatoria X se distribuye normalmente : en este caso, el estimador de máxima verosimilitud de la matriz de covarianza es ligeramente diferente de la estimación insesgada y se da por

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

A continuación se ofrece una derivación de este resultado. Claramente, la diferencia entre el estimador insesgado y el estimador de máxima verosimilitud disminuye para valores grandes de n .

En el caso general, la estimación no sesgada de la matriz de covarianza proporciona una estimación aceptable cuando los vectores de datos en el conjunto de datos observados están todos completos: es decir, no contienen elementos faltantes . Un enfoque para estimar la matriz de covarianza es tratar la estimación de cada varianza o covarianza por pares por separado y utilizar todas las observaciones para las que ambas variables tienen valores válidos. Suponiendo que los datos faltantes faltan al azar, esto da como resultado una estimación para la matriz de covarianza que es no sesgada. Sin embargo, para muchas aplicaciones, esta estimación puede no ser aceptable porque no se garantiza que la matriz de covarianza estimada sea semidefinida positiva. Esto podría llevar a correlaciones estimadas que tengan valores absolutos que sean mayores que uno y/o una matriz de covarianza no invertible.

Al estimar la covarianza cruzada de un par de señales que son estacionarias en sentido amplio , las muestras faltantes no necesitan ser aleatorias (por ejemplo, el submuestreo por un factor arbitrario es válido). ^{[ cita requerida ]}

Estimación de máxima verosimilitud para la distribución normal multivariante

Un vector aleatorio X ∈ R ^p (un "vector columna" p ×1) tiene una distribución normal multivariada con una matriz de covarianza no singular Σ precisamente si Σ ∈ R ^{p × p} es una matriz definida positiva y la función de densidad de probabilidad de X es

${\displaystyle 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)}$

donde μ ∈ R ^{p ×1} es el valor esperado de X . La matriz de covarianza Σ es el análogo multidimensional de lo que en una dimensión sería la varianza , y

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

normaliza la densidad para que se integre a 1. ${\displaystyle f(x)}$

Supongamos ahora que X ₁ , ..., X _n son muestras independientes e idénticamente distribuidas de la distribución anterior. Con base en los valores observados x ₁ , ..., x _n de esta muestra , deseamos estimar Σ.

Pinitos

La función de verosimilitud es:

${\displaystyle {\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)}$

Se demuestra con bastante facilidad que la estimación de máxima verosimilitud del vector medio μ es el vector de " media de la muestra ":

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

Consulte la sección sobre estimación en el artículo sobre la distribución normal para obtener más detalles; el proceso aquí es similar.

Dado que la estimación no depende de Σ, podemos simplemente sustituirla por μ en la función de verosimilitud , obteniendo ${\displaystyle {\bar {x}}}$

${\displaystyle {\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),}$

y luego buscar el valor de Σ que maximice la verosimilitud de los datos (en la práctica es más fácil trabajar con log  ). ${\displaystyle {\mathcal {L}}}$

La traza de una matriz 1 × 1

Ahora llegamos al primer paso sorprendente: considerar el escalar como la traza de una matriz 1×1. Esto hace posible utilizar la identidad tr( AB ) = tr( BA ) siempre que A y B sean matrices con una forma tal que existan ambos productos. Obtenemos ${\displaystyle (x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})}$

${\displaystyle {\begin{aligned}{\mathcal {L}}({\overline {x}},\Sigma )&\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\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}}}$

dónde

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

${\displaystyle S}$A veces se denomina matriz de dispersión y es definida positiva si existe un subconjunto de los datos que consiste en observaciones afínmente independientes (lo cual asumiremos). ${\displaystyle p}$

Utilizando el teorema espectral

Del teorema espectral del álgebra lineal se deduce que una matriz simétrica positiva definida S tiene una única raíz cuadrada simétrica positiva definida S ^1/2 . Podemos utilizar nuevamente la "propiedad cíclica" de la traza para escribir

${\displaystyle \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).}$

Sea B = S ^1/2 Σ ⁻¹ S ^1/2 . Entonces la expresión anterior se convierte en

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

La matriz definida positiva B se puede diagonalizar y luego resolver el problema de encontrar el valor de B que maximice

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

Dado que la traza de una matriz cuadrada es igual a la suma de los valores propios ( "traza y valores propios" ), la ecuación se reduce al problema de encontrar los valores propios λ ₁ , ..., λ _p que maximizan

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

Este es solo un problema de cálculo y obtenemos λ _i = n para todo i. Por lo tanto, supongamos que Q es la matriz de vectores propios, entonces

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

es decir, n veces la matriz identidad p × p .

Pasos finales

Por fin lo conseguimos

${\displaystyle \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}},}$

es decir, la "matriz de covarianza de muestra" p × p

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

es el estimador de máxima verosimilitud de la "matriz de covarianza de la población" Σ. En este punto, estamos utilizando una X mayúscula en lugar de una x minúscula porque la consideramos "como un estimador en lugar de como una estimación", es decir, como algo aleatorio cuya distribución de probabilidad nos resultaría útil conocer. Se puede demostrar que la matriz aleatoria S tiene una distribución Wishart con n − 1 grados de libertad. ^[5] Es decir:

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

Derivación alternativa

Se puede realizar una derivación alternativa del estimador de máxima verosimilitud mediante fórmulas de cálculo matricial (véase también diferencial de un determinante y diferencial de la matriz inversa ). También se verifica el hecho mencionado anteriormente sobre la estimación de máxima verosimilitud de la media. Reescriba la verosimilitud en forma logarítmica utilizando el truco de la traza:

${\displaystyle \ln {\mathcal {L}}(\mu ,\Sigma )=\operatorname {constant} -{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].}$

El diferencial de esta verosimilitud logarítmica es

${\displaystyle 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].}$

Naturalmente, se descompone en la parte relacionada con la estimación de la media y la parte relacionada con la estimación de la varianza. La condición de primer orden para el máximo, , se satisface cuando los términos que se multiplican por y son idénticos a cero. Suponiendo que (la estimación de máxima verosimilitud de) no es singular, la condición de primer orden para la estimación del vector de media es ${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$ ${\displaystyle d\mu }$ ${\displaystyle d\Sigma }$ ${\displaystyle \Sigma }$

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

lo que conduce al estimador de máxima verosimilitud

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

Esto nos permite simplificar

${\displaystyle \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}$

como se definió anteriormente. Entonces los términos que involucran pueden combinarse como ${\displaystyle d\Sigma }$ ${\displaystyle d\ln L}$

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

La condición de primer orden se cumplirá cuando el término entre corchetes sea cero (valor matricial). Al multiplicar previamente este último por y dividir por se obtiene ${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$ ${\displaystyle \Sigma }$ ${\displaystyle n}$

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

lo cual por supuesto coincide con la derivación canónica dada anteriormente.

Dwyer ^[6] señala que la descomposición en dos términos como aparece arriba es "innecesaria" y deriva el estimador en dos líneas de trabajo. Obsérvese que puede no ser trivial demostrar que dicho estimador derivado es el único maximizador global de la función de verosimilitud.

Estimación de la matriz de covarianza intrínseca

Expectativa intrínseca

Dada una muestra de n observaciones independientes x ₁ ,..., x _n de una variable aleatoria gaussiana p -dimensional de media cero X con covarianza R , el estimador de máxima verosimilitud de R viene dado por

${\displaystyle {\hat {\mathbf {R} }}={1 \over n}\sum _{i=1}^{n}x_{i}x_{i}^{\mathrm {T} }.}$

El parámetro pertenece al conjunto de matrices definidas positivas , que es una variedad de Riemann , no un espacio vectorial , por lo tanto, las nociones habituales de espacio vectorial de expectativa , es decir, " ", y sesgo del estimador deben generalizarse a las variedades para dar sentido al problema de la estimación de la matriz de covarianza. Esto se puede hacer definiendo la expectativa de un estimador con valor de variedad con respecto al punto con valor de variedad como ${\displaystyle R}$ ${\displaystyle \mathrm {E} [{\hat {\mathbf {R} }}]}$ ${\displaystyle {\hat {\mathbf {R} }}}$ ${\displaystyle R}$

${\displaystyle \mathrm {E} _{\mathbf {R} }[{\hat {\mathbf {R} }}]\ {\stackrel {\mathrm {def} }{=}}\ \exp _{\mathbf {R} }\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

dónde

${\displaystyle \exp _{\mathbf {R} }({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\exp \left(\mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

${\displaystyle \exp _{\mathbf {R} }^{-1}({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\left(\log \mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

son el mapa exponencial y el mapa exponencial inverso, respectivamente, "exp" y "log" denotan la matriz exponencial ordinaria y el logaritmo matricial , y E[·] es el operador de expectativa ordinaria definido en un espacio vectorial, en este caso el espacio tangente de la variedad. ^[1]

Sesgo de la matriz de covarianza de la muestra

El campo vectorial de sesgo intrínseco del estimador SCM se define como ${\displaystyle {\hat {\mathbf {R} }}}$

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=\exp _{\mathbf {R} }^{-1}\mathrm {E} _{\mathbf {R} }\left[{\hat {\mathbf {R} }}\right]=\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

El sesgo intrínseco del estimador viene dado por . ${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})}$

Para variables aleatorias gaussianas complejas , se puede demostrar ^[1] que este campo vectorial de sesgo es igual a

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=-\beta (p,n)\mathbf {R} }$

dónde

${\displaystyle \beta (p,n)={\frac {1}{p}}\left(p\log n+p-\psi (n-p+1)+(n-p+1)\psi (n-p+2)+\psi (n+1)-(n+1)\psi (n+2)\right)}$

y ψ(·) es la función digamma . El sesgo intrínseco de la matriz de covarianza de la muestra es igual a

${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})=e^{-\beta (p,n)}\mathbf {R} }$

y el SCM es asintóticamente imparcial cuando n → ∞.

De manera similar, la ineficiencia intrínseca de la matriz de covarianza de la muestra depende de la curvatura riemanniana del espacio de matrices definidas positivas.

Estimación de la contracción

Si el tamaño de la muestra n es pequeño y el número de variables consideradas p es grande, los estimadores empíricos anteriores de covarianza y correlación son muy inestables. En concreto, es posible proporcionar estimadores que mejoren considerablemente la estimación de máxima verosimilitud en términos de error cuadrático medio. Además, para n  <  p (el número de observaciones es menor que el número de variables aleatorias), la estimación empírica de la matriz de covarianza se vuelve singular , es decir, no se puede invertir para calcular la matriz de precisión .

Como alternativa, se han sugerido muchos métodos para mejorar la estimación de la matriz de covarianza. Todos estos enfoques se basan en el concepto de contracción, que está implícito en los métodos bayesianos y en los métodos de máxima verosimilitud penalizada , y es explícito en el enfoque de contracción de tipo Stein .

Una versión simple de un estimador de contracción de la matriz de covarianza está representada por el estimador de contracción de Ledoit-Wolf. ^{[7] [8] [9] [10]} Se considera una combinación convexa del estimador empírico ( ) con algún objetivo elegido adecuado ( ), por ejemplo, la matriz diagonal. Posteriormente, se selecciona el parámetro de mezcla ( ) para maximizar la precisión esperada del estimador contraído. Esto se puede hacer mediante validación cruzada o utilizando una estimación analítica de la intensidad de la contracción. Se puede demostrar que el estimador regularizado resultante ( ) supera al estimador de máxima verosimilitud para muestras pequeñas. Para muestras grandes, la intensidad de la contracción se reducirá a cero, por lo tanto, en este caso, el estimador de contracción será idéntico al estimador empírico. Aparte de una mayor eficiencia, la estimación de la contracción tiene la ventaja adicional de que siempre es positiva definida y bien condicionada. ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle \delta }$ ${\displaystyle \delta A+(1-\delta )B}$

Se han propuesto varios objetivos de contracción:

1. la matriz identidad , escalada por la varianza promedio de la muestra ;
2. el modelo de índice único ;
3. el modelo de correlación constante, donde se conservan las varianzas de la muestra, pero se supone que todos los coeficientes de correlación por pares son iguales entre sí;
4. la matriz de dos parámetros, donde todas las varianzas son idénticas y todas las covarianzas son idénticas entre sí (aunque no idénticas a las varianzas);
5. la matriz diagonal que contiene las varianzas de muestra en la diagonal y ceros en el resto;
6. La matriz identidad . ^[8]

El estimador de contracción se puede generalizar a un estimador de contracción de múltiples objetivos que utiliza varios objetivos simultáneamente. ^[11] El software para calcular un estimador de contracción de covarianza está disponible en R (paquetes corpcor ^[12] y ShrinkCovMat ^[13] ), en Python ( biblioteca scikit-learn [1]) y en MATLAB . ^[14]

Véase también

• Propagación de la incertidumbre
• Media de la muestra y covarianza de la muestra
• Componentes de la varianza

Referencias

1. Smith, Steven Thomas (mayo de 2005). "Covarianza, subespacio y límites intrínsecos de Cramér–Rao". IEEE Trans. Signal Process . 53 (5): 1610–1630. doi :10.1109/TSP.2005.845428. S2CID  2751194.
2. Robust Statistics , Peter J. Huber , Wiley, 1981 (publicado nuevamente en edición de bolsillo, 2004)
3. "Estadísticas aplicadas modernas con S", William N. Venables, Brian D. Ripley , Springer, 2002, ISBN 0-387-95457-0 , ISBN 978-0-387-95457-8 , página 336  
4. Devlin, Susan J. ; Gnanadesikan, R.; Kettenring, JR (1975). "Estimación robusta y detección de valores atípicos con coeficientes de correlación". Biometrika . 62 (3): 531–545. doi :10.1093/biomet/62.3.531.
5. KV Mardia , JT Kent y JM Bibby (1979) Análisis multivariado , Academic Press .
6. Dwyer, Paul S. (junio de 1967). "Algunas aplicaciones de las derivadas matriciales en el análisis multivariante". Revista de la Asociación Estadounidense de Estadística . 62 (318): 607–625. doi :10.2307/2283988. JSTOR  2283988.
7. O. Ledoit y M. Wolf (2004a) "Un estimador bien condicionado para matrices de covarianza de gran dimensión Archivado el 5 de diciembre de 2014 en Wayback Machine " Journal of Multivariate Analysis 88 (2): 365—411.
8. A. Touloumis (2015) "Estimadores de matriz de covarianza de contracción de tipo Stein no paramétricos en entornos de alta dimensión" Computational Statistics & Data Analysis 83 : 251—261.
9. O. Ledoit y M. Wolf (2003) "Estimación mejorada de la matriz de covarianza de los retornos de acciones con una aplicación a la selección de carteras Archivado el 5 de diciembre de 2014 en Wayback Machine " Journal of Empirical Finance 10 (5): 603—621.
10. O. Ledoit y M. Wolf (2004b) "Cariño, encogí la matriz de covarianza de muestra Archivado el 5 de diciembre de 2014 en Wayback Machine ." The Journal of Portfolio Management 30 (4): 110—119.
11. T. Lancewicki y M. Aladjem (2014) "Estimación de contracción de múltiples objetivos para matrices de covarianza", IEEE Transactions on Signal Processing , Volumen: 62, Número 24, páginas: 6380-6390.
12. corpcor: Estimación eficiente de la covarianza y la correlación (parcial), CRAN, 16 de septiembre de 2021{{citation}}: CS1 maint: location missing publisher (link)
13. ShrinkCovMat: Estimadores de la matriz de covarianza de contracción, CRAN, 30 de julio de 2019{{citation}}: CS1 maint: location missing publisher (link)
14. Código MATLAB para objetivos de contracción: identidad escalada, modelo de índice único, modelo de correlación constante, matriz de dos parámetros y matriz diagonal.