Stochastic Information Processing

by Daniel Pollithy

This is a summary of the lecture Stochastic Information Processing which I attended in the winter semester 2019/2020. The content has overlaps with LMA (localization of mobile agents), i.e. Kalman Filter. In SI it is taught how to convert a generative model into a probabilistic one and how to confront parametric integrals by means of approximations.

SI

Unfortunately, I wrote this in German. But the derivations could be interesting nevertheless.

Statische Systeme

Statische Systeme sind besonders einfache wertediskrete Systeme, in denen die Ausgabe y nur von der Eingabe u abhängt. Wobei u und y diskrete Zufallsvariablen sind. Ein statisches System S stellt eine bedingte Wahrscheinlichkeitsverteilung dar: $P(y=i|u=j)$. Von der Verteilung u zur Verteilung y. Da beide diskreten ZV auf endlichen Mengen definiert sind, können wir S als stochastische Matrix A darstellen. Diese enthält per Konvention $(P(y=i|u=j))_{i,j}$, also pro Zeile die Verteilung für ein konkretes u. A enthält nur nicht negative Einträge, jede Zeilensumme ist eins. Wenn zusätzlich auch noch die Spaltensummen ein ergeben, dann spricht man von einer doppelt-stochastischen Matrix.

Möchte man nun y aus u errechnen, so kann man so vorgehen:

\[P(y=i) \stackrel{(1)}{=} \sum_{j=1}^{J}{P(y=i, u=j)} \stackrel{(2)}{=} \sum_{j=1}^{J}{ P(y=i|u=j) \cdot P(u=j)} \stackrel{(3)}{=} A^T \cdot P(u)\]

(1) Satz der totalen W’keit. (2) Konditionierung. (3) Summe als Matrix.

Man sieht also, dass das System nicht deterministisch ist. Unsicherheiten werden weiter propagiert.

Dynamische Systeme (diskret)

Der Ausgang $\boldsymbol{y_k}$ hängt von Eingabe $\boldsymbol{u_k}$ und Zustand $\boldsymbol{x_k}$ ab. Der Zustand fasst die gesamte Historie zusammen. Man unterteilt dyn. Systeme in die Systemabbildung (dynamischer Teil) und die statische Messabbildung. Diskrete ZV $\boldsymbol{x_k}$ ist Markov-Kette erster Ordnung, falls gilt:

Systemabbildung

\[P(\boldsymbol{x_{k+1}}|\boldsymbol{x_{k}}, ... \boldsymbol{x_{0}}, \boldsymbol{u_{k}}, ...\boldsymbol{x_{0}}) = P(\boldsymbol{x_{k+1}}|\boldsymbol{x_{k}})\]

(1. Markoweigenschaft)

Man sagt: $\boldsymbol{x_{k+1}}$ ist bedingt unabhängig von $\boldsymbol{x_{k-1}}$ … $\boldsymbol{x_{0}}$. Wobei die Bedingung ist, dass $\boldsymbol{x_{k}}$ bekannt ist.

Achtung! $\boldsymbol{x_{k+1}}$ ist nicht unabhängig von $\boldsymbol{x_{0}}$. Nur bedingt unabhängig. Es folgt also, dass die zukünftige Entwicklung bedingt unabhängig von alten Zuständen ist, falls der aktuelle Zustand bekannt ist. Man betrachte zum Beispiel einen fallenden Ball. Kennt man den Zustand (Position, Geschwindigkeit, Beschleunigung, etc.) des Balles, so kann man seinen nächsten Zustand ohne Wissen über die Vergangenheit prädizieren.

Messabbildung

Der Zustand ist meist nicht verfügbar (Wortwahl ist Absicht). Man spricht von einer “latenten Variable”. Stattdessen kann man meist modellieren, wie $y_k$ sein müsste, wenn man einen gewissen Zustand $x_k$ und einen gewissen Eingang $u_k$ hätte. Diese Abbildung kann man ebenfalls als bedingte W’keit formulieren: $P(y_k=j|x_k=i, u_k=m)$. Der Einfluss von $u_k$ wird als “Durchgriff” bezeichnet, aber oft weggelassen. Wenn $y_k$ bedingt unabhängig von allen anderen ZV ist, falls $x_k$ gegeben ist, so nennt man dies die 2. Markoweigenschaft.

Die diskrete Messabbildung wird als Messmatrix B(i,j) bezeichnet (analog zu A.).

Hidden Markov Model (HMM)

Mit folgenden Zutaten kann nun also ein HMM beschrieben werden:

  1. Initialer Zustandsvektor
  2. Systemabbildung
  3. Messabbildung

Prädiktion, Filterung und Glättung

Die Aufgabe eines Schätzers in einem solchen System läge zum Beispiel darin, gegeben die Messungen bis zu einem Zeitpunkt t1, einen einzelnen oder mehrere Zustände zu einem beliebigen Zeitpunkt t2 zu schätzen.

