PRML Chapter 13 Estimation la plus probable Implémentation Python du modèle de Markov caché

Au chapitre 13, le modèle de Markov caché est présenté comme un modèle de gestion des données de série. Cette fois, nous implémenterons l'estimation la plus probable avec le modèle de Markov caché (HMM). J'ai fait plusieurs fois pour corriger les paramètres et estimer l'état caché à l'aide de l'algorithme avant-arrière, mais je n'ai jamais estimé les paramètres, c'est donc une bonne occasion de le faire. J'ai décidé de le voir. À propos, HMM est la base des filtres de Kalman et des filtres à particules.

Modèle de Markov caché

Dans le modèle de Markov caché, considérez les données de série (par exemple, le recto et le verso du lancement de pièces 100 fois). Et, on dit qu'il y a une variable latente (par exemple, une pièce dont le verso (recto) est facile à sortir a été utilisée) derrière ces données de séries observées. Le modèle graphique du modèle de Markov caché est illustré dans la figure ci-dessous (à partir de la figure 13.5 de PRML). $ \ {x_i \} $ est la variable observée (les deux côtés de la pièce), et $ \ {z_i \} $ est la variable latente (biais de pièce). La variable latente dépend uniquement de la variable latente précédente (par exemple, après l'utilisation d'une pièce qui est facile à apparaître dans le tableau, les pièces biaisées vers la même table sont réutilisées), et la variable observée dépend uniquement de la variable latente à ce moment. (Exemple: si des pièces biaisées vers l'avant sont utilisées, l'avant est facile à sortir). Lorsqu'une variable latente ne dépend que de la variable latente précédente de cette manière, elle est appelée chaîne de Markov du premier ordre.

Regardons de plus près en utilisant des formules mathématiques. La variable latente $ z_n $ ne dépend que de $ z_ {n-1} $, et la relation est représentée par la matrice de probabilité de transition $ A $. Chaque composante de la matrice $ A $ représente la probabilité de transition de la variable latente. Le composant $ A_ {ij} $ est $ j $ ($ z_) lorsque $ z_ {n-1} $ est dans l'état $ i $ ($ z_ {n-1, i} = 1 $) et $ z_n $ est dans l'état $ j $ ($ z_ {n, j} = 1 $). Cependant, supposons que la variable latente utilise la méthode de codage par paire K. Si vous les écrivez sous forme de formules,

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

Ce sera. Cependant, comme la première variable latente $ z_1 $ n'a pas de variable latente devant elle, le vecteur de probabilité K-dimensionnel $ \ pi $ est utilisé.

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

ça ira.

Si la variable observée $ x $ est une variable discrète qui prend D états discrets (par exemple, une variable discrète qui a deux états à l'avant et à l'arrière d'une pièce), la distribution conditionnelle des variables observées est

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

Ce sera. Cependant, la méthode de codage paire K est également utilisée pour $ x $.

Estimation la plus probable de HMM

Dans le modèle de Markov caché présenté ci-dessus, trois paramètres $ \ pi, A, \ mu $ sont apparus. Ici, ces trois paramètres sont très probablement estimés à partir des données de la série observée $ X = \ {x_n \} _ {n = 1} ^ N $. Si ces trois paramètres sont exprimés collectivement par $ \ theta $, la fonction de vraisemblance à maximiser est

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

Ce sera. Cependant, $ Z = \ {z_n \} _ {n = 1} ^ N $. Le nombre dans le cas de Z est la Nième puissance du nombre d'états K de $ z $ (par exemple: s'il y a deux biais de pièces et que la pièce est lancée 100 fois, $ 2 ^ {100} \ approx10 ^ {30} ($ Street), il ne peut pas être calculé directement. Par conséquent, nous utiliserons une méthode pour réduire la quantité de calcul en profitant du fait d'être un modèle de Markov caché du premier ordre dans le cadre de l'algorithme EM.

L'algorithme EM itère le calcul de formule et la maximisation suivants pour $ \ theta $.

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

Étape M

Soit $ \ gamma (z_n) $ la distribution postérieure périphérique de la variable latente $ z_n $, et $ \ xi (z_ {n-1}, z_n) $ la distribution postérieure simultanée de variables latentes consécutives.

\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}

En utilisant ces derniers, si vous réécrivez la formule à maximiser,

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)

Si cela est maximisé pour les paramètres $ \ pi, A, \ mu $, alors les paramètres utiliseront la méthode du multiplicateur de Lagrange.

\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}

Je peux le trouver.

Étape E

Dans cette étape, nous allons calculer les deux distributions de probabilités postérieures $ \ gamma (z_n), \ xi (z_ {n-1}, z_n) $ pour les variables latentes utilisées à l'étape M. La méthode s'appelle l'algorithme avant-arrière et, comme son nom l'indique, la distribution de probabilité postérieure peut être calculée efficacement en effectuant deux calculs de l'avant et de l'arrière pour les données de la série.

Vers l'avant

Récursivement $ \ alpha (z_n) = p (x_1, \ dots, x_n, z_n) $ avec $ \ alpha (z_1) = p (x_1, z_1) = p (x_1 | z_1) p (z_1) $ Calculer.

\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}

En arrière

