Como criar uma matriz de covariância arbitrária

21

Por exemplo, em R, a MASS::mvrnorm()função é útil para gerar dados para demonstrar várias coisas nas estatísticas. É necessário um Sigmaargumento obrigatório, que é uma matriz simétrica que especifica a matriz de covariância das variáveis. Como eu criaria uma matriz n × simétrica n×ncom entradas arbitrárias?

rsl
fonte
3
Eu acho que essa pergunta se beneficiaria de ser editada para focar em "como criar uma matriz de covariância arbitrária" e menos no aspecto da codificação. Certamente há uma questão estatística subjacente aqui, como demonstrado pela resposta.
Silverfish 31/05

Respostas:

22

Criar um matriz A com valores arbitráriosn×nA

e então use como sua matriz de covariância. Σ=ATA

Por exemplo

n <- 4  
A <- matrix(runif(n^2)*2-1, ncol=n) 
Sigma <- t(A) %*% A
Henry
fonte
Da mesma forma Sigma <- A + t(A),.
rsl 31/05
6
@MoazzemHossen: Sua sugestão irá produzir uma matriz simétrica, mas pode não ser sempre semidefinite positivo (por exemplo, a sua sugestão poderia produzir uma matriz com valores próprios negativos) e por isso não pode ser apropriado como uma matriz de covariância
Henry
Sim, notei que R retorna erro no caso de minha maneira sugerida produzir matriz inadequada.
rsl
4
Observe que, se você preferir uma matriz de correlação para melhor interpretabilidade, existe a função ? Cov2cor , que pode ser aplicada posteriormente.
gung - Restabelece Monica
1
@ B11b: Você precisa que sua matriz de covariância seja semi-definida positiva. Que iria colocar alguns limites sobre os valores de covariância, queridos não totalmente óbvias quando n>2
Henry
24

Eu gosto de ter controle sobre os objetos que eu crio, mesmo quando eles podem ser arbitrários.

Considere-se, então, que todos os possíveis matrizes de covariâncian×n pode ser expressa na formaΣ

Σ=P Diagonal(σ1,σ2,,σn) P

onde é uma matriz ortogonal e σ 1σ 2σ nP .σ1σ2σn0

Geometricamente, isso descreve uma estrutura de covariância com uma variedade de componentes principais de tamanhos . Estes componentes apontam em direcções das linhas de P . Veja as figuras em Compreendendo a análise de componentes principais, vetores próprios e valores próprios para exemplos com n = 3 . Configurando o σ iσiPn=3σi definirá as magnitudes das covariâncias e seus tamanhos relativos, determinando assim qualquer forma elipsoidal desejada. As linhas de orientam os eixos da forma como você preferir.P

Um benefício algébrico e computacional dessa abordagem é que, quando , Σσn>0Σ é prontamente invertido (que é uma operação comum em matrizes de covariância):

Σ1=P Diagonal(1/σ1,1/σ2,,1/σn) P.

Não se preocupe com as direções, mas apenas com os intervalos de tamanhos dos ? Tudo bem: você pode gerar facilmente uma matriz ortogonal aleatória. Apenas envolva n 2 iid valores normais padrão em uma matriz quadrada e ortogonalize-a. Certamente funcionará (desde que nσin2n não seja enorme). A decomposição QR fará isso, como neste código

n <- 5
p <- qr.Q(qr(matrix(rnorm(n^2), n)))

Isso funciona porque o n distribuição multinormal variável assim gerada é "elíptica": é invariável sob todas as rotações e reflexões (através da origem). Assim, todas as matrizes ortogonais são geradas uniformemente, conforme discutido em Como gerar pontos uniformemente distribuídos na superfície da esfera unitária 3-d? .

Uma maneira rápida de obter de P eo σ i , uma vez que você tenha especificado ou criou, usos e exploits 's re-uso de matrizes em operações aritméticas, como neste exemplo, com σ = ( σ 1 , ... , σ 5 ) = ( 5 , 4 , 3 , 2 , 1 ) :ΣPσicrossprodRσ=(σ1,,σ5)=(5,4,3,2,1)

Sigma <- crossprod(p, p*(5:1))

Como verificação, a decomposição do Valor Singular deve retornar e P . Você pode inspecioná-lo com o comandoσP

svd(Sigma)

O inverso Sigma, é claro, é obtido apenas mudando a multiplicação por σ em uma divisão:

Tau <- crossprod(p, p/(5:1))

Você pode verificar isso através da visualização zapsmall(Sigma %*% Tau), que deve ser o matriz identidade. Um inverso generalizado (essencial para os cálculos de regressão) é obtida pela substituição de qualquer σ i0 por 1 / σ i , exactamente como acima, mas mantendo os zeros entre o σ i como eram.n×nσi01/σiσi

whuber
fonte
P
1
Vale a pena mencionar que os valores singulares svd(Sigma)serão reordenados - isso me confundiu por um minuto.
precisa saber é
1

Você pode simular matrizes definidas positivas aleatórias a partir da distribuição Wishart usando a função "rWishart" do pacote amplamente utilizado "stats".

n <- 4
rWishart(1,n,diag(n))
Carlos Llosa
fonte
1

Existe um pacote especificamente para isso, clusterGeneration (escrito entre outros por Harry Joe, um grande nome nesse campo).

Existem duas funções principais:

  • genPositiveDefMat gerar uma matriz de covariância, 4 métodos diferentes
  • rcorrmatrix : gerar uma matriz de correlação

Exemplo rápido:

library(clusterGeneration)
#> Loading required package: MASS
genPositiveDefMat("unifcorrmat",dim=3)
#> $egvalues
#> [1] 15.408962  5.673916  1.228842
#> 
#> $Sigma
#>          [,1]     [,2]     [,3]
#> [1,] 6.714871 1.643449 6.530493
#> [2,] 1.643449 6.568033 2.312455
#> [3,] 6.530493 2.312455 9.028815
genPositiveDefMat("eigen",dim=3)
#> $egvalues
#> [1] 8.409136 4.076442 2.256715
#> 
#> $Sigma
#>            [,1]       [,2]      [,3]
#> [1,]  2.3217300 -0.1467812 0.5220522
#> [2,] -0.1467812  4.1126757 0.5049819
#> [3,]  0.5220522  0.5049819 8.3078880

Criado em 2019-10-27 pelo pacote reprex (v0.3.0)

Por fim, observe que uma abordagem alternativa é fazer uma primeira tentativa do zero e depois usar Matrix::nearPD()para tornar sua matriz positiva.

Matifou
fonte