PRML Kapitel 13 Wahrscheinlichste Schätzung Python-Implementierung des Hidden-Markov-Modells

In Kapitel 13 wird das Hidden-Markov-Modell als Modell für den Umgang mit Seriendaten vorgestellt. Dieses Mal werden wir die wahrscheinlichste Schätzung mit dem Hidden Markov Model (HMM) implementieren. Ich habe mehrmals versucht, die Parameter zu korrigieren und den verborgenen Zustand mithilfe des Vorwärts-Rückwärts-Algorithmus zu schätzen, aber ich habe die Parameter nie geschätzt, daher ist dies eine gute Gelegenheit, dies zu tun. Ich beschloss es zu sehen. HMM ist übrigens die Basis für Kalman-Filter und Partikelfilter.

Verstecktes Markov-Modell

Berücksichtigen Sie im Hidden-Markov-Modell die Seriendaten (z. B. die Vorder- und Rückseite des 100-maligen Werfens von Münzen). Und es wird gesagt, dass sich hinter diesen beobachteten Seriendaten eine latente Variable befindet (zum Beispiel wurde eine Münze verwendet, deren Rückseite (Vorderseite) leicht herauszukommen ist). Das grafische Modell des Hidden-Markov-Modells ist in der folgenden Abbildung dargestellt (aus PRML Abbildung 13.5). $ \ {x_i \} $ ist die beobachtete Variable (beide Seiten der Münze) und $ \ {z_i \} $ ist die latente Variable (Münzvorspannung). Die latente Variable hängt nur von der vorherigen latenten Variablen ab (z. B. nachdem eine Münze verwendet wurde, die leicht in der Tabelle angezeigt wird, werden Münzen, die auf dieselbe Tabelle vorgespannt sind, erneut verwendet), und die beobachtete Variable hängt zu diesem Zeitpunkt nur von der latenten Variablen ab. (Beispiel: Wenn nach vorne vorgespannte Münzen verwendet werden, ist die Vorderseite leicht herauszukommen.) Wenn eine latente Variable auf diese Weise nur von der vorherigen latenten Variablen abhängt, wird sie als Markov-Kette erster Ordnung bezeichnet.

Schauen wir uns die mathematischen Formeln genauer an. Die latente Variable $ z_n $ hängt nur von $ z_ {n-1} $ ab, und die Beziehung wird durch die Übergangswahrscheinlichkeitsmatrix $ A $ dargestellt. Jede Komponente der Matrix $ A $ repräsentiert die Übergangswahrscheinlichkeit der latenten Variablen. Die Komponente $ A_ {ij} $ ist $ j $ ($ z_), wenn sich $ z_ {n-1} $ im Zustand $ i $ ($ z_ {n-1, i} = 1 $) und $ z_n $ im Zustand $ j $ ($ z_) befindet {n, j} = 1 $). Angenommen, die latente Variable verwendet das Paar-K-Codierungsverfahren. Wenn Sie diese als Formeln schreiben,

p(z_n|z_{n-1},A) = \prod_{i=1}^K\prod_{j=1}^KA_{ij}^{z_{n-1,i},z_{nj}}

Es wird sein. Da jedoch die erste latente Variable $ z_1 $ keine latente Variable vor sich hat, wird der K-dimensionale Wahrscheinlichkeitsvektor $ \ pi $ verwendet.

p(z_1|\pi) = \prod_{i=1}^K\pi_i^{z_{1i}}

Wird besorgt.

Wenn die beobachtete Variable $ x $ eine diskrete Variable ist, die D diskrete Zustände annimmt (z. B. eine diskrete Variable, die zwei Zustände auf der Vorder- und Rückseite einer Münze hat), ist die bedingte Verteilung der beobachteten Variablen

p(x_n|z_n,\mu) = \prod_{i=1}^D\prod_{k=1}^K\mu_{ik}^{x_{ni}z_{nk}}

Es wird sein. Die Paar-K-Codierungsmethode wird jedoch auch für $ x $ verwendet.

HMM wahrscheinlichste Schätzung

