Gere números aleatórios distribuídos normalmente com matriz de covariância definida não positiva

15

Estimei a matriz de covariância da amostra e obtive uma matriz simétrica. Com C , eu gostaria de criar rn distribuído normal com n variáveis, mas, portanto, preciso da decomposição de C de Cholesky . O que devo fazer se C não for positivo definitivo?CCnCC

Klaus
fonte
11
Qual é a diferença com esta pergunta stackoverflow.com/questions/17295627/… ?
dickoa 9/07/2013
11
Matrizes semidefinidas positivas têm várias raízes quadradas (veja a explicação no final de stats.stackexchange.com/a/71303/919 , por exemplo). Você não precisa necessariamente daquele produzido pela decomposição de Cholesky. Aí reside o cerne do problema: encontre um método para calcular raízes quadradas que funcione mesmo quando a matriz é singular. @amoeba O título sugere que sua interpretação está correta.
whuber

Respostas:

8

Os questão diz respeito a como gerar variates aleatória de uma distribuição normal com uma multivariada (possivelmente) singular matriz de covariância . Esta resposta explica uma maneira que funcionará para qualquer matriz de covariância. Ele fornece uma implementação que testa sua precisão.CR


Análise algébrica da matriz de covariância

Como é uma matriz de covariância, é necessariamente simétrica e semidefinida positiva. Para completar a informação de fundo, vamos μ ser o vetor de meio desejado.Cμ

Como é simétrico, sua Decomposição de Valor Singular (SVD) e sua composição automática terão automaticamente a formaC

C=VD2V

para algumas matrizes ortogonais e matrizes diagonais D 2 . Em geral, os elementos diagonais de D 2 não são negativos (o que implica que todos têm raízes quadradas reais: escolha os positivos para formar a matriz diagonal D ). As informações que temos sobre C dizem que um ou mais desses elementos diagonais são zero - mas isso não afetará nenhuma das operações subseqüentes nem impedirá que o SVD seja calculado.VD2D2DC

Gerando valores aleatórios multivariados

Let tem uma distribuição normal multivariada padrão: cada componente tem média zero, variância unidade e todos os covariâncias são zero: a sua matriz de covariância é a identidade que eu . Então a variável aleatória Y = V D X tem matriz de covariânciaXIY=VDX

Cov(Y)=E(YY)=E(VDXXDV)=VDE(XX)DV=VDIDV=VD2V=C.

Por conseguinte, a variável aleatória tem uma distribuição normal com média multivariada μ e matriz de covariância C .μ+YμC

Código de exemplo e computação

O Rcódigo a seguir gera uma matriz de covariância de determinadas dimensões e classificações, analisa-a com o SVD (ou, em código comentado, com uma composição de eigend), usa essa análise para gerar um número especificado de realizações de (com vetor médio 0 ) e, em seguida, compara a matriz de covariância desses dados com a matriz de covariância pretendida, tanto numérica como graficamente. Como mostrado, gera 10 ,Y0 realizações em que a dimensão de Y é de 100 e a patente de C é 50 . A saída é10,000Y100C50

        rank           L2 
5.000000e+01 8.846689e-05 

Ou seja, a classificação dos dados também é e a matriz de covariância estimada a partir dos dados está à distância50 de C --que está próximo. Como uma verificação mais detalhada, os coeficientes de C são plotados contra os de sua estimativa. Todos estão próximos da linha da igualdade:8×105CC

Figura

O código é exatamente paralelo à análise anterior e, portanto, deve ser autoexplicativo (mesmo para não Rusuários, que podem emular no seu ambiente de aplicativo favorito). Uma coisa que revela é a necessidade de cautela ao usar algoritmos de ponto flutuante: as entradas de podem ser facilmente negativas (mas minúsculas) devido à imprecisão. Tais entradas precisam ser zeradas antes de calcular a raiz quadrada para encontrar D em si.D2D

n <- 100         # Dimension
rank <- 50
n.values <- 1e4  # Number of random vectors to generate
set.seed(17)
#
# Create an indefinite covariance matrix.
#
r <- min(rank, n)+1
X <- matrix(rnorm(r*n), r)
C <- cov(X)
#
# Analyze C preparatory to generating random values.
# `zapsmall` removes zeros that, due to floating point imprecision, might
# have been rendered as tiny negative values.
#
s <- svd(C)
V <- s$v
D <- sqrt(zapsmall(diag(s$d)))
# s <- eigen(C)
# V <- s$vectors
# D <- sqrt(zapsmall(diag(s$values)))
#
# Generate random values.
#
X <- (V %*% D) %*% matrix(rnorm(n*n.values), n)
#
# Verify their covariance has the desired rank and is close to `C`.
#
s <- svd(Sigma <- cov(t(X)))
(c(rank=sum(zapsmall(s$d) > 0), L2=sqrt(mean(Sigma - C)^2)))