Dabei unterscheidet man:

  • Prädiktion: t2 > t1
  • Filterung: t2 = t1 (Man prädiziert von (t2-1) zum aktuellen Zeitschritt und korriegiert dann seine Schätzung mit der ebenfalls bekannten Messung)
  • Glättung: t2 < t1.
    • “Fixed-Lag smoothing” ähnlich zu Filterung bloß kleiner Lag.
    • “Fixed-Point Smoothing”: Der Zeitpunkt, der geschätzt werden soll, steht fest. Es kommen aber immer noch mehr Daten hinein.
    • “Fixed-interval smoothing”: Man schätzt alle Zustände bis zum t2=t1-1, verwendet aber nur die Messungen, die bei den aktuellen Zeitpunkten zur Verfügung standen.
    • “Fixed-interval prediction”: Man hat nur die Messungen in einem festen Interval, es kommt nicht neues mehr dazu. Diese feste Anzahl an Punkten verwendet man dann um beliebige Prädiktionen zu machen.

Wonham Filter

Komplettes rekursives, diskretes Filter unter Verwendung von Einschrittprädiktion.

  • Dynamikmodell: $x_{t+1} = A^T x_t$
  • Messmodell: $y_t = B^T x_t$

Aus einer Messung auf den Zustand zu schließen, funktioniert wie folgt:

\[P(x_k|y_{0:k}) \stackrel{(1)}{=} P(x_k|y_k, y_{0:{k-1}}) \\ \stackrel{(2)}{=} \frac{P(y_k|x_k, y_{0:{k-1}}) P(x_k|y_{0:{k-1}})}{P(y_k|y_{0:{k-1}})} \\ \stackrel{(3)}{=} \frac{P(y_k|x_k, y_{0:{k-1}}) P(x_k|y_{0:{k-1}})}{\sum_{i}{P(y_k, x_k=i|y_{0:{k-1}})}} \\ \stackrel{(4)}{=} \frac{P(y_k|x_k, y_{0:{k-1}}) P(x_k|y_{0:{k-1}})}{\sum_{i}{P(y_k|x_k=i,y_{0:{k-1}}) P(x_k=i|y_{0:{k-1}})}} \\ \stackrel{(5)}{=} \frac{P(y_k|x_k) P(x_k|y_{0:{k-1}})}{\sum_{i}{P(y_k|x_k=i) P(x_k=i|y_{0:{k-1}})}} \\\]

(1) Nur andere Schreibweise

(2) Bayes Regel

(3) Satz der totalen W’keit

(4) Konditionierung

(5) 2. Markovannahme

Für ein konkret beobachtetes Symbol y’ bedeutet das, dass man aus der Matrix B die beobachtete Spalte entnimmt. Diese enthält die Likelihood von y’ unter allen möglichen Zuständen. Diesen multipliziert man dann mit vom Systemmodell prädizierten Zuständen und normiert. So erhält man einen durch die Beobachtungen verbesserte Schätzung über den Zustand.


# Beispiel: Wie viele Sterne (1-5) hat ein Buch wirklich verdient?
#           Böse Bots bewerten Bücher entweder mit 5 oder 1 Stern.
#           Menschen tendieren zu "nicht extremen" Bewertungen.
#              Also lieber 4 als 5 Sterne geben, obwohl das Buch sehr gut war.
#              Und wenn ein Buch schon gut bewertet wird, dann ist man umso kritischer.


# Starte gleichverteilt
x = np.array([0.2, 0.2, 0.2, 0.2, 0.2])[:, np.newaxis]

A = np.identity(5)

B = np.array([
    [0.60, 0.20, 0.1, 0.01, 0.09],  # Messungen, wenn das Buch 5 Sterne hätte
    [0.10, 0.70, 0.1, 0.01, 0.09],  # Messungen, wenn das Buch 4 Sterne hätte
    [0.09, 0.01, 0.8, 0.01, 0.09],  # Messungen, wenn das Buch 3 Sterne hätte
    [0.09, 0.01, 0.1, 0.70, 0.10],  # Messungen, wenn das Buch 2 Sterne hätte
    [0.09, 0.01, 0.1, 0.20, 0.60],  # Messungen, wenn das Buch 1 Sterne hätte
])

# Wie viele Sterne sollen wir auf unserer Website darstellen, wenn wir folgende
#    Bewertungen erhalten haben
star_ratings= [5, 1, 3, 5, 3, 4, 4, 1, 3, 4, 2]

# Rekursive Filterung

print("Kunde bewertet: \t\t 1 Stern \t 2 Sterne \t 3 Sterne \t 4 Sterne \t 5 Sterne")
print("-"*110)

for star_rating in star_ratings:    
    # step prediction
    x = A.T @ x

    # filter.
    # select the column of B for our observed measurement
    B_col = B[:, star_rating-1]

    x = (np.diag(B_col) @ x) / (np.ones(5) @  (np.diag(B_col) @ x))  

    print("{} Sterne \t\t->\t {:.2f} \t\t {:.2f} \t\t {:.2f} \t\t {:.2f} \t\t {:.2f}".format(
        star_rating, x[0,0], x[1,0], x[2,0], x[3,0], x[4,0])
         )

