Estimativa de matrizes de covariância - Estimation of covariance matrices

Em estatística , às vezes a matriz de covariância de uma variável aleatória multivariada não é conhecida, mas deve ser estimada . A estimativa de matrizes de covariância trata então da questão de como aproximar a matriz de covariância real com base em uma amostra da distribuição multivariada . Casos simples, onde as observações são completas, podem ser tratados usando a matriz de covariância de amostra . A matriz de covariância de amostra (SCM) é um estimador não enviesado e eficiente da matriz de covariância se o espaço das matrizes de covariância for visto como um cone convexo extrínseco em R ^{p × p} ; no entanto, medido usando a geometria intrínseca de matrizes definidas positivas , o SCM é um estimador enviesado e ineficiente. Além disso, se a variável aleatória tem distribuição normal , a matriz de covariância da amostra tem distribuição de Wishart e uma versão em escala ligeiramente diferente dela é a estimativa de máxima verossimilhança . Casos envolvendo dados perdidos requerem considerações mais profundas. Outra questão é a robustez para outliers , para os quais as matrizes de covariância de amostra são altamente sensíveis.

As análises estatísticas de dados multivariados frequentemente envolvem estudos exploratórios da maneira como as variáveis ​​mudam em relação umas às outras e isso pode ser seguido por modelos estatísticos explícitos envolvendo a matriz de covariância das variáveis. Assim, a estimativa de matrizes de covariância diretamente a partir de dados observacionais desempenha dois papéis:

• fornecer estimativas iniciais que podem ser usadas para estudar as inter-relações;
• para fornecer estimativas de amostra que podem ser usadas para verificação de modelo.

Estimativas de matrizes de covariância são necessárias nos estágios iniciais da análise de componentes principais e análise fatorial , e também estão envolvidas em versões de análise de regressão que tratam as variáveis ​​dependentes em um conjunto de dados, juntamente com a variável independente como o resultado de uma amostra aleatória .

Estimativa em um contexto geral

Dada uma amostra que consiste em n observações independentes x ₁ , ..., x _n de um vetor aleatório p- dimensional X ∈ R ^{p × 1} (a p × 1 vetor coluna), um estimador imparcial do ( p × p ) matriz de covariância

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

é a matriz de covariância de amostra

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

onde é a i -ésima observação do vetor aleatório p- dimensional, e o vetor ${\ displaystyle x_ {i}}$

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

é a média da amostra . Isso é verdade independentemente da distribuição da variável aleatória X , contanto, é claro, que os meios teóricos e as covariâncias existam. A razão para o fator n  - 1 em vez de n é essencialmente a mesma que a razão para o mesmo fator aparecer em estimativas não enviesadas de variâncias e covariâncias amostrais , que se relacionam ao fato de que a média não é conhecida e é substituída pela amostra significa (ver a correção de Bessel ).

Nos casos em que se sabe que a distribuição da variável aleatória X está dentro de uma certa família de distribuições, outras estimativas podem ser derivadas com base nessa suposição. Uma instância bem conhecida é quando a variável aleatória X é normalmente distribuída : neste caso, o estimador de máxima verossimilhança da matriz de covariância é ligeiramente diferente da estimativa não enviesada e é dado por

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

Uma derivação desse resultado é fornecida abaixo. Claramente, a diferença entre o estimador não enviesado e o estimador de máxima verossimilhança diminui para n grande .

No caso geral, a estimativa não enviesada da matriz de covariância fornece uma estimativa aceitável quando os vetores de dados no conjunto de dados observado estão todos completos: ou seja, eles não contêm elementos ausentes . Uma abordagem para estimar a matriz de covariância é tratar a estimativa de cada variância ou covariância par a par separadamente e usar todas as observações para as quais ambas as variáveis ​​têm valores válidos. Assumindo que os dados faltantes estão faltando aleatoriamente, isso resulta em uma estimativa para a matriz de covariância que é imparcial. No entanto, para muitas aplicações, esta estimativa pode não ser aceitável porque a matriz de covariância estimada não é garantida como semidefinida positiva. Isso pode levar a correlações estimadas com valores absolutos maiores do que um e / ou uma matriz de covariância não invertível.

Ao estimar a covariância cruzada de um par de sinais estacionários de sentido amplo , as amostras ausentes não precisam ser aleatórias (por exemplo, a subamostragem por um fator arbitrário é válida).

Estimativa de máxima verossimilhança para a distribuição normal multivariada

Um vetor aleatório X ∈ R ^p (a p × 1 "vetor coluna") tem uma distribuição normal multivariada com uma matriz de covariância não singular Σ precisamente se Σ ∈ R ^{p × p} é uma matriz definida positiva e a função de densidade de probabilidade de X é

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

onde u ∈ R ^{p × 1} é o valor esperado de x . A matriz de covariância Σ é o análogo multidimensional do que em uma dimensão seria a variância , e

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

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

