Estimation des paramètres LogLikelihood pour le filtre de Kalman gaussien linéaire

13

J'ai écrit du code qui peut effectuer le filtrage de Kalman (en utilisant un certain nombre de filtres de type Kalman différents [Information Filter et al.]) Pour l'analyse linéaire de l'espace d'état gaussien pour un vecteur d'état à n dimensions. Les filtres fonctionnent très bien et j'obtiens une belle sortie. Cependant, l'estimation des paramètres via l'estimation de loglik vraisemblance me prête à confusion. Je ne suis pas statisticien mais physicien, alors soyez gentil.

Considérons le modèle linéaire de l'espace d'état gaussien

yt=Ztαt+ϵt,
αt+1=Ttαt+Rtηt,

où est notre vecteur d'observation, notre vecteur d'état au pas de temps t α tytαtt . Les quantités en gras sont les matrices de transformation du modèle d'espace d'état qui sont fixées en fonction des caractéristiques du système considéré. Nous avons aussi

η tN I D ( 0 , Q t ) , α 1N I D ( a 1 , P 1 ) .

ϵtNID(0,Ht),
ηtNID(0,Qt),
α1NID(a1,P1).

. Maintenant, j'ai dérivé et implémenté la récursivité du filtre de Kalman pour ce modèle générique d'espace d'état en devinant les paramètres initiaux et les matrices de variance H 1 et Q 1t=1,,nH1Q1 je peux produire des tracés comme

Filtre de Kalman

où les points sont les niveaux d'eau du Nil pour janvier sur 100 ans, la ligne est l'état estimé de Kalamn et les lignes en pointillés sont les niveaux de confiance à 90%.

Or, pour cet ensemble de données 1D, les matrices et Q t ne sont que des scalaires σ ϵ et σ η respectivement. Alors maintenant, je veux obtenir les paramètres corrects pour ces scalaires en utilisant la sortie du filtre de Kalman et la fonction loglikelihoodHtQtσϵση

logL(Yn)=np2log(2π)12t=1n(log|Ft|+vtTFt1vt)

est l'erreur d'état et F t est la variance d'erreur d'état. Maintenant, voici où je suis confus. Du filtre de Kalman, j'ai toutes les informations dont j'ai besoin pour calculer L , mais cela ne me semble pas plus proche de pouvoir calculer la probabilité maximale de σ ϵ et σ η . Ma question est de savoir comment calculer la probabilité maximale de σ ϵ et σ ηvtFtLσϵσησϵση utilisant l'approche loglik vraisemblance et l'équation ci-dessus? Une panne algorithmique serait comme une bière fraîche pour moi en ce moment ...

Merci pour votre temps.


Remarque. Pour le cas 1D et H t = σ 2 η . Il s'agit du modèle univarié au niveau local.Ht=σϵ2Ht=ση2

MoonKnight
la source

Réponses:

11

Lorsque vous exécutez le filtre de Kalman comme vous l'avez fait, avec des valeurs données de etσ 2 η , vous obtenez une séquence d'innovations νtet leurs covariancesFt, vous pouvez donc calculer la valeur delogL(Yn)σϵ2ση2νtFtlogL(Yn) utilisant la formule que vous donnez.

En d'autres termes, vous pouvez considérer le filtre de Kalman comme un moyen de calculer une fonction implicite de et σ 2 ησϵ2ση2 . La seule chose que vous devez faire ensuite est de regrouper ce calcul dans une fonction ou un sous-programme et de gérer cette fonction dans une routine d'optimisation - comme optimdans R. Cette fonction doit accepter comme entrées et σ 2 η et renvoyer log L ( Y n ) .σϵ2ση2logL(Yn)

Certains packages dans R (par exemple dlm) le font pour vous (voir par exemple la fonction dlmMLE).

F. Tusell
la source
Merci pour votre réponse. J'apprécie que je semble avoir tous les composants nécessaires pour calculer explicitement la loglik vraisemblance, mais toutes les références que j'ai semblent suggérer que j'utilise et σ η comme inconnues dans la fonction loglik vraisemblance et maximiser cela en utilisant une méthode de type Newton? C'est ce qui me déroute; "la probabilité logicielle est maximisée numériquement par rapport au vecteur d'état inconnu" - comment? σϵση
MoonKnight
σϵσηlogL(Yn)νtFtlogL(Yn)σϵση
F. Tusell
1
Il se trouve que j'ai un code détaillé (en R) montrant comment le faire précisément pour les données du Nil. Je l'utilise comme illustration pour mes élèves. Malheureusement, il est en espagnol, mais j'espère que le code est assez clair (et je peux traduire les commentaires sinon). Vous pouvez récupérer cet exemple à partir de et.bs.ehu.es/~etptupaf/N4.html .
F. Tusell
C'est extrêmement utile. Merci beaucoup pour votre temps. Votre commentaire a beaucoup aidé! Parfois, il est difficile de "voir le bois pour les arbres" et avoir quelque chose de simple expliqué explicitement est tout ce qui est nécessaire ... Merci encore.
MoonKnight
Je voudrais également demander si je pourrais jeter un œil à la page où vous passez par la récursion de lissage d'état. Votre lissage est meilleur que le mien et je ne sais pas pourquoi!? J'ai essayé de le trouver sur votre site Web mais je ne trouve pas la page requise ...
MoonKnight