Kunde bewertet: 		 1 Stern 	 2 Sterne 	 3 Sterne 	 4 Sterne 	 5 Sterne
--------------------------------------------------------------------------------------------------------------
5 Sterne 		->	 0.09 		 0.09 		 0.09 		 0.10 		 0.62
1 Sterne 		->	 0.40 		 0.07 		 0.06 		 0.07 		 0.40
3 Sterne 		->	 0.28 		 0.05 		 0.34 		 0.05 		 0.28
5 Sterne 		->	 0.11 		 0.02 		 0.13 		 0.02 		 0.72
3 Sterne 		->	 0.06 		 0.01 		 0.54 		 0.01 		 0.38
4 Sterne 		->	 0.01 		 0.00 		 0.06 		 0.08 		 0.85
4 Sterne 		->	 0.00 		 0.00 		 0.00 		 0.25 		 0.74
1 Sterne 		->	 0.00 		 0.00 		 0.00 		 0.25 		 0.74
3 Sterne 		->	 0.00 		 0.00 		 0.02 		 0.25 		 0.73
4 Sterne 		->	 0.00 		 0.00 		 0.00 		 0.54 		 0.46
2 Sterne 		->	 0.00 		 0.00 		 0.00 		 0.54 		 0.46

Wertekontinuierliche Systeme

Wertekontinuierliche lineare Systeme

Linearität eines System ($A:x\rightarrow y$) besteht dann, wenn folgende zwei Eigenschaften vorliegen:

  1. Skalierung: $A(c \cdot x) = c \cdot A(x)$
  2. Superposition: $A(x_1 + x_2) = A(x_1) + A(x_2)$

In der Realität existieren lineare Systeme nur in gewissen Grenzen.

Statisches System

$ y_k = f(u_k) $, wobei $u_k \in R^{P}$ und $y_k \in R^{M}$

Für $f$ linear kann man schreiben: $y_k = A_k u_k$ mit $A \in R^{MxP}$.

Dieses generative Modell wird stochastisch, wenn $u_k$ eine W’keitsverteilung ist oder man endogenes Rauschen hinzufügt (also Rauschen im System). Dafür erweitert man $u_k$ zu $u_{k}^{*} = [u_k, w_k]$, wobei $w_k$ das Rauschen zusammenfasst.

Die Unsicherheiten in $u_k$ und $y_k$ beschreibt man oft mit den ersten beiden Momenten dieser vektoriellen ZV (Erwartungswert und Kovarianzmatrix).

Wir definieren $\hat{u}_k = E[u_k]$ und $C_k^u = Cov[u_k]$

Wichtige Erinnerung

Die Linearität des Erwartungswertes liefert für ein lineares System: $E[y_k] = E[A_k u_k] = A_k E[u_k] = A_k \hat{u}_k$.

Und die Kovarianz: $Cov[y_k] = Cov[y_k] = E[(E[y_k] - y_k)^2] = E[ (\hat{y}_k - y_k)^2 ] = E[ (A \hat{u}_k - A u_k)^2 ] = E[(A (\hat{u}_k - u_k))^2] = E[A ( \hat{u}_k - u_k)^2 A^T] = A E[ ( \hat{u}_k - u_k)^2 ] A^T = A Cov_k^u A^T$

Dynamische Systeme

Lineare Entwicklung des Zustandes: $x_{k+1} = A_k x_k + B_k u_k$ (Markovkette 1. Ordnung, zeitvariant, ohne k wäre es zeitinvariant).

Additives Rauschen: $u_k = \hat{u}_k + w_k$, wobei $w_k$ ein Zufallsvektor mit $E[w_k]=0$ und $Cov[w_k]=C_k^w$

Messmodell: Zustand nicht verfügbar. Achtung: Verfügbarkeit heißt, dass man eine Größe nicht messen kann, wobei Beobachtbarkeit bedeutet, dass man mit Messungen auf die Größe schließen kann. Der verdeckte Zustand ist nicht verfügbar, aber beobachtbar. Beispiel: Position eines Autos in Foto ist verfügbar, Geschwindigkeit nur über mehrere Frames beobachtbar.

Lineare Messabbilung: $y_k = H_k x_k + v_k$ wobei $v_k$ wieder additives Rauschen ist.

Kalman Filter

Prädiktionsschritt

$x_{k+1} = A_k x_k + B_k (\hat{u}_k + w_k)$

$E[x_{k+1}] = E[A_k x_k + B_k (\hat{u}_k + w_k)] = E[A_k x_k] + E[B_k (\hat{u}_k + w_k)] = A_k \hat{x}_k + B_k E[\hat{u}_k + w_k] = A_k \hat{x}_k + B_k E[\hat{u}_k + w_k] = A_k \hat{x}_k + B_k \hat{u}$

Notiz: $E[w_k] = 0$ per Definition. Wenn er nicht null wäre, aber bekannt, so könnte man ihn von $w_k$ abziehen. Wenn er nicht bekannt wäre, dann müsste man diesen ebenfalls schätzen.

