LogLikelihood-Parameterschätzung für linearen Gaußschen Kalman-Filter

13

Ich habe einen Code geschrieben, der die Kalman-Filterung (unter Verwendung einer Reihe verschiedener Kalman-Filter [Information Filter et al.]) Für die lineare Gaußsche Zustandsraumanalyse für einen n-dimensionalen Zustandsvektor durchführen kann. Die Filter funktionieren sehr gut und ich bekomme eine schöne Ausgabe. Die Parameterschätzung über die Loglikelihood-Schätzung verwirrt mich jedoch. Ich bin kein Statistiker, sondern ein Physiker. Seien Sie also bitte freundlich.

Betrachten wir das lineare Gaußsche Zustandsraummodell

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

wobei yt unser Beobachtungsvektor ist, αt unser Zustandsvektor zum Zeitpunkt t . Die fett gedruckten Größen sind die Transformationsmatrizen des Zustandsraummodells, die gemäß den Merkmalen des betrachteten Systems eingestellt werden. Wir haben auch

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

wo . Jetzt habe ich die Rekursion für den Kalman-Filter für dieses generische Zustandsraummodell abgeleitet und implementiert, indem ich die Anfangsparameter und Varianzmatrizen H 1 und Q 1 erraten habe, wie ich Diagramme erstellen kannt=1,,nH1Q1

Kalman Filter

Wo die Punkte die Nilwasserstände für Jan. über 100 Jahre sind, ist die Linie der geschätzte Kalamn-Zustand, und die gestrichelten Linien sind die 90% -Konfidenzniveaus.

Für diesen 1D-Datensatz sind die Matrizen und Q t nur Skalare σ ϵ bzw. σ η . Jetzt möchte ich die richtigen Parameter für diese Skalare mithilfe der Ausgabe des Kalman-Filters und der Loglikelihood-Funktion ermittelnHtQtσϵση

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

Dabei ist der Zustandsfehler und F t die Zustandsfehlervarianz. Hier bin ich verwirrt. Aus dem Kalman-Filter habe ich alle Informationen, die ich brauche, um L zu berechnen, aber dies scheint mir nicht näher zu kommen, um die maximale Wahrscheinlichkeit von σ ϵ und σ η berechnen zu können . Meine Frage ist, wie ich die maximale Wahrscheinlichkeit von σ ϵ und σ η unter Verwendung des Loglikelihood-Ansatzes und der obigen Gleichung berechnen kann. Eine algorithmische Panne wäre für mich im Moment wie ein kaltes Bier ...vtFtLσϵσησϵση

Vielen Dank für Ihre Zeit.


Hinweis. Für den 1D-Fall ist und H t = σ 2 η . Dies ist das univariate Modell auf lokaler Ebene.Ht=σϵ2Ht=ση2

Mond Ritter
quelle

Antworten:

11

Wenn Sie den Kalman-Filter so ausführen, wie Sie es getan haben, erhalten Sie mit gegebenen Werten von und σ 2 η eine Folge von Innovationen ν t und deren Kovarianzen F t , daher können Sie den Wert von log L ( Y n ) mit berechnen die Formel, die Sie geben.σϵ2ση2νtFtlogL(Yn)

Mit anderen Worten, Sie können den Kalman-Filter als einen Weg betrachten, eine implizite Funktion von und σ 2 η zu berechnen . Das einzige, was Sie dann tun müssen, ist, diese Berechnung in eine Funktion oder ein Unterprogramm zu packen und diese Funktion für eine Optimierungsroutine zu verarbeiten - wie in R. Diese Funktion sollte als Eingaben σ 2 ϵ und σ 2 η akzeptieren und log L zurückgeben ( Y n ) .σϵ2ση2optimσϵ2ση2logL(Yn)

Einige Pakete in R (zB dlm) erledigen dies für Sie (siehe zB Funktion dlmMLE).

F. Tusell
quelle
Danke für deine Antwort. Ich schätze, dass ich anscheinend alle Komponenten habe, die für die explizite Berechnung der Loglikelihood erforderlich sind. Alle Referenzen scheinen jedoch darauf hinzudeuten, dass ich und σ η als Unbekannte in der Loglikelihood-Funktion verwende und diese mithilfe einer Newton-Methode maximiere. Das ist es, was mich verwirrt. "Loglikelihood wird numerisch in Bezug auf den unbekannten Zustandsvektor maximiert" - wie? σϵση
MoonKnight
Die Berechnung der Wahrscheinlichkeit ist nicht so explizit, da und σ η nicht explizit im Ausdruck von log L ( Y n ) vorkommen . Sie beeinflussen vielmehr die Wahrscheinlichkeit über ν t und F t . Daher müssen Sie den Kalman-Filter ausführen, um log L ( Y n ) für jedes Wertepaar von σ ϵ und σ η zu berechnen . Sobald Sie dies in Form einer Funktion codiert haben, können Sie es in eine Newton-Maximierungsfunktion (oder eine andere universelle Maximierungsfunktion) umwandeln.σϵσηlogL(Yn)νtFtlogL(Yn)σϵση
F. Tusell
1
Ich habe zufällig einen detaillierten Code (in R), der zeigt, wie dies genau für die Nil-Daten gemacht wird. Ich benutze es als Illustration für meine Schüler. Leider ist es in Spanisch, aber ich hoffe, der Code ist ziemlich klar (und ich kann die Kommentare übersetzen, wenn nicht). Sie können dieses Beispiel unter et.bs.ehu.es/~etptupaf/N4.html abrufen .
F. Tusell
Das ist enorm hilfreich. Vielen Dank für Ihre Zeit. Ihr Kommentar hat sehr geholfen! Manchmal ist es schwierig, "den Wald vor lauter Bäumen zu sehen", und nur eine einfache Erklärung ist erforderlich ... Nochmals vielen Dank.
MoonKnight
Ich möchte auch fragen, ob ich einen Blick auf die Seite werfen könnte, auf der Sie die Statusglättungsrekursion durchgehen. Deine Glättung sieht besser aus als meine und ich bin mir nicht sicher warum !? Ich habe versucht, es von Ihrer Website zu finden, aber ich kann die erforderliche Seite nicht finden ...
MoonKnight