Implémentation de l'estimation des paramètres HMM en python. Comme manuel, j'ai utilisé ["Reconnaissance de formes continue / facile à comprendre"](https://www.amazon.co.jp/ Suite / Reconnaissance de formes facile à comprendre - Introduction à l'apprentissage sans enseignant - Ishii-Kenichiro / dp / 427421530X).
Structure de cet article
Hidden Markov Model
Une caractéristique dans laquelle l'état actuel est déterminé de manière probabiliste en fonction de l'état il y a une fois est appelée ** propriété de Markov **. Le processus de satisfaction de la propriété Markov est appelé le ** processus de Markov **. Des exemples de tels événements comprennent les cours des actions, les signaux vocaux et les langues. Dans le modèle de Markov, un modèle dans lequel seule la série de symboles de sortie peut être observée et la série d'états ne peut pas être observée. Il s'appelle ** Modèle de Markov caché **.
HMM sera expliqué en utilisant le lancer de dés comme exemple. Pensez à lancer $ c $ seed dice avec $ m $ seeds $ v_1, v_2, \ cdots, v_m $ $ \ omega_1, \ omega_2, \ cdots, \ omega_c $ $ n $ fois, Le type de dés lancés à l'heure $ t $ est représenté par $ s_t $, et le nombre de dés lancés et observés est représenté par $ x_t $. $ s_t $ et $ x_t $ sont exprimés comme suit.
\begin{eqnarray}
\left\{\begin{array}{ll}
s_t \in \bigl\{\omega_1, \omega_2, \cdots, \omega_c\bigr\} \\
x_t \in \bigl\{v_1, v_2, \cdots, v_m\bigr\} \\
\end{array} \right.
\end{eqnarray}
Le dé $ s_t $ sorti au $ t $ ème temps est déterminé stochastiquement par les dés $ s_ {t-1} $ sortis au temps $ (t -1) $. La probabilité est représentée par $ a (s_ {t -1}, s_t) $. Lorsque le dé $ (t -1) $ sorti est $ \ omega_i $, la probabilité que le dé sorti $ t $ soit $ \ omega_j $ est Il s'exprime comme suit.
a(\omega_i, \omega_j) = P(s_t = \omega_j \ | \ s_{t - 1} = \omega_i) \ \ \ \ \ (i, j = 1, 2, \cdots, c) \tag{1}
La probabilité que $ x_t $ soit observé en lançant les dés $ s_t $ au $ t $ ème temps est représentée par $ b (s_t, x_t) $. La probabilité que $ v_k $ soit observé en lançant les dés $ \ omega_j $ est exprimée comme suit.
b(\omega_j, v_k) = P(x_t = v_k \ | \ s_t = \omega_j) \ \ \ \ \ (j = 1, 2, \cdots, c) \ \ (k = 1, 2, \cdots, m) \tag{2}
Pour plus de simplicité, $ a (\ omega_i, \ omega_j) $ est abrégé en $ a_ {ij} $, $ b (\ omega_j, v_k) $ est abrégé en $ b_ {jk} $, et $ a_ {ij} $ est abrégé en *. * La probabilité de transition ** et $ b_ {jk} $ sont appelées ** probabilité de sortie **. Avec HMM, vous pouvez obtenir des informations sur les ** yeux de sortie ** (= résultats d'observation) en lançant un dé, L'information sur le ** type de dés ** (= état) qui produit le jet ne peut pas être obtenue.
Les symboles pour exprimer le modèle de Markov sont résumés ci-dessous.
Lorsque les paramètres $ \ boldsymbol A, \ boldsymbol B, \ boldsymbol \ rho $ sont connus Trouver la probabilité $ P (\ boldsymbol x) $ d'obtenir le résultat de l'observation $ \ boldsymbol x $ s'appelle ** evaluation **.
\begin{align}
P(\boldsymbol x) &= \sum_{\boldsymbol s} P(\boldsymbol x, \boldsymbol s) \tag{3} \\
&= \sum_{s_1} \sum_{s_2} \cdots \sum_{s_n} P(\boldsymbol x, s_1s_2 \cdots s_n) \\
P(\boldsymbol x, \boldsymbol s) &= \prod_{t = 1}^n a(s_{t - 1}, s_t)b(s_t, x_t) \tag{4}
\end{align}
Si $ \ boldsymbol x $ et $ \ boldsymbol s $ peuvent être observés en même temps, $ P (\ boldsymbol x, \ boldsymbol s) $ peut être facilement calculé. D'un autre côté, si seulement $ \ boldsymbol x $ peut être observé, la quantité de calcul de marginalisation pour tous les états possibles $ \ boldsymbol s $ sera énorme. Nous allons introduire ** algorithme avant ** et ** algorithme arrière ** comme algorithmes pour un calcul efficace.
forward algorithm
La probabilité que le résultat de l'observation $ \ boldsymbol x $ soit obtenu et que le dé $ \ omega_i $ soit sorti au $ t $ ème temps est la suivante.
\begin{align}
& P(\boldsymbol x, s_t = \omega_i) \\
&= P(x_1x_2 \cdots x_t, s_t = \omega_i) \cdot P(x_{t + 1}x_{t + 2} \cdots x_n \ | \ s_t = \omega_i) \\
&= \alpha_t(i) \cdot \beta_t(i) \tag{5}
\end{align}
$ \ Alpha_t (i) $ dans la formule $ (5) $ est la probabilité que $ x_1x_2 \ cdots x_t $ soit obtenu et que le dé $ \ omega_i $ soit sorti au $ t $ ème temps. $ \ Alpha_t (i) $ peut être calculé récursivement comme suit.
\begin{align}
& \alpha_t(j) = \Biggl[\sum_{i = 1}^c \alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{6} \\
& (t = 2, 3, \cdots, n) \ \ (j = 1, 2, \cdots, c)
\end{align}
Cependant, $ \ alpha_1 (i) $ est comme suit.
\begin{align}
\alpha_1(i) &= P(x_1, s_1 = \omega_i) \\
&= \rho_i b(\omega_i, x_1) \ \ \ \ \ (i = 1, 2, \cdots, c) \tag{7}
\end{align}
L'équation suivante est obtenue en définissant $ t = n $ dans l'équation $ (6) $.
\alpha_n(i) = P(x_1x_2 \cdots x_n, s_n = \omega_i) \tag{8}
L'équation suivante est obtenue en marginalisant l'équation $ (8) $ par rapport à $ s_n $.
\begin{align}
P(\boldsymbol x) &= P(x_1x_2 \cdots x_n) \\
&= \sum_{i = 1}^c P(x_1x_2 \cdots x_n, s_n = \omega_i) \\
&= \sum_{i = 1}^c \alpha_n(i) \tag{9}
\end{align}
Sur la base de ce qui précède, la procédure de traitement de l'algorithme de transfert est indiquée ci-dessous. < forward algorithm >
backward algorithm
$ \ Beta_t (i) $ dans la formule $ (5) $ est soumis à la condition que le dé $ \ omega_i $ soit sorti au $ t $ ème temps. C'est la probabilité que le résultat de l'observation après cela soit $ x_ {t + 1} x_ {t + 2} \ cdots x_n $. $ \ Beta_t (i) $ peut être calculé récursivement comme suit.
\begin{align}
& \beta_t(i) = \sum_{j = 1}^c a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \tag{10} \\
& (t = 1, 2, \cdots, (n - 1)) \ \ (i = 1, 2, \cdots, c)
\end{align}
Cependant, $ \ beta_n (i) $ est comme suit.
\beta_n(i) = 1 \ \ \ \ \ (i = 1, 2, \cdots, c) \tag{11}
L'équation suivante est obtenue en définissant $ t = 1 $ dans l'équation $ (10) $.
\beta_1(i) = P(x_2x_3 \cdots x_n \ | \ s_1 = \omega_i) \tag{12}
A partir de la formule $ (12) $, $ P (x_1x_2 \ cdots x_n, s_1 = \ omega_i) $ peut être exprimé comme suit.
\begin{align}
& P(x_1x_2 \cdots x_n, s_1 = \omega_i) \\
&= P(s_1 = \omega_i)P(x_1 \ | \ s_1 = \omega_i)P(x_2x_3 \cdots x_n \ | \ s_1 = \omega_i) \\
&= \rho_i b(\omega_i, x_1) \beta_1(i) \tag{13}
\end{align}
L'équation suivante est obtenue en marginalisant l'équation $ (13) $ par rapport à $ s_n $.
\begin{align}
P(\boldsymbol x) &= P(x_1x_2 \cdots x_n) \\
&= \sum_{i = 1}^c P(x_1x_2 \cdots x_n, s_1 = \omega_i) \\
&= \sum_{i = 1}^c \rho_i b(\omega_i, x_1) \beta_1(i) \tag{14}
\end{align}
Sur la base de ce qui précède, la procédure de traitement de l'algorithme arrière est illustrée ci-dessous. < backward algorithm >
Lorsque les paramètres $ \ boldsymbol A, \ boldsymbol B, \ boldsymbol \ rho $ sont inconnus L'obtention d'un paramètre inconnu à partir du résultat de l'observation $ \ boldsymbol x $ s'appelle ** estimation **. Nous définissons $ \ gamma_t (i) $ et $ \ xi_t (i, j) $ comme suit en préparation pour dériver la solution de l'estimation.
\begin{align}
\gamma_t(i) &\overset{\rm{def}}{=} P(s_t = \omega_i \ | \ \boldsymbol x) \tag{15} \\
\xi_t(i) &\overset{\rm{def}}{=} P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \tag{16}
\end{align}
De plus, la relation suivante est valable pour $ \ gamma_t (i) $ et $ \ xi_t (i, j) $.
\begin{align}
\sum_{j = 1}^c \xi_t(i, j) &= \sum_{j = 1}^c P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \\
&= P(s_t = \omega_i \ | \ \boldsymbol x) \\
&= \gamma_t(i) \tag{17}
\end{align}
Ici, $ P (\ boldsymbol x, s_t = \ omega_i, s_ {t + 1} = \ omega_j) $ est exprimé comme suit.
\begin{align}
&P(\boldsymbol x, s_t = \omega_i, s_{t + 1} = \omega_j) \\
&= P(x_1 \cdots x_t, s_t = \omega_i) \cdot P(s_{t + 1} = \omega_j \ | \ s_t = \omega_i) \cdot P(x_{t + 1} \ | \ s_{t + 1} = \omega_j) \cdot P(x_{t + 2} \cdots x_n \ | \ s_{t + 1} = \omega_j) \\
&= \alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \tag{18}
\end{align}
De ce qui précède, l'équation suivante peut être dérivée.
\begin{align}
\gamma_t(i) &= P(s_t = \omega_i \ | \ \boldsymbol x) \\
&= P(\boldsymbol x, s_t = \omega_i) \bigl/ P(\boldsymbol x) \\
&= \alpha_t(i)\beta_t(i) \Bigl/ \sum_{i = 1}^c \alpha_n(i) \tag{19} \\
\xi_t(i, j) &= P(s_t = \omega_i, s_{t + 1} = \omega_j \ | \ \boldsymbol x) \\
&= P(\boldsymbol x, s_t = \omega_i, s_{t + 1} = \omega_j) \bigl/ P(\boldsymbol x) \\
&= \alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\beta_{t + 1}(j) \Bigl/ \sum_{i = 1}^c \alpha_n(i) \tag{20}
\end{align}
Baum-Welch algorithm
Sur la base des résultats ci-dessus, les paramètres de HMM sont estimés en utilisant l'estimation la plus probable. La méthode d'estimation expliquée cette fois est appelée ** algorithme de Baum-Welch **. Les paramètres $ \ boldsymbol A, \ boldsymbol B, \ boldsymbol \ rho $ sont résumés comme suit.
\boldsymbol \theta = (\boldsymbol A, \boldsymbol B, \boldsymbol \rho) \tag{21}
Appliquez l'algorithme EM et maximisez la fonction $ Q $ suivante.
\begin{align}
Q(\boldsymbol \theta^0, \boldsymbol \theta) &= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \log P(\boldsymbol x, \boldsymbol s; \boldsymbol \theta) \\
&= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \bigl\{\log P(\boldsymbol s; \boldsymbol \theta) + \log P(\boldsymbol x \ | \ \boldsymbol s; \boldsymbol \theta)\bigr\} \\
&= \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \bigl\{\log P(s_1) + \sum_{t = 1}^{n - 1} \log a(s_t, s_{t + 1}) + \sum_{t = 1}^n \log b(s_t, x_t) \bigr\} \tag{22}
\end{align}
Ici, il est défini comme suit.
\begin{align}
Q(\boldsymbol \theta^0, \boldsymbol \rho) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \log P(s_1) \tag{23} \\
Q(\boldsymbol \theta^0, \boldsymbol A) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \sum_{t = 1}^{n - 1} \log a(s_t, s_{t + 1}) \tag{24} \\
Q(\boldsymbol \theta^0, \boldsymbol B) &\overset{\rm{def}}{=} \sum_{\boldsymbol s} P(\boldsymbol s \ | \ \boldsymbol x; \boldsymbol \theta^0) \sum_{t = 1}^n \log b(s_t, x_t) \tag{25}
\end{align}
Maximisez les équations $ (23), (24), (25) $ pour $ \ boldsymbol \ rho, \ boldsymbol A, \ boldsymbol B $, respectivement. Les paramètres estimés sont les suivants. (La dérivation est omise)
\begin{align}
\hat{a}_{ij} &= \sum_{t = 1}^{n - 1} \xi_t(i, j) \Bigl/ \sum_{t = 1}^{n - 1} \gamma_t(i) \tag{26} \\
\hat{b}_{jk} &= \sum_{t = 1}^n \delta(x_t, v_k) \gamma_t(i) \Bigl/ \sum_{t = 1}^n \gamma_t(i) \tag{27} \\
\hat{\rho}_i &= \gamma_1(i) \tag{28}
\end{align}
La procédure de traitement de l'algorithme Baum-Welch est illustrée ci-dessous. < Baum-Welch algorithm >
L'algorithme de Baum-Welch utilise l'algorithme avant et l'algorithme arrière pour calculer $ \ alpha_t (i) et \ beta_t (i) $. Si la série devient longue, un sous-débordement se produira et il ne sera pas possible de calculer $ \ alpha_t (i) et \ beta_t (i) $. Nous allons introduire ** la mise à l'échelle ** comme contre-mesure contre le sous-débit. Introduisez le nouveau $ \ hat \ alpha_t (i), c_t $.
\begin{align}
& \hat \alpha_t(i) \overset{\rm{def}}{=} P(s_t = \omega_i \ | \ x_1 \cdots x_t) = \frac{P(x_1 \cdots x_t, s_t = \omega_i)}{P(x1 \cdots x_t)} \tag{29} \\
& c_t = P(x_t \ | \ x_1 \cdots x_{t - 1}) \tag{30}
\end{align}
Ici, $ P (x_1 \ cdots x_t) $ est exprimé comme suit.
\begin{align}
P(x_1 \cdots x_t) &= P(x_1) \cdot P(x_2 \ | \ x_1) \cdot \cdots \cdot P(x_{t - 1} \ | \ x_1 \cdots x_{t - 2}) \cdot P(x_t \ | \ x_1 \cdots x_{t - 1}) \\
&= c_1 \cdot c_2 \cdot \cdots \cdot c_{t - 1} \cdot c_t \\
&= \prod_{m = 1}^t c_m \tag{31}
\end{align}
À partir de là, la relation suivante tient.
\begin{align}
\hat \alpha_t(i) &= \frac{\alpha_t(i)}{P(x_1 \cdots x_t)} \\
&= \frac{\alpha_t(i)}{\prod_{m = 1}^t c_m} \\
&= \frac{\alpha_t(i)}{c_t \cdot \prod_{m = 1}^{t - 1} c_m} \\
&= \frac{\alpha_t(i)}{c_t \cdot P(x_1 \cdots x_{t - 1})} \tag{32}
\end{align}
L'expression $ (6) $ peut être réécrite comme suit.
\begin{align}
\alpha_t(j) &= \Biggl[\sum_{i = 1}^c \alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
&= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) P(x_1 \cdots x_{t - 1}) a_{ij} \Biggr]b(\omega_j, x_t) \\
\frac{\alpha_t(j)}{P(x_1 \cdots x_{t - 1})} &= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
c_t \hat\alpha_t(j)&= \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{33}
\end{align}
A partir de la formule $ (33) $, $ c_t $ peut être calculé comme suit.
\begin{align}
\sum_{j = 1}^c c_t \hat\alpha_t(j) &= \sum_{j = 1}^c \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \\
c_t &= \sum_{j = 1}^c \Biggl[\sum_{i = 1}^c \hat\alpha_{t - 1}(i) a_{ij} \Biggr]b(\omega_j, x_t) \tag{34}
\end{align}
De plus, $ \ hat \ beta_t (i) $ est exprimé comme suit en utilisant $ c_t $ obtenu par l'algorithme de transfert.
c_{t + 1} \hat\beta_t(i) = \sum_{j = 1}^c a_{ij}b(\omega_j, x_{t + 1})\hat\beta_{t + 1}(j) \tag{35}
Avec $ \ hat \ alpha_t (i), \ hat \ beta_t (i) $, $ \ gamma_t (i), \ xi_t (i, j) $ peuvent être réécrits comme suit.
\begin{align}
\gamma_t(i) &= \hat\alpha_t(i)\hat\beta_t(i) \tag{36} \\
\xi_t(i) &= \frac{\hat\alpha_t(i)a_{ij}b(\omega_j, x_{t + 1})\hat\beta_{t + 1}(j)}{c_{t + 1}} \tag{37}
\end{align}
Les paramètres sont mis à jour en remplaçant les équations $ (36), (37) $ dans les équations $ (26), (27), (28) $.
Implémentation de l'estimation des paramètres à l'aide de l'algorithme Baum-Welch en python.
Le code d'implémentation peut être trouvé ici.
Basé sur le lancer de dés, le nombre d'états est de 3 $ espèces ($ c = 3 $) de $ \ omega_1, \ omega_2, \ omega_3 $,
Le symbole de sortie est l'espèce $ 2 $ ($ m = 2 $) de $ v_1 et v_2 $.
Cela équivaut à sélectionner et à lancer différents dés d'espèces $ 3 $ et à observer des lancers impairs ($ v_1
Les vrais paramètres sont indiqués ci-dessous.
\boldsymbol A = \begin{pmatrix}
0.1 & 0.7 & 0.2 \\
0.2 & 0.1 & 0.7 \\
0.7 & 0.2 & 0.1
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.9 & 0.1 \\
0.6 & 0.4 \\
0.1 & 0.9
\end{pmatrix}, \ \ \
\boldsymbol \rho = \begin{pmatrix}
1.0 & 0.0 & 0.0
\end{pmatrix}
Les valeurs initiales d'estimation ont été fixées comme suit.
\boldsymbol A = \begin{pmatrix}
0.15 & 0.60 & 0.25 \\
0.25 & 0.15 & 0.60 \\
0.60 & 0.25 & 0.15
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.5 & 0.5 \\
0.5 & 0.5 \\
0.5 & 0.5
\end{pmatrix}, \ \ \
\boldsymbol \rho = \begin{pmatrix}
1.0 & 0.0 & 0.0
\end{pmatrix}
Le résultat de l'estimation est le suivant.
\boldsymbol A = \begin{pmatrix}
0.11918790 & 0.63863200 & 0.24218010 \\
0.13359797 & 0.07275209 & 0.79364994 \\
0.72917405 & 0.18715812 & 0.08366783
\end{pmatrix}, \ \ \
\boldsymbol B = \begin{pmatrix}
0.92503797 & 0.07496203 \\
0.56298958 & 0.43701042 \\
0.14059588 & 0.85940412
\end{pmatrix}
Le processus d'augmentation de la vraisemblance logarithmique est le suivant.
Enfin, le processus de convergence de $ b_ {11}, b_ {21}, b_ {31} $ est illustré ci-dessous. La fine ligne noire représente la vraie valeur.
J'ai pu implémenter l'estimation des paramètres HMM en python. Si vous avez des erreurs, veuillez les signaler.
Recommended Posts