$Cov[x_{k+1}] = Cov[A_k x_k + B_k (\hat{u}_k + w_k)] = A_k C_k^x A_k^T + B_k C_k^w B_k^T$

Trick:

Um die Kovarianz zu berechnen, definiert man $A_k^{} = [A_k, B_k]$ und $x_k^{} = [x_k, u_k]^T$.

Rauschen sei unabhängig von Zustand. Dann sieht $Cov[x_k^{*}]$ wie folgt aus: $\begin{bmatrix}C_k^x & 0\0 & C_k^w\end{bmatrix}$

Dann geht:

\[Cov[x_{k+1}^{*}] = E[(x_k^{*} - \hat{x_k^{*}})^2] = E[(x_k^{*} - \hat{x_k^{*}}) (x_k^{*} - \hat{x_k^{*}})^T]\] \[= E[A_k^{*} \cdot \begin{bmatrix} x_k - \hat{x_k} \\ u_k - \hat{u_k} \end{bmatrix} \cdot [(x_k - \hat{x_k})^T , (u_k - \hat{u_k})^T] \cdot {A_k^{*}}^T ]\] \[= A_k^{*} \cdot E[\begin{bmatrix} x_k - \hat{x_k} \\ u_k - \hat{u_k} \end{bmatrix} \cdot [(x_k - \hat{x_k})^T , (u_k - \hat{u_k})^T] ] \cdot {A_k^{*}}^T\]

Da die Kovarianz mittelwertbefreit ist, folgt:

\[= A_k^{*} \cdot Cov[ \begin{bmatrix} x_k \\ u_k \end{bmatrix} ] \cdot {A_k^{*}}^T\] \[= A_k^{*} \cdot \begin{bmatrix}C_k^x & 0\\0 & C_k^w\end{bmatrix} \cdot {A_k^{*}}^T\] \[= [A_k, B_k] \cdot \begin{bmatrix}C_k^x & 0\\0 & C_k^w\end{bmatrix} \cdot \begin{bmatrix}A_k^T \\ B_k^T \end{bmatrix}\] \[= A_k C_k^x A_k^T + B_k C_k^w B_k^T\]

Kalman Filterschritt

Lineare Messabbildung: $y_k = H_k x_k + v_k$ “Was für eine Messung erwarte ich?”

Die Kombination von Vorwissen und Messung ermöglicht es, auf einen höherdimensionalen verdeckten Zustand zu schließen:

$x_k^p = K_k^1 x_k + K_k^2 y_k$ (Der Filterschritt ist eine gewichtete Linearkombination)

Gesucht ist das BLUE filter (Best linear unbiased estimator). Linear ist er ja schon.

Unbiased: $E[x_k] = \tilde{x_k}$ Wobei $\tilde{x_k}$ die Grundwahrheit ist.

\[E[x_k^p] = E[K_k^1 x_k + K_k^2 y_k] = K_k^1 E[x_k] + K_k^2 E[y_k] = K_k^1 E[x_k] + K_k^2 H_k E[x_k] + E[v_k] = K_k^1 E[x_k] + K_k^2 H_k E[y_k] = K_k^1 \tilde{x_k} + K_k^2 H_k \tilde{x_k}\] \[K_k^1 \tilde{x_k} + K_k^2 H_k \tilde{x_k} = \tilde{x_k}\] \[K_k^1 + K_k^2 H_k = I\] \[K_k^1 = I - K_k^2 H_k\]

=> Dadurch kann man die Linearkombination mit nur einer Matrix $K_k$ ausdrücken:

\[x_k^p = (I - K_k^2 H_k) x_k + K_k^2 y_k = (I - K_k H_k) x_k + K_k y_k\]

Die Frage ist nun, wie bestimmen wir das optimale K? Wir suchen das K, für das die prädizierte Kovarianz $C_k^e$ minimal wird.

\[C_k^e = Cov[x_k^p] = (I - K_k H_k) C_k^x (I - K_k H_k)^T + K_k C_k^y K_k^T\]

$C_k^e$ ist abhängig von $K_k$. Daher schreiben wir: $C_k^p(K_k)$

$C_k^p(K_k)$ ist positiv definit und symmetrisch, das heißt, es handelt sich um eine Kovarianzmatrix. Wie minimiert man das?

Man könnte $det(C_k^p(K_k))$ oder $trace(C_k^p(K_k))$ minimieren. Beide sind proportional zum Flächeninhalt der Kovarianzellipse. Die Determinante ist 1/4 des umschließenden Rechteckes. Die Spur hat die Größe des größten Quaders, der in die Ellipse hineinpasst.

Wenn man die Determinante minimiert, dann ist es theoretisch möglich, dass eine Richtung 0 wird. Dadurch wird der Flächeninhalt null und die Covarianz ist minimal. Das ist unschön, weil man die Unsicherheit in alle Richtungen minimieren will.