In dem oben eingeführten Hidden-Markov-Modell sind drei Parameter $ \ pi, A, \ mu $ erschienen. Hier werden diese drei Parameter höchstwahrscheinlich aus den beobachteten Seriendaten $ X = \ {x_n \} _ {n = 1} ^ N $ geschätzt. Wenn diese drei Parameter zusammen als $ \ theta $ ausgedrückt werden, ist die zu maximierende Wahrscheinlichkeitsfunktion

p(X|\theta) = \sum_Zp(X,Z|\theta)

Es wird sein. $ Z = \ {z_n \} _ {n = 1} ^ N $. Die Zahl im Fall von Z ist die N-te Potenz der Anzahl der Zustände K von $ z $ (zum Beispiel: Wenn es zwei Vorspannungen von Münzen gibt und die Münze 100 Mal geworfen wird, $ 2 ^ {100} \ ca. 10 ^ {30} ($ Street) kann nicht direkt berechnet werden. Daher werden wir eine Methode verwenden, um den Rechenaufwand zu reduzieren, indem wir die Tatsache ausnutzen, dass es sich um ein Hidden-Markov-Modell erster Ordnung im Rahmen des EM-Algorithmus handelt.

Der EM-Algorithmus wiederholt die folgende Formelberechnung und -maximierung für $ \ theta $.

Q(\theta,\theta^{old}) = \sum_Zp(Z|X,\theta^{old})\ln p(X,Z|\theta)

M Schritt

Sei $ \ gamma (z_n) $ die periphere posteriore Verteilung der latenten Variablen $ z_n $ und $ \ xi (z_ {n-1}, z_n) $ die gleichzeitige posteriore Verteilung aufeinanderfolgender latenter Variablen.

\begin{align}
\gamma(z_n) &= p(z_n|X,\theta^{old})\\
\xi(z_{n-1},z_n) &= p(z_{n-1},z_n|\theta^{old})
\end{align}

Wenn Sie diese verwenden, um die zu maximierende Formel neu zu schreiben,

Q(\theta,\theta^{old}) = \sum_{k=1}^K\gamma(z_{1k})\ln\pi_k + \sum_{n=2}^N\sum_{j=1}^K\sum_{k=1}^K\xi(z_{n-1,j},z_{n,k})\ln A_{jk} + \sum_{n=1}^N\sum_{k=1}^K \gamma(z_{nk})\ln p(x_n|z_n,\mu)

Wenn dies für die Parameter $ \ pi, A, \ mu $ maximiert ist, verwenden die Parameter die Lagrange-Multiplikatormethode.

\begin{align}
\pi_k &= {\gamma(z_{1k})\over\sum_{k=1}^K\gamma(z_{1j})}\\
A_{jk} &= {\sum_{n=2}^N\xi(z_{n-1,j},z_{nk})\over\sum_{l=1}^K\sum_{n=2}^N\xi(z_{n-1,j},z_{nl})}\\
\mu_{ik} &= {\sum_{n=1}^N\gamma(z_{nk})x_{ni}\over\sum_{n=1}^N\gamma(z_{nk})}
\end{align}

Ich finde es.

E Schritt

In diesem Schritt berechnen wir die beiden hinteren Wahrscheinlichkeitsverteilungen $ \ gamma (z_n), \ xi (z_ {n-1}, z_n) $ für die in Schritt M verwendeten latenten Variablen. Die Methode wird als Vorwärts-Rückwärts-Algorithmus bezeichnet, und wie der Name schon sagt, kann die posteriore Wahrscheinlichkeitsverteilung effizient berechnet werden, indem zwei Berechnungen von vorne und hinten für die Seriendaten durchgeführt werden.

Nach vorne

Rekursiv $ \ alpha (z_n) = p (x_1, \ dots, x_n, z_n) $ mit $ \ alpha (z_1) = p (x_1, z_1) = p (x_1 | z_1) p (z_1) $ Berechnung.

\begin{align}
\alpha(z_n) &= p(x_1,\dots,x_n,z_n)\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_n,z_n,z_{n-1})\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_n|z_n,z_{n-1})p(z_n|z_{n-1})p(z_{n-1})\\
&= \sum_{z_{n-1}}p(x_1,\dots,x_{n-1}|z_{n-1})p(x_n|z_n)p(z_n|z_{n-1})p(z_{n-1})\\
&= p(x_n|z_n)\sum_{z_{n-1}}\alpha(z_{n-1})p(z_n|z_{n-1})
\end{align}

Rückwärts

\beta(z_N)=p(x_N|z_N)Wie\beta(z_n)=p(x_{n+1},\dots,x_N|z_n)Wird rekursiv berechnet.

\begin{align}
\beta(z_n)&=p(x_{n+1},\dots,x_N|z_n)\\
 &= \sum_{z_{n+1}}p(z_{n+1},x_{n+1},\dots,x_N|z_n)\\
&= \sum_{z_{n+1}}p(x_{n+1},\dots,x_N|z_n,z_{n+1})p(z_{n+1}|z_n)\\
&= \sum_{z_{n+1}}p(x_{n+2},\dots,x_N|z_{n+1})p(x_{n+1}|z_{n+1})p(z_{n+1}|z_n)\\
&= \sum_{z_{n+1}}\beta(z_{n+1})p(x_{n+1}|z_{n+1})p(z_{n+1}|z_n)
\end{align}

Ex-post-Wahrscheinlichkeitsverteilung

Von vorne und hinten oben

\begin{align}
\gamma(z_n) &= {\alpha(z_n)\beta(z_n)\over p(X)}\\
\xi(z_{n-1},z_n) &= {\alpha(z_{n-1})p(x_n|z_n)p(z_n|z_{n-1})\beta(z_n)\over p(X)}
\end{align}

Die posteriore Wahrscheinlichkeitsverteilung kann wie in erhalten werden.

Code

import Verwenden Sie wie gewohnt Numpy und Matplotlib.

import matplotlib.pyplot as plt
import numpy as np

Verstecktes Markov-Modell

class HiddenMarkovModel(object):

    def __init__(self, n_states_hidden, n_states_observe):
        #Anzahl der latenten variablen Zustände
        self.n_states_hidden = n_states_hidden

        #Anzahl der beobachteten variablen Zustände
        self.n_states_observe = n_states_observe

        #Initialisierung des Parameters pi der Verteilung der anfänglichen latenten Variablen
        self.initial = np.ones(n_states_hidden) / n_states_hidden

        #Initialisierung der Übergangswahrscheinlichkeitsmatrix A.
        self.transition = np.ones((n_states_hidden, n_states_hidden)) / (2 * n_states_hidden)
        self.transition += np.eye(n_states_hidden) * 0.5

        #Initialisierung des Parameters mu der Verteilung der beobachteten Variablen
        self.observation = np.random.rand(n_states_observe, n_states_hidden)
        self.observation /= np.sum(self.observation, axis=0, keepdims=True)

    # pi,A,Höchstwahrscheinlich Schätzung von mu
    def fit(self, sequence, iter_max=100):
        #EM-Schritt wiederholen
        for i in xrange(iter_max):
            params = np.hstack((self.transition.ravel(), self.observation.ravel()))

            #E Schritt
            p_hidden, p_transition = self.expectation(sequence)

            #M Schritt
            self.maximization(sequence, p_hidden, p_transition)

            #Überprüfen Sie, ob es konvergiert hat
            if np.allclose(params, np.hstack((self.transition.ravel(), self.observation.ravel()))):
                break

    #E Schritt
    def expectation(self, sequence):
        N = len(sequence)
        forward = np.zeros(shape=(N, self.n_states_hidden))
        # alpha(z_1)=p(x_1,z_1)Berechnung
        forward[0] = self.initial * self.observation[sequence[0]]
        backward = np.zeros_like(forward)
        # beta(z_N)=p(x_N|z_N)Berechnung
        backward[-1] = self.observation[sequence[-1]]

        #Nach vorne
        for i in xrange(1, len(sequence)):
            #PRML-Formel(13.36)
            forward[i] = self.transition.dot(forward[i - 1]) * self.observation[sequence[i]]

        #Rückwärts
        for j in xrange(N - 2, -1, -1):
            #PRML-Formel(13.38)
            backward[j] = (self.observation[sequence[j + 1]] * backward[j + 1]).dot(self.transition)

        #PRML-Formel(13.33)Latente Variable z_posteriores Wahrscheinlichkeitsverteilungsgamma von n(z_n)Berechnung
        p_hidden = forward * backward
        p_hidden /= np.sum(p_hidden, axis=-1, keepdims=True)

        #PRML-Formel(13.43)Gleichzeitige posteriore Wahrscheinlichkeitsverteilung aufeinanderfolgender latenter Variablen xi(z_{n-1},z_n)Berechnung
        p_transition = self.transition * (self.observation[sequence[1:]] * backward[1:])[:, :, None] * forward[:-1, None, :]
        p_transition /= np.sum(p_transition, axis=(1, 2), keepdims=True)

        return p_hidden, p_transition

    #M Schritt
    def maximization(self, sequence, p_hidden, p_transition):
        #PRML-Formel(13.18)Aktualisieren Sie die Parameter für die Verteilung der anfänglichen latenten Variablen
        self.initial = p_hidden[0] / np.sum(p_hidden[0])

        #PRML-Formel(13.19)Aktualisieren Sie die Übergangswahrscheinlichkeitsmatrix
        self.transition = np.sum(p_transition, axis=0) / np.sum(p_transition, axis=(0, 2))
        self.transition /= np.sum(self.transition, axis=0, keepdims=True)

        #PRML-Formel(13.23)Aktualisierung der Parameter des Beobachtungsmodells
        x = p_hidden[:, None, :] * (np.eye(self.n_states_observe)[sequence])[:, :, None]
        self.observation = np.sum(x, axis=0) / np.sum(p_hidden, axis=0)

Daten

#Erstellen Sie Münzwurfdaten
def create_toy_data(sample_size=100):

    #Werfen Sie Münzen nach Voreingenommenheit
    def throw_coin(bias):
        if bias == 1:
            #Tabelle ist 0.Wirf Münzen mit einer Wahrscheinlichkeit von 8
            return np.random.choice(range(2), p=[0.2, 0.8])
        else:
            #Die Rückseite ist 0.Wirf Münzen mit einer Wahrscheinlichkeit von 8
            return np.random.choice(range(2), p=[0.8, 0.2])

    #Bias der ersten Münze
    bias = np.random.uniform() > 0.5
    coin = []
    cheats = []
    for i in xrange(sample_size):
        coin.append(throw_coin(bias))
        cheats.append(bias)

        # 0.Ändern Sie die Vorspannung von Münzen mit einer Wahrscheinlichkeit von 01
        bias = bias + np.random.choice(range(2), p=[0.99, 0.01])
        bias = bias % 2
    coin = np.asarray(coin)

    return coin, cheats

Ganzer Code

import matplotlib.pyplot as plt
import numpy as np


class HiddenMarkovModel(object):

    def __init__(self, n_states_hidden, n_states_observe):
        self.n_states_hidden = n_states_hidden
        self.n_states_observe = n_states_observe
        self.initial = np.ones(n_states_hidden) / n_states_hidden
        self.transition = np.ones((n_states_hidden, n_states_hidden)) / (2 * n_states_hidden)
        self.transition += np.eye(n_states_hidden) * 0.5
        self.observation = np.random.rand(n_states_observe, n_states_hidden)
        self.observation /= np.sum(self.observation, axis=0, keepdims=True)

    def fit(self, sequence, iter_max=100):
        for i in xrange(iter_max):
            params = np.hstack((self.transition.ravel(), self.observation.ravel()))
            p_hidden, p_transition = self.expectation(sequence)
            self.maximization(sequence, p_hidden, p_transition)
            if np.allclose(params, np.hstack((self.transition.ravel(), self.observation.ravel()))):
                break

    def expectation(self, sequence):
        N = len(sequence)
        forward = np.zeros(shape=(N, self.n_states_hidden))
        forward[0] = self.initial * self.observation[sequence[0]]
        backward = np.zeros_like(forward)
        backward[-1] = self.observation[sequence[-1]]
        for i in xrange(1, len(sequence)):
            forward[i] = self.transition.dot(forward[i - 1]) * self.observation[sequence[i]]
        for j in xrange(N - 2, -1, -1):
            backward[j] = (self.observation[sequence[j + 1]] * backward[j + 1]).dot(self.transition)
        p_hidden = forward * backward
        p_hidden /= np.sum(p_hidden, axis=-1, keepdims=True)
        p_transition = self.transition * (self.observation[sequence[1:]] * backward[1:])[:, :, None] * forward[:-1, None, :]
        p_transition /= np.sum(p_transition, axis=(1, 2), keepdims=True)

        return p_hidden, p_transition

    def maximization(self, sequence, p_hidden, p_transition):
        self.initial = p_hidden[0] / np.sum(p_hidden[0])
        self.transition = np.sum(p_transition, axis=0) / np.sum(p_transition, axis=(0, 2))
        self.transition /= np.sum(self.transition, axis=0, keepdims=True)
        x = p_hidden[:, None, :] * (np.eye(self.n_states_observe)[sequence])[:, :, None]
        self.observation = np.sum(x, axis=0) / np.sum(p_hidden, axis=0)


def create_toy_data(sample_size=100):

    def throw_coin(bias):
        if bias == 1:
            return np.random.choice(range(2), p=[0.2, 0.8])
        else:
            return np.random.choice(range(2), p=[0.8, 0.2])

    bias = np.random.uniform() > 0.5
    coin = []
    cheats = []
    for i in xrange(sample_size):
        coin.append(throw_coin(bias))
        cheats.append(bias)
        bias = bias + np.random.choice(range(2), p=[0.99, 0.01])
        bias = bias % 2
    coin = np.asarray(coin)

    return coin, cheats


def main():
    coin, cheats = create_toy_data(200)

    hmm = HiddenMarkovModel(2, 2)
    hmm.fit(coin, 100)
    p_hidden, _ = hmm.expectation(coin)

    plt.plot(cheats)
    plt.plot(p_hidden[:, 1])
    for i in xrange(0, len(coin), 2):
        plt.annotate(str(coin[i]), (i - .75, coin[i] / 2. + 0.2))
    plt.ylim(-0.1, 1.1)
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

Wenn Sie den obigen Code ausführen, erhalten Sie möglicherweise das in der folgenden Abbildung gezeigte Ergebnis. Die horizontale Achse zeigt die Anzahl der Münzwürfe und die vertikale Achse zeigt die Wahrscheinlichkeit. Die blaue Linie zeigt die Vorspannung der tatsächlich verwendeten Münzen (zwei Arten von Münzen, eine, die leicht herauszukommen ist, und die andere, die leicht herauszukommen ist). Wird genutzt. Die grüne Linie stellt die Wahrscheinlichkeit dar, dass sich die zu diesem Zeitpunkt geworfene Münze in Zustand 1 befindet (ich weiß nicht, welche Münze sich in Zustand 1 befindet). 0,1 in der Grafik repräsentieren die Vorder- und Rückseite der tatsächlich beobachteten Münze. Es ist jedoch schwierig zu lesen, wenn alle angezeigt werden, sodass es nur alle zwei Male angezeigt wird.

Abhängig vom Anfangswert kann der EM-Algorithmus in der lokalen optimalen Lösung stecken bleiben, so dass das obige Ergebnis möglicherweise nicht immer erhalten wird.

Am Ende

Wenn die Parameter auch in Variablen geändert werden, passt sie häufig in eine schlechte lokale optimale Lösung. Daher kann es sinnvoll sein, anstelle der wahrscheinlichsten Schätzung Vorwissen einzuführen und die maximale posteriore Wahrscheinlichkeitsschätzung durchzuführen.