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?
15
Respostas:
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.C
R
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
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.V D2 D2 D C
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ânciaX I Y=VDX
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
OY 0 realizações em que a dimensão de Y é de 100 e a patente de C é 50 . A saída é10,000 Y 100 C 50
R
có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 ,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×10−5 C C
O código é exatamente paralelo à análise anterior e, portanto, deve ser autoexplicativo (mesmo para nãoD2 D
R
usuá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.fonte
Solução Método A :
No MATLAB, o código seria
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:
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.
fonte
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
fonte
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.
fonte
Você pode obter os resultados da função nearPD no pacote Matrix em R. Isso fornecerá uma matriz com valor real de volta.
fonte