Daher verwenden wir die Projektion der Matrix auf alle möglichen Einheitsvektoren e. Diese Projektion P(K) liefert ein skalares Gütemaß, welches wir minimieren.

\(P(K_k) = e^T C_k^e(K_k) e\).

Ableitungsregeln für Matrizen:

Matrix C:

  • $ d/dC (a^T C b) = a^T b $ ($a^T b$ ist diadisches Produkt, also eine singuläre Matrix. Rang=1. det=0. nicht invertierbar.)

  • $ d/dK (a^T K C K^T b) = ab^T K C^T + ba^T K C $ (KCK^T ist ein lineare Propagation, evtl. Basiswechsel. a^T … b ist eine Projektion.)

  • wenn a=b=e und C symm., dann: $ d/dK (e^T K C K^T e) = 2 ee^T KC $

jetzt können wir $P(K_k)$ ableiten:

\[\frac{d}{dK_k} e^T C_k^e(K_k) e =\] \[\frac{d}{dK_k} e^T ((I - K_k H_k) C_k^x (I - K_k H_k)^T + K_k C_k^y K_k^T) e =\] \[\frac{d}{dK_k} e^T ( (I C_k^x I^T) - (I C_k^x H_k^T K_k^T) - (K_k H_k C_k^x I^T) + (K_k H_k C_k^x H_k^T K_k^T) + K_k C_k^y K_k^T) e =\]

Erster Summand fällt weg, weil er kein K enthält.

\[\frac{d}{dK_k} e^T ( - (I C_k^x H_k^T K_k^T) - (K_k H_k C_k^x I^T) + (K_k H_k C_k^x H_k^T K_k^T) + K_k C_k^y K_k^T) e =\]

Projektion auf alle Summanden anwenden:

\[\frac{d}{dK_k} - e^T (I C_k^x H_k^T K_k^T) e - e^T (K_k H_k C_k^x I^T) e + e^T (K_k H_k C_k^x H_k^T K_k^T) e + e^T ( K_k C_k^y K_k^T) e\]
  1. Summand (Mit $D^T = C_k^x H_k^T$ und) \(\frac{d}{dK_k} - e^T (I C_k^x H_k^T K_k^T) e =\)

    \[(\frac{d}{dK_k^T} - (e^T (D^T K_k^T) e))^T = - (D ee^T)^T = - (H_k C_k^x ee^T)^T\]
  • 2 Summand: \(- ee^T (H_k C_k^x)^T\)

  • 3 Summand: \(2 ee^T K_k H_k C_k^x H_k^T\)

  • 4 Summand: \(2 ee^T K_k C_k^y\)

Daraus folgt:

\[\frac{d}{dK_k} - e^T (I C_k^x H_k^T K_k^T) e - e^T (K_k H_k C_k^x I^T) e + e^T (K_k H_k C_k^x H_k^T K_k^T) e + e^T ( K_k C_k^y K_k^T) e\] \[= - (H_k C_k^x ee^T)^T - ee^T (H_k C_k^x)^T + 2 ee^T K_k H_k C_k^x H_k^T + 2 ee^T K_k C_k^y\]

Setzen wir Null:

\[- (H_k C_k^x ee^T)^T - ee^T (H_k C_k^x)^T + 2 ee^T K_k H_k C_k^x H_k^T + 2 ee^T K_k C_k^y \stackrel{!}{=} 0\] \[= ee^T [-{C_k^x}^T H_k^T - {C_k^x}^T H_k^T + 2 K_k H_k C_k^x H_k^T + 2 K_k C_k^y ]\] \[= ee^T [-2{C_k^x}^T H_k^T + 2 K_k H_k C_k^x H_k^T + 2 K_k C_k^y ]\]

Wir multiplizieren beidseitig mit $({ee^T})^{-1}$:

\[-2{C_k^x}^T H_k^T + 2 K_k H_k C_k^x H_k^T + 2 K_k C_k^y = 0\] \[+ 2 K_k H_k C_k^x H_k^T + 2 K_k C_k^y = +2{C_k^x}^T H_k^T\]

Teilen durch 2:

\[K_k H_k C_k^x H_k^T + K_k C_k^y = {C_k^x}^T H_k^T\]

Klammern $K_k$ aus:

\[K_k (H_k C_k^x H_k^T + C_k^y) = {C_k^x}^T H_k^T\]

Und multiplizieren rechtseitig mit $(H_k C_k^x H_k^T + C_k^y)^{-1}$, um $K_k$ alleine zu stellen:

\[K_k = (C_k^x)^T H_k^T (H_k C_k^x H_k^T + C_k^y)^{-1}\]

Das Transponierte der Kovarianzmatrix ist immernoch die Kovarianzmatrix, da sie symmetrisch ist:

\[K_k = C_k^x H_k^T (H_k C_k^x H_k^T + C_k^y)^{-1}\]

$K_k$ heißt Kalman Gain und ist optimal, wenn:

  1. das Filter linear ist und
  2. die Zufallsvariablen normalverteilt sind (zum Beispiel additives, normalverteiltes Rauschen)

