Como garantir propriedades da matriz de covariância ao ajustar o modelo normal multivariado usando a máxima verossimilhança?

22

Suponha que eu tenha o seguinte modelo

yi=f(xi,θ)+εi

onde , é um vetor de variáveis ​​explicativas, são os parâmetros da função não linear e , onde é naturalmente matriz.yiRKxiθfεiN(0,Σ)ΣK×K

O objetivo é o usual para estimar e . A escolha óbvia é o método de máxima verossimilhança. A probabilidade de log para este modelo (assumindo que temos uma amostra ) pareceθΣ(yi,xi),i=1,...,n

l(θ,Σ)=n2log(2π)n2logdetΣi=1n(yif(xi,θ))Σ1(yf(xi,θ)))

Agora isso parece simples, a probabilidade de log é especificada, inserida em dados e usa algum algoritmo para otimização não linear. O problema é como garantir que Σ seja definitivo positivo. Usar, por exemplo, optimem R (ou qualquer outro algoritmo de otimização não linear) não garante que Σ seja definitivo positivo.

Portanto, a questão é como garantir que Σ permaneça positivo definitivamente? Eu vejo duas soluções possíveis:

  1. Reparametrize Σ como RR que R é uma matriz triangular superior ou simétrica. Então Σ sempre será positivo-definido e R pode ser irrestrito.

  2. Use a probabilidade do perfil. Derive as fórmulas para e . Comece com e itere \ hat {\ Sigma} _j = \ hat \ Sigma (\ hat \ theta_ {j-1}) , \ hat {\ theta} _j = \ hat \ theta (\ hat \ Sigma_ {j -1}) até convergência.θ^(Σ)Σ^(θ)θ0Σ^j=Σ^(θ^j1)θ^j=θ^(Σ^j1)

Existe alguma outra maneira e quanto a essas duas abordagens, elas funcionarão, são padrão? Isso parece um problema bastante comum, mas a pesquisa rápida não me deu nenhuma dica. Eu sei que a estimativa bayesiana também seria possível, mas no momento eu não gostaria de me engajar nela.

mpiktas
fonte
Eu tenho o mesmo problema em um algoritmo de Kalman, mas o problema é muito mais complicado e não é tão fácil de usar o truque de Hamilton. Gostaria de saber então se uma coisa mais simples a fazer seria simplesmente usar . Dessa forma, forço o código a não dar erro e não altero a solução. Isso também tem o benefício de forçar esse termo a ter o mesmo sinal que a parte final da probabilidade. Alguma ideia? log(detΣ+1)
econ_pipo 12/09

Respostas:

6

Supondo que, ao construir a matriz de covariância, você cuide automaticamente do problema de simetria, sua probabilidade de log será quando não for definido positivamente devido ao termo no modelo certo? Para evitar um erro numérico se eu pré-calcularia e, se não for positivo, torne a probabilidade do log igual a -Inf, caso contrário continue. Você deve calcular o determinante de qualquer maneira, para que isso não esteja lhe custando nenhum cálculo extra. Σlogdet Σdet Σ<0det Σ

Macro
fonte
5

Como se vê, você pode usar a probabilidade máxima do perfil para garantir as propriedades necessárias. Você pode provar que para um dado θ , l ( θ , Σ ) é maximizada porθ^l(θ^,Σ)

Σ^=1ni=1nε^iε^i,

Onde

ε^i=yif(xi,θ^)

Então é possível mostrar que

i=1n(yif(xi,θ^))Σ^1(yf(xi,θ^)))=const,

portanto, precisamos apenas maximizar

lR(θ,Σ)=n2logdetΣ^.

Naturalmente, neste caso, satisfará todas as propriedades necessárias. As provas são idênticas para o caso em que f é linear, o que pode ser encontrado na Análise de séries temporais por JD Hamilton, página 295; portanto, eu as omiti.Σf

mpiktas
fonte
3

Uma alternativa para a parametrização da matriz de covariância é em termos de valores próprios e p ( p - 1 ) / 2 ângulos "dados" θ i j .λ1,...,λpp(p1)/2θij

Ou seja, podemos escrever

Σ=GTΛG

onde é ortonormal eG

Λ=diag(λ1,...,λp)

com .λ1...λp0

Enquanto isso, pode ser parametrizado exclusivamente em termos de p ( p - 1 ) / 2 ângulos, θ i j , em que i = 1 , 2 , . . . , P - 1 e j = i , . . . , p - 1. [1]Gp(p1)/2θiji=1,2,...,p1j=i,...,p1

(detalhes a serem adicionados)

[1]: Hoffman, Raffenetti, Ruedenberg. "Generalização de ângulos de Euler para matrizes ortogonais em dimensões N". J. Math. Phys. 13, 528 (1972)

charles.y.zheng
fonte
GΣyi vector and the model function f(xi,θ) so that the errors are independent, then applying OLS to each of the rotated components (I think).
probabilityislogic
2

Along the lines of charles.y.zheng's solution, you may wish to model Σ=Λ+CC, where Λ is a diagonal matrix, and C is a Cholesky factorization of a rank update to Λ. You only then need to keep the diagonal of Λ positive to keep Σ positive definite. That is, you should estimate the diagonal of Λ and the elements of C instead of estimating Σ.

shabbychef
fonte
Can below diagonal elements in this settings be anything I want as long as the diagonal is positive? When simulate matrices this way in numpy not all of them are positive definite.
sztal
Λ is a diagonal matrix.
shabbychef