Suponha agora que X ₁ , ..., X _n são amostras independentes e distribuídas de forma idêntica da distribuição acima. Com base nos valores observados x ₁ , ..., x _n desta amostra , desejamos estimar Σ.

Primeiros passos

A função de verossimilhança é:

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

É bastante facilmente mostrado que a estimativa de máxima verossimilhança do vetor médio μ é o vetor " média da amostra ":

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

Consulte a seção sobre estimativa no artigo sobre distribuição normal para obter detalhes; o processo aqui é semelhante.

Uma vez que a estimativa não depende de Σ, podemos apenas substituí-la por μ na função de verossimilhança , obtendo ${\ 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),}$

e então busque o valor de Σ que maximize a probabilidade dos dados (na prática é mais fácil trabalhar com log  ). ${\ displaystyle {\ mathcal {L}}}$

O traço de uma matriz 1 × 1

Agora chegamos ao primeiro passo surpreendente: considere o escalar como o traço de uma matriz 1 × 1. Isso torna possível usar a identidade tr ( AB ) = tr ( BA ) sempre que A e B são matrizes formadas de modo que ambos os produtos existam. Nós temos ${\ displaystyle (x_ {i} - {\ overline {x}}) ^ {\ mathrm {T}} \ Sigma ^ {- 1} (x_ {i} - {\ overline {x}})}$

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

Onde

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

${\ displaystyle S}$às vezes é chamada de matriz de dispersão e é definida positivamente se existir um subconjunto dos dados consistindo em observações afinamente independentes (o que iremos supor). ${\ displaystyle p}$

Usando o teorema espectral

Segue-se do teorema espectral da álgebra linear que uma matriz simétrica positiva-definida S tem uma única raiz quadrada simétrica positiva-definida S ^1/2 . Podemos novamente usar a "propriedade cíclica" do traço para escrever

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

Seja B = S ^1/2 Σ ⁻¹ S ^1/2 . Então a expressão acima se torna

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

A matriz positiva-definida B pode ser diagonalizada, e então o problema de encontrar o valor de B que maximiza

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

Como o traço de uma matriz quadrada é igual à soma dos autovalores ( "traço e autovalores" ), a equação se reduz ao problema de encontrar os autovalores λ ₁ , ..., λ _p que maximizam

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

Este é apenas um problema de cálculo e obtemos λ _i = n para todo i. Assim, suponha que Q é a matriz de vetores próprios, então

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

ou seja, n vezes o p × p matriz identidade.

Passos finais

Finalmente conseguimos

${\ displaystyle \ Sigma = S ^ {\ frac {1} {2}} B ^ {- 1} S ^ {\ frac {1} {2}} = S ^ {\ frac {1} {2}} \ esquerda ({\ frac {1} {n}} I_ {p} \ direita) S ^ {\ frac {1} {2}} = {\ frac {S} {n}},}$

ou seja, o p × p "matriz de covariância de amostra"

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

é o estimador de máxima verossimilhança da "matriz de covariância populacional" Σ. Neste ponto, estamos usando um X maiúsculo em vez de um x minúsculo porque o estamos pensando "como um estimador e não como uma estimativa", isto é, como algo aleatório cuja distribuição de probabilidade poderíamos lucrar sabendo. Pode-se mostrar que a matriz aleatória S tem uma distribuição de Wishart com n - 1 graus de liberdade. Isso é:

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

Derivação alternativa

Uma derivação alternativa do estimador de máxima verossimilhança pode ser realizada por meio de fórmulas de cálculo de matriz (ver também diferencial de um determinante e diferencial da matriz inversa ). Verifica também o fato mencionado sobre a estimativa de máxima verossimilhança da média. Reescreva a probabilidade no formulário de registro usando o truque de rastreamento:

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

O diferencial dessa probabilidade logarítmica é

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

Naturalmente, ele se divide na parte relacionada à estimativa da média e na parte relacionada à estimativa da variância. A condição de primeira ordem para máximo ,, é satisfeita quando os termos se multiplicam e são idênticos a zero. Supondo que (a estimativa de máxima verossimilhança de) não seja singular, a condição de primeira ordem para a estimativa do vetor médio é ${\ 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,}$

o que leva ao estimador de máxima verossimilhança

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

Isso 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}$

conforme definido acima. Em seguida, os termos que envolvem em podem ser combinados 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).}$

A condição de primeira ordem será mantida quando o termo entre colchetes for (com valor de matriz) zero. A pré-multiplicação do último por e a divisão por dá ${\ displaystyle d \ ln {\ mathcal {L}} (\ mu, \ Sigma) = 0}$${\ displaystyle \ Sigma}$${\ displaystyle n}$

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

que naturalmente coincide com a derivação canônica dada anteriormente.