Somit ergibt sich:

\[E[x_k^e] = (I - K_k H_k) E[x_k^p] + K_k E[y_k]\]

Und:

\(C_k^e = Cov[x_k^p] = (I - K_k H_k) C_k^x (I - K_k H_k)^T + K_k C_k^y K_k^T\) \(...\) \(C_k^e = C_k^p - K_k H_k C_k^p\)

Für intuitive Anschauungen eignet es sich die Feedback Form zu betrachten:

\[\hat{x_k^e} = \hat{x_k^p} + K_k (\hat{y_k} - H_k \hat{x_k^p})\]

Hieran lässt sich erkennen, dass der geschätzte Zustand nach einem Messupdate eine Kombination des priori Zustandes mit dem neuen Wissen aus der Messung ist. Wenn das Kalman Gain hoch ist, dann wird die Messung mit geringer Unsicherheit gesehen und somit werden starke Abweichungen der Messung von der erwarteten Messung stark auf den neuen Zustand einfließen. Wenn das Messrauschen unendlich groß wäre, dann wäre $K_k=0.$

Allgemeine Systeme

Statische Systeme

Generatives Modell (Werte abgebildet auf Werte): \(\underline{\mathbf{y}_k} = a(\underline{\mathbf{u}_k})\)

Jetzt sind aber $u_k$ und $y_k$ durch die W’keitsverteilungsdichtefunktionen (kurz “Dichten”) $f_k^u(\underline{\mathbf{u}_k})$ und $f_k^y(\underline{\mathbf{y}_k})$

Würde man nun einfach trotzdem mit Erwartungswert und Kovarianz weiterrechnen, anstatt die Dichten zu betrachten, dann wäre das eine Linearisierung, wie sie beim EKF vorgenommen wird.

Gesucht ist also die Dichte $f_k^y(\underline{\mathbf{y}_k})$ zur gegebenen Dichte $f_k^u(\underline{\mathbf{u}_k})$. Ein Abbildung von einer Dichte zu einer anderen Dichte nennen wir probabilistisches Modell.

Dynamische Systeme

Diese zerfallen wieder in das Beobachtungsmodell und die Systemabbildung.

Ein “einfacheres” nichtlineares Systemmodell: \(x_{k+1} = a_k(x_k, \hat{u}_k, w_k)\)

Es wäre schwieriger, wenn $x_{x+1}$ auch auf der rechten Seite vorkäme, dann spricht man von einem impliziten System. Dieser Fall wurde in LMA behandelt.)

Für additives Rauschen kommt man leicht vom generativen zu einem probabilistischen Modell:

\[x_{k+1} = a(x_k, \hat{u}_k) + w_k\]

Das Systemrauschen \(w_k\) ist hier additiv und wird beschrieben durch die Dichte $f_k^w(w_k)$.

Zusätzlich gehen wir davon aus, dass das Systemrauschen:

  1. Gaussverteil ist mit bekannten Parametern.
  2. $f_k^w(w_k)$ für unterschiedliche k unabhängig voneinander sind. Im Falle von Gaußverteilung reicht unkorreliert, da in diesem Fall unabhängig daraus folgt. Dann gilt für die gemeinsame Verteilung zweier Rauschdichten: $f_{1,2}^w(w_1, w_2) = f_1^w(w_1) \cdot f_2^w(w_2)$

Bei der Messabbildung kann genauso verfahren werden: \(y_k = h_k(x_k) + v_k\)

Was für Dichten können wir verwenden?

  • Kontinuierliche analytische Dichten
    • Gaußdichte
    • GMM
    • Exponentialdichte (Exponenten komplexer machen)
    • Edgeworth-Entwicklungen (Gaußdichte mit Polynom multipliziert. Nichts globales, aber lokal bekommt man ein paar Hügel rein).
  • Diskrete Dichten auf kontinuierlichen Domänen
    • Dirac-Mixture Densities

Blick auf Gaußdichte

2d: \(f_{x,y}(x,y) = N(\begin{bmatrix}x\\y\end{bmatrix} - \begin{bmatrix}\hat{x}\\\hat{y}\end{bmatrix}, \begin{bmatrix}\sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{bmatrix})\)

Wobei $\rho \in [-1, +1]$ als Korrelationskoeffizient bezeichnet wird. Wenn dieser =+1 ist, dann lieft eine deterministische Abhängigkeit vor.

Man schreibt manchmal: \(f_{x,y}(x,y) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{(1-\rho^2)}} \exp(-\frac{1}{2}Q(x,y))\)

mit \(Q(x,y) = \frac{1}{1-\rho} \Big( \frac{(x-\hat{x})^2}{\sigma_x^2} + \frac{(y-\hat{y})^2}{\sigma_y^2} - 2 \rho \frac{x-\hat{x}}{\sigma_x} \frac{y-\hat{y}}{\sigma_y} \Big)\)

