Comment garantir les propriétés de la matrice de covariance lors de l'ajustement d'un modèle normal multivarié en utilisant le maximum de vraisemblance?

22

Supposons que j'ai le modèle suivant

yi=f(xi,θ)+εi

yiRK , xi est un vecteur de variables explicatives, θ est les paramètres de la fonction non linéaire f et εiN(0,Σ) , où Σ est naturellement la matrice K×K

Le but est l'habituel d'estimer θ et Σ . Le choix évident est la méthode du maximum de vraisemblance. Log-vraisemblance pour ce modèle (en supposant que nous avons un échantillon ) ressemble(yi,xi),i=1,...,n

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

Maintenant, cela semble simple, la log-vraisemblance est spécifiée, insérée dans les données et utilise un algorithme pour l'optimisation non linéaire. Le problème est de savoir comment s'assurer que est défini positif. Utiliser par exemple dans R (ou tout autre algorithme d'optimisation non linéaire) ne me garantira pas que \ Sigma est défini positif.ΣΣoptimΣ

Donc, la question est de savoir comment s'assurer que Σ reste définitivement défini? Je vois deux solutions possibles:

  1. Reparametrise Σ as RRR est une matrice triangulaire supérieure ou symétrique. Alors Σ sera toujours positif-défini et R peut être sans contrainte.

  2. Utilisez la vraisemblance du profil. Dérivez les formules pour θ^(Σ) et Σ^(θ) . Commencez par θ0 et itérez Σ^j=Σ^(θ^j1) , θ^j=θ^(Σ^j1) jusqu'à convergence.

Y a-t-il une autre manière et qu'en est-il de ces 2 approches, vont-elles fonctionner, sont-elles standard? Cela semble un problème assez standard, mais la recherche rapide ne m'a donné aucun pointeur. Je sais que l'estimation bayésienne serait également possible, mais pour le moment je ne voudrais pas m'y engager.

mpiktas
la source
J'ai le même problème dans un algorithme de Kalman, mais le problème est beaucoup plus compliqué et pas aussi facile à utiliser l'astuce Hamilton. Je me demande alors si une chose plus simple à faire consisterait simplement à utiliser . De cette façon, je force le code à ne pas donner d'erreur et ne change pas la solution. Cela a également l'avantage de forcer ce terme à avoir le même signe que la dernière partie de la probabilité. Des idées? log(detΣ+1)
econ_pipo

Réponses:

6

En supposant que dans la construction de la matrice de covariance, vous êtes automatiquement en prenant soin de la question de symétrie, votre log-vraisemblance sera lorsque Σ est pas définie positive en raison du log d e t Σ terme dans le bon modèle? Pour éviter une erreur numérique si d e t Σ < 0, je précalculerais d e t Σ et, s'il n'est pas positif, je fais en sorte que la probabilité de log soit égale à -Inf, sinon continuez. Vous devez quand même calculer le déterminant, donc cela ne vous coûte aucun calcul supplémentaire. Σlogdet Σdet Σ<0det Σ

Macro
la source
5

Il s'avère que vous pouvez utiliser la probabilité maximale du profil pour garantir les propriétés nécessaires. Vous pouvez prouver que pour donné θ , l ( θ , Σ ) est maximisée parθ^l(θ^,Σ)

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

ε^i=yif(xi,θ^)

Il est alors possible de montrer que

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

il nous suffit donc de maximiser

lR(θ,Σ)=n2logdetΣ^.

Naturellement, dans ce cas, satisfera toutes les propriétés nécessaires. Les preuves sont identiques pour le cas où f est linéaire qui peut être trouvé dans l' analyse des séries temporelles par JD Hamilton page 295, donc je les ai omises.Σf

mpiktas
la source
3

Une alternative pour la paramétrisation de la matrice de covariance est en termes de valeurs propres et p ( p - 1 ) / 2 angles "Givens" θ i j .λ1,...,λpp(p1)/2θij

Autrement dit, nous pouvons écrire

Σ=GTΛG

est orthonormé, etG

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

avec .λ1...λp0

Pendant ce temps, peut être paramétré de manière unique en termes de p ( p - 1 ) / 2 angles, θ i j , où i = 1 , 2 , . . . , P - 1 et j = i , . . . , p - 1. [1]Gp(p1)/2θiji=1,2,...,p1j=i,...,p1

(détails à ajouter)

[1]: Hoffman, Raffenetti, Ruedenberg. "Généralisation des angles d'Euler aux matrices orthogonales à N dimensions". J. Math. Phys. 13, 528 (1972)

charles.y.zheng
la source
La matrice est en fait orthogonale, car Σ est une matrice symétrique. C'est l'approche que j'allais recommander - Fondamentalement, cela revient à faire tourner le vecteur y i et la fonction de modèle f ( x i , θ ) afin que les erreurs soient indépendantes, puis à appliquer l'OLS à chacune des composantes tournées (je pense). GΣyif(xi,θ)
probabilitéislogic
2

Dans la lignée de la solution de charles.y.zheng, vous pouvez modéliser , où Λ est une matrice diagonale et C est une factorisation de Cholesky d'une mise à jour de rang vers Λ . Il vous suffit alors de garder la diagonale de Λ positif pour garder Σ positif défini. Autrement dit, vous devez estimer la diagonale de Λ et les éléments de C au lieu d'estimer Σ .Σ=Λ+CCΛCΛΛΣΛCΣ

shabbychef
la source
Les éléments sous la diagonale dans ces paramètres peuvent-ils être tout ce que je veux tant que la diagonale est positive? Lorsque vous simulez des matrices de cette manière en numpy, toutes ne sont pas définies positives.
sztal
Λ