plot(as.vector(C), as.vector(Sigma), col="#00000040",
     xlab="Intended Covariances",
     ylab="Estimated Covariances")
abline(c(0,1), col="Gray")
whuber
fonte
2
+1, mas quando você diz "indefinido" em sua primeira frase, o que exatamente você quer dizer? Eu verifiquei na Wikipedia e ele diz que semidefinido positivo não é indefinido, ou seja, indefinido significa que C tem valores próprios positivos e negativos. É isso o que você quer dizer aí?
Ameba diz Reinstate Monica
2
@amoeba Sim, isso foi um deslize. Obrigado por perceber. "Indefinido" significa que a assinatura da matriz possui sinais positivos e negativos, enquanto "semidefinido" significa que a assinatura possui apenas um sinal.
whuber
6

Solução Método A :

  1. Se C não for simétrico, simetrize-o. D <- 0.5(C+CT)
  2. D+(mmin(eigenvalue(D)))I

No MATLAB, o código seria

D = 0.5 * (C + C');
D =  D + (m - min(eig(CD)) * eye(size(D));

Método da Solução B : Formule e resolva um SDP Convexo (Programa Semidefinito) para encontrar a matriz D mais próxima de C de acordo com a norma frobenius de sua diferença, de modo que D seja definido positivamente, tendo especificado o valor próprio mínimo m.

Usando CVX no MATLAB, o código seria:

n = size(C,1);
cvx_begin
variable D(n,n)
minimize(norm(D-C,'fro'))
D -m *eye(n) == semidefinite(n)
cvx_end

Comparação de métodos de solução : além de simetrizar a matriz inicial, o método de solução A ajusta (aumenta) apenas os elementos diagonais em uma quantidade comum e mantém os elementos fora da diagonal inalterados. O método de solução B encontra a matriz definida positiva mais próxima (à matriz original) com o autovalor mínimo especificado, no sentido da norma frobenius mínima da diferença entre a matriz definida positiva D e a matriz original C, que é baseada nas somas de diferenças quadráticas de todos os elementos de D - C, para incluir os elementos fora da diagonal. Portanto, ajustando elementos fora da diagonal, isso pode reduzir a quantidade pela qual os elementos diagonais precisam ser aumentados e os elementos diagoanais não são necessariamente todos aumentados na mesma quantidade.

Mark L. Stone
fonte
2

Eu começaria pensando no modelo que você está estimando.

Se uma matriz de covariância não for semi-definida positiva, isso pode indicar que você tem um problema de colinearidade em suas variáveis, o que indicaria um problema com o modelo e não deve ser necessariamente resolvido por métodos numéricos.

Se a matriz não for positiva semidefinida por razões numéricas, existem algumas soluções que podem ser lidas aqui

johneric
fonte
11
A suposição é que o modelo é um modelo misto linear. E, neste caso, não é relevante encontrar um modelo correto para os dados, mas os dados são dados como exemplo para algum cálculo. Agora existe a possibilidade de você obter uma matriz semidefinida não positiva como estimativa para a covariância. Então, o que fazer a partir daí, se eu quiser descobrir a covariância da população distribuída normal de onde os dados vêm. Que a amostra é distribuída normalmente é a suposição.
Klaus
1

Uma maneira seria calcular a matriz a partir de uma decomposição de autovalor. Agora vou admitir que não conheço muito a matemática por trás desses processos, mas, a partir de minha pesquisa, parece proveitoso olhar para este arquivo de ajuda:

http://stat.ethz.ch/R-manual/R-patched/library/Matrix/html/chol.html

e alguns outros comandos relacionados em R.

Além disso, confira 'nearPD' no pacote Matrix.

Desculpe, não pude ter mais ajuda, mas espero que minha pesquisa possa ajudar a empurrá-lo na direção certa.

Frank P.
fonte
Oi, obrigado pelos links. No que diz respeito à decomposição do valor de eigen, essa decomposição não ajuda, porque a partir daí você obtém valores de autovalor complexos para matriz de raiz quadrada, mas eu preciso revelar matriz com valor.
Klaus
1

Você pode obter os resultados da função nearPD no pacote Matrix em R. Isso fornecerá uma matriz com valor real de volta.

library(Matrix)
A <- matrix(1, 3,3); A[1,3] <- A[3,1] <- 0
n.A <- nearPD(A, corr=T, do2eigen=FALSE)
n.A$mat

# 3 x 3 Matrix of class "dpoMatrix"
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.7606899 0.1572981
# [2,] 0.7606899 1.0000000 0.7606899
# [3,] 0.1572981 0.7606899 1.0000000
Dr. Mike
fonte
Para os usuários de R .. isso pode não ser uma versão ruim do "pobre homem" (com menos controle) do Método B da solução em minha resposta.
Mark L. Stone
Concordo que isso não é o ideal, mas às vezes funciona.
Dr. Mike