Maßgrößen für Kovarianzen:

  • die Determinante von Cov det(C) ist proportional zum Volumen der Kovarianz-Körpers, da die Determinante das Produkt der Eigenwerte ist. Somit ist im 2d-Fall zum Beispiel 4*det(C) die Fläche des umschließenden Rechtecks.
  • Spur (Summe der Quadrate der Diagonale): Nach Pythagoras ist die Spur einer diagonalen Kovarianzmatrix das Quadrat der der Länge … -> Wurzel … Vorteil: Spur invariant ggü. unitärer Rotation. Selbst wenn ein Eigenwert mal 0 sein sollte, dann ist die Spur immernoch nicht null.

Moment Matching

Berechne die Verteilung einer Population anhand ihrer Momente. Angenommen wir haben eine Stichprobe aus der Population genommen. Nun könnten wir zum Beispiel die empirischen 1. und 2. Momente berechnen. Empirischer Mittelwert und empirische Varianz. Nun wählen wir eine beliebige Verteilungsfunktion. Jetzt können wir Parameter für die Verteilung bestimmen, die empirischen Momente liefern. Diese Methode ist besser per Hand zu rechnen und kann zum Beispiel auch als Initialisierung für eine Maximum-Likelihood-Schätzung verwendet werden.

Theorem: Für eine gegebene Dichte $\tilde{f}(x)$ führt die Wahl von $m=E_{\tilde{f}}[x]$ und $\sigma^2=Cov_{\tilde{f}}[x]$, zu einer Minimierung der Kullback-Leibler-Divergenz zwischen $\tilde{f}$ und der Normalverteilung.

\[KL = \int_{R} { \tilde{f}(x) log\big( \frac{\tilde{f}(x)}{f(x)} \big) } dx\]

Offene Frage: Warum $\tilde{f}(x)$ als Vorfaktor?

Das i-te zentrale Momente einer Gaußdichte lautet:

\[C_i = E[(x-x)^{i}]\]

Ein zentrales Moment ist das gleiche, wie ein “normales” Moment, bloß, dass man die Dichte um Null zentriert hat.

ToDo: Beschreibe, wie man Momente generell berechnen kann.

GMM: $f(x) = \sum_{i=1}^{L}{w_i N(\hat{x}_i, \sigma_i)}$

  • Anzahl Moden ist ungleich Anzahl Komponenten -> Maximum bestimmen nicht trivial
  • Wir akzeptieren nur $w_i > 0$. Theoretisch geht es aber auch anders. Und Summe der Gewichte muss 1 sein.
  • Ein GMM mit L Komponenten hat 3L-1 Freiheitsgrade

Diracsche Deltafunktion

$\delta(x-x_0) = $ …

  • nicht definiert, für $x=x_0$
  • 0, für alles andere

Mit $\int_{R}{\delta(x-x_0) dx} = 1$

Sie besitzt die Ausblendeigenschaft:

  • $\int f(x) \delta(x-x_0) dx = f(x_0)$
  • Auch: $f(x) \cdot \delta(x-x_0) = f(x_0) \cdot \delta(x-x_0)$
    • Fall $x=x_0$: $f(x_0) = f(x_0)$
    • Sonst: 0

Allgemein gilt:

\[\delta(f(x)) = \sum_{i=1}^{N}{\frac{1}{|f'(x_i)|} \delta(x-x_i) }\]

$x_i$ sind die N Nullstellen der Funktion f(x). $f(x) \ne 0$

Beweis dazu: https://math.stackexchange.com/questions/276583/dirac-delta-function-of-a-function.

(Eine sehr umfassende Sammlung von Eigenschaften der diracschen Delta-Funktion: https://mathworld.wolfram.com/DeltaFunction.html)

Beispiele generatives Modell zu bedingter Dichte:

$y = 3 * x$ —–> $f(y x) = \delta(y-3x)$
$y = a(x) + w$ —–> $f(y x,w) = \delta(y-a(x)-w)$

Bei den bedingten Dichten handelt es sich also um dirac distributionen, die nur dort ungleich Null sind, wo die Funktion wahr ist.

Additives Rauschen

\[x_{k+1} = a_k(x_k) + w_k\] \[f^{x_{k+1}}(x_{k+1}| x_k, w_k) = \delta(a_k(x_k) + w_k - x_{k+1})\] \[f^{x_{k+1}}(x_{k+1}| x_k) = \int{f^{x_{k+1}, w_k}(x_{k+1}, w_k| x_k) d w_k}\] \[f^{x_{k+1}}(x_{k+1}| x_k) = \int{f^{x_{k+1}}(x_{k+1}|x_k, w_k) \cdot f^w(w_k|x_k) d w_k}\]

Hier nehmen wir jetzt an, dass das Rauschen unabhängig von den Daten ist.

\[f^{x_{k+1}}(x_{k+1}| x_k) = \int{f^{x_{k+1}}(x_{k+1}|x_k, w_k) \cdot f^w(w_k) d w_k}\]

Jetzt Dirac einsetzen:

\[f^{x_{k+1}}(x_{k+1}| x_k) = \int{\delta(a_k(x_k) + w_k - x_{k+1}) \cdot f^w(w_k) d w_k}\]

Für welches $w_k$ ist der Dirac-Stoß ungleich Null? (vgl. $\delta(w_k - w_{k,0})$)

\[a_k(x_k) + w_k - x_{k+1} = 0\] \[w_k = x_{k+1} - a_k(x_k)\]

Durch Siebeigenschaft des Dirac folgt:

\[f^{x_{k+1}}(x_{k+1}| x_k) = f^w(x_{k+1} - a_k(x_k) )\]

Was eine bequeme Form für additives Rauschen ist. Wenn das Rauschen zum Beispiel normalverteilt ist, dann ist die Transitionsdichte für ein gegebens x_k ebenfalls normalverteilt.

Multiplikatives Rauschen

$ x_{k+1} = x_k \cdot w_k $ mit $w_k \sim f_k^w(w_k)$

\[f(x_{k+1}|x_k, w_k) = \delta(g(w_k))\] \[g(w_k) = x_k w_k - x_{k+1}\]

Finde Nullstellen von g:

\[x_k w_k - x_{k+1} = 0\] \[x_k w_k = x_{k+1}\] \[w_k = \frac{x_{k+1}}{x_k}\] \[\frac{d}{d w_k} g(w_k) = x_k\]

Daraus folgt:

\[\delta(g(w_k)) = \frac{1}{|x_k|} \delta(w_{k} - \frac{x_{k+1}}{x_k})\] \[f(x_{k+1}|x_k) = \int{f(x_{k+1}|x_k, w_k) f_k^w(w_k) d w_k}\] \[f(x_{k+1}|x_k) = \int{\delta(g(w_k)) f_k^w(w_k) d w_k}\] \[f(x_{k+1}|x_k) = \int{\frac{1}{|x_k|} \delta(w_{k} - \frac{x_{k+1}}{x_k}) f_k^w(w_k) d w_k}\] \[f(x_{k+1}|x_k) = \frac{1}{|x_k|} f_k^w(\frac{x_{k+1}}{x_k})\]

Visualisierung der Transitionsdichte

y = a(x) + w

y = x + N(0, 1)

\[f(y|x) = \int{f(y|x,w) \cdot f(w) dw} = f^w(y - a(x)) = f^w(y - x)\] \[f(y|x) = N(y - x)\]
from scipy.stats import norm
%matplotlib notebook

step = 0.5

x = np.arange(0, 10, step)
y = np.arange(0, 10, step*0.1)
z = np.arange(0, 10, step)

X, Y = np.meshgrid(x, y)

f = norm.pdf(Y-X, 0.1)



fig = plt.figure()
ax = fig.gca(projection='3d')

surf = ax.plot_surface(X, Y, f, cmap=cm.coolwarm,
                       linewidth=0, antialiased=True)

# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.xlabel("$x_k$")
plt.ylabel("$y_k$")
ax.set_zlabel("$P(y|x)$")

plt.show()
<IPython.core.display.Javascript object>

Dieser Plot ist so zu sehen, dass für jedes gegeben $\hat{x}_k$ eine Normalverteilung vorliegt. Die Verbundwahrscheinlichkeit von $P(x,y)$ im Diskreten oder $f^{x,y}(x,y)$ im Kontinuierlichen hängt von $x_k$ ab. Im Diskreten kann man folgenden Satz nutzen: $P(X \cap Y) = P(X|Y) \cdot P(Y) = P(Y|X) \cdot P(X)$.

Chapman-Kolmogorov-Integral (CKI)

Betrachten wir einen stochastischen Prozess $(x)_k$ bei dem wir die 1. Markov-Eigenschaft annehmen.

\[f(x_{k+1}) = \int_{R}{ f(x_{k+1}|x_k) \cdot f(x_k) d x }\]

Wir berechnen die Transitionsdichte zunächst weiterhin mit der dirachsen Delta-Funktion. Für den Trivialfall, dass $f(x_k)$ deterministisch (also dirac-verteilt) ist, löst sich das Intergral noch von selbst auf.

Das Hauptproblem ist, dass es sich beim CKI um ein Parameterintegral handelt. Das heißt, es ist nicht nur abhängig von der Integriervariable $x_k$, sondern auch noch von anderen “freien” Variablen. In diesem Fall $x_{k+1}$.

Beobachtungen einfließen lassen

Gegeben einer Beobachtung $\hat{y}$ müssen wir uns entscheiden, was für eine Art von Schätzer wir haben wollen:

Maximum-Likelihood: $$arg_x max f(\hat{y} x)$$
Maximum-A-Posteriori: $$arg_x max f(x \hat{y})$$

Angenommen wir haben das folgende generative Modell: $y = x^2 + v$ Dann erhalten wir \(f(y|x,v) = \delta(y-x^2-v)\)

Interessieren wir uns für $f(x \hat{y})$:
\[f(x|\hat{y}) = \frac{f(x, \hat{y})}{f(\hat{y})} = \frac{f(\hat{y}|x) f(x)}{f(\hat{y})}\]