Dwyer aponta que a decomposição em dois termos, como aparece acima, é "desnecessária" e deriva o estimador em duas linhas de trabalho. Observe que pode não ser trivial mostrar que tal estimador derivado é o único maximizador global para a função de verossimilhança.

Estimativa da matriz de covariância intrínseca

Expectativa intrínseca

Dada uma amostra de n observações independentes x ₁ , ..., x _n de uma variável aleatória Gaussiana de média zero p- dimensional X com covariância R , o estimador de máxima verossimilhança de R é dado por

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

O parâmetro pertence ao conjunto de matrizes definidas-positivas , que é uma variedade Riemanniana , não um espaço vetorial , portanto, as noções usuais de espaço vetorial de expectativa , ou seja, " ", e o viés do estimador deve ser generalizado para variedades para dar sentido ao problema de estimação da matriz de covariância. Isso pode ser feito definindo a expectativa de um estimador de valor múltiplo em relação ao ponto de valor múltiplo 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]}$

Onde

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

são o mapa exponencial e o mapa exponencial inverso, respectivamente, "exp" e "log" denotam a matriz exponencial comum e o logaritmo da matriz , e E [·] é o operador de expectativa comum definido em um espaço vetorial, neste caso o espaço tangente de o múltiplo.

Viés da matriz de covariância da amostra

O campo vetorial de polarização intrínseca do estimador SCM é definido para ser ${\ 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} }}\certo]}$

A tendência intrínseca do estimador é então dada por . ${\ displaystyle \ exp _ {\ mathbf {R}} \ mathbf {B} ({\ hat {\ mathbf {R}}})}$

Para variáveis ​​aleatórias gaussianas complexas , este campo de vetor de polarização pode ser mostrado como igual

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

Onde

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

e ψ (·) é a função digamma . A tendência intrínseca da matriz de covariância da amostra é igual a

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

e o SCM é assintoticamente imparcial como n → ∞.

Da mesma forma, a ineficiência intrínseca da matriz de covariância da amostra depende da curvatura Riemanniana do espaço de matrizes definidas positivas.

Estimativa de encolhimento

Se o tamanho da amostra n for pequeno e o número de variáveis ​​consideradas p for grande, os estimadores empíricos de covariância e correlação acima são muito instáveis. Especificamente, é possível fornecer estimadores que melhoram consideravelmente na estimativa de máxima verossimilhança em termos de erro quadrático médio. Além disso, para n  <  p (o número de observações é menor que o número de variáveis ​​aleatórias), a estimativa empírica da matriz de covariância torna-se singular , ou seja, não pode ser invertida para calcular a matriz de precisão .

Como alternativa, vários métodos têm sido sugeridos para melhorar a estimativa da matriz de covariância. Todas essas abordagens baseiam-se no conceito de encolhimento. Isso está implícito nos métodos bayesianos e nos métodos de máxima verossimilhança penalizados e explícito na abordagem de encolhimento do tipo Stein .

Uma versão simples de um estimador de contração da matriz de covariância é representada pelo estimador de contração de Ledoit-Wolf. Considera-se uma combinação convexa do estimador empírico ( ) com algum alvo escolhido adequado ( ), por exemplo, a matriz diagonal. Subsequentemente, o parâmetro de mistura ( ) é selecionado para maximizar a precisão esperada do estimador encolhido. Isso pode ser feito por validação cruzada ou usando uma estimativa analítica da intensidade de encolhimento. O estimador regularizado resultante ( ) pode ser mostrado para superar o estimador de máxima verossimilhança para pequenas amostras. Para grandes amostras, a intensidade de contração será reduzida a zero, portanto, neste caso, o estimador de contração será idêntico ao estimador empírico. Além de aumentar a eficiência, a estimativa de encolhimento tem a vantagem adicional de ser sempre positiva e bem condicionada. ${\ displaystyle A}$${\ displaystyle B}$${\ displaystyle \ delta}$${\ displaystyle \ delta A + (1- \ delta) B}$

Vários alvos de redução foram propostos:

1. a matriz de identidade , escalonada pela variância média da amostra ;
2. o modelo de índice único ;
3. o modelo de correlação constante, em que as variâncias da amostra são preservadas, mas todos os coeficientes de correlação de pares são considerados iguais entre si;
4. a matriz de dois parâmetros, onde todas as variâncias são idênticas e todas as covariâncias são idênticas umas às outras (embora não idênticas às variâncias);
5. a matriz diagonal contendo as variações da amostra na diagonal e zeros em todas as outras partes;
6. a matriz de identidade .

O estimador de encolhimento pode ser generalizado para um estimador de encolhimento de múltiplos alvos que utiliza vários alvos simultaneamente. O software para calcular um estimador de redução de covariância está disponível em R (packages corpcor e ShrinkCovMat ), em Python (biblioteca scikit-learn ) e em MATLAB .

Veja também

• Propagação de incerteza
• Média da amostra e covariância da amostra
• Componentes de variância

Referências