\beta(z_N)=p(x_N|z_N)Comme\beta(z_n)=p(x_{n+1},\dots,x_N|z_n)Est calculé récursivement.

\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}

Distribution de probabilité ex post

De l'avant et de l'arrière au-dessus

\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}

La distribution de probabilité postérieure peut être obtenue comme indiqué dans.

code

import Utilisez numpy et matplotlib comme d'habitude.

import matplotlib.pyplot as plt
import numpy as np

Modèle de Markov caché

class HiddenMarkovModel(object):

    def __init__(self, n_states_hidden, n_states_observe):
        #Nombre d'états de variables latentes
        self.n_states_hidden = n_states_hidden

        #Nombre d'états variables observés
        self.n_states_observe = n_states_observe

        #Initialisation du paramètre pi de la distribution de la variable latente initiale
        self.initial = np.ones(n_states_hidden) / n_states_hidden

        #Initialisation de la matrice de probabilité de transition A
        self.transition = np.ones((n_states_hidden, n_states_hidden)) / (2 * n_states_hidden)
        self.transition += np.eye(n_states_hidden) * 0.5

        #Initialisation du paramètre mu de distribution des variables observées
        self.observation = np.random.rand(n_states_observe, n_states_hidden)
        self.observation /= np.sum(self.observation, axis=0, keepdims=True)

    # pi,A,Estimation la plus probable de mu
    def fit(self, sequence, iter_max=100):
        #Répétez l'étape EM
        for i in xrange(iter_max):
            params = np.hstack((self.transition.ravel(), self.observation.ravel()))

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

            #Étape M
            self.maximization(sequence, p_hidden, p_transition)

            #Vérifiez s'il a convergé
            if np.allclose(params, np.hstack((self.transition.ravel(), self.observation.ravel()))):
                break

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

        #Vers l'avant
        for i in xrange(1, len(sequence)):
            #Formule PRML(13.36)
            forward[i] = self.transition.dot(forward[i - 1]) * self.observation[sequence[i]]

        #En arrière
        for j in xrange(N - 2, -1, -1):
            #Formule PRML(13.38)
            backward[j] = (self.observation[sequence[j + 1]] * backward[j + 1]).dot(self.transition)

        #Formule PRML(13.33)Variable latente z_distribution de probabilité postérieure gamma de n(z_n)Calculer
        p_hidden = forward * backward
        p_hidden /= np.sum(p_hidden, axis=-1, keepdims=True)

        #Formule PRML(13.43)Distribution de probabilité postérieure simultanée de variables latentes consécutives xi(z_{n-1},z_n)Calculer
        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

    #Étape M
    def maximization(self, sequence, p_hidden, p_transition):
        #Formule PRML(13.18)Mettre à jour les paramètres de distribution des variables latentes initiales
        self.initial = p_hidden[0] / np.sum(p_hidden[0])

        #Formule PRML(13.19)Mettre à jour la matrice de probabilité de transition
        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)

        #Formule PRML(13.23)Mise à jour des paramètres du modèle d'observation
        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)

Les données

#Créer des données de lancer de pièces
def create_toy_data(sample_size=100):

    #Jetez des pièces en fonction du biais
    def throw_coin(bias):
        if bias == 1:
            #La table est 0.Jetez des pièces avec une probabilité de 8
            return np.random.choice(range(2), p=[0.2, 0.8])
        else:
            #Le dos est 0.Jetez des pièces avec une probabilité de 8
            return np.random.choice(range(2), p=[0.8, 0.2])

    #Biais de la première pièce
    bias = np.random.uniform() > 0.5
    coin = []
    cheats = []
    for i in xrange(sample_size):
        coin.append(throw_coin(bias))
        cheats.append(bias)

        # 0.Changer le biais des pièces avec une probabilité de 01
        bias = bias + np.random.choice(range(2), p=[0.99, 0.01])
        bias = bias % 2
    coin = np.asarray(coin)

    return coin, cheats

Code entier

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()

résultat

Si vous exécutez le code ci-dessus, vous pouvez obtenir le résultat indiqué dans la figure ci-dessous. L'axe horizontal montre le nombre de lancers de pièces et l'axe vertical montre la probabilité. La ligne bleue montre le biais des pièces réellement utilisées (deux types de pièces, l'une facile à sortir et l'autre facile à sortir). Lorsqu'elle est à 0, elle est biaisée vers l'arrière, et lorsqu'elle est égale à 1, elle est biaisée vers l'avant. Est utilisé. La ligne verte représente la probabilité que la pièce lancée à ce moment-là soit dans l'état 1 (je ne sais pas quelle pièce est dans l'état 1). Les 0 et 1 dans le graphique représentent le recto et le verso des pièces réellement observées. Cependant, il est difficile de lire si tout est affiché, il ne s'affiche donc qu'une fois toutes les deux fois.

En fonction de la valeur initiale, l'algorithme EM peut rester bloqué dans la solution optimale locale, de sorte que le résultat ci-dessus peut ne pas toujours être obtenu.

À la fin

Si les paramètres sont également modifiés en variables, cela s'intègre assez souvent dans une mauvaise solution optimale locale, il peut donc être bon d'introduire des connaissances préalables au lieu de l'estimation la plus probable et d'effectuer l'estimation de probabilité postérieure maximale.