Algoritmo Forward-Backward¶
En este ejemplo, planteamos una forma eficiente de programar el procedimiento de avance a partir de operaciones con arreglos de numpy. En primer lugar, llamamos a la librería de numpy, la cual utilizaremos para crear y operar con los arreglos.
import numpy as np
import pandas as pd
from ForwardBackward import forward, backward, forward_backward
Modelo Oculto de Márkov¶
Supongamos que tenemos un Modelo Oculto de Markov $HMM = (O,S,A,B,\Pi)$, donde el alfabeto de emisiones se define por:
$$S = \{D, NN, V\}$$
Y el alfabeto de observaciones está determinado como:
$$O = \{la, nina, garza, pasa\}$$
Para determinar estos símbolos, indexaremos los elementos de cada alfabeto a partir de diccionarios (en este caso hace falta únicamente indexar el alfabeto $\Sigma$):
O = {'la': 0, 'nina': 1, 'garza':2, 'pasa':3}
S = {0:'DA', 1:'N', 2:'V'}
Las probabilidades se pueden determinar de la siguiente forma:
Pi = np.array([4/6, 1/6, 1/6])
A = np.array([[1/6,1/6,1/5],[4/6,1/6,1/5],[1/6,3/6,1/5]])
B = np.array([[4/7,1/7,1/6],[1/7,2/7,1/6],[1/7,2/7,1/6],[1/7,2/7,3/6]])
pd.DataFrame(data=A, columns=S.values(), index=S.values())
| DA | N | V | |
|---|---|---|---|
| DA | 0.166667 | 0.166667 | 0.2 |
| N | 0.666667 | 0.166667 | 0.2 |
| V | 0.166667 | 0.500000 | 0.2 |
pd.DataFrame(data=Pi, columns=['Start'], index=S.values())
| Start | |
|---|---|
| DA | 0.666667 |
| N | 0.166667 |
| V | 0.166667 |
pd.DataFrame(data=B, columns=S.values(), index=O.keys())
| DA | N | V | |
|---|---|---|---|
| la | 0.571429 | 0.142857 | 0.166667 |
| nina | 0.142857 | 0.285714 | 0.166667 |
| garza | 0.142857 | 0.285714 | 0.166667 |
| pasa | 0.142857 | 0.285714 | 0.500000 |
Estas probabilidades, junto con los alfabetos de observaciones y emisiones nos dan el Modelo Oculto de Markov $HMM$. Ahora supongamos que queremos determinar la probabilidad de una cadena de observaciones $d\in\Delta^*$, definida como sigue:
obs = 'la nina pasa la garza'
obs = obs.split()
print(obs)
['la', 'nina', 'pasa', 'la', 'garza']
Forward Procedure¶
Utilizaremos el procedimiento de avance para determinar dicha probabilidad. Realizaremos los pasos de inicialización, inducción y terminación (en cada paso, guardamos las variables en una matriz):
(1) Inicialización. La inicialización está dada por el almacenamiento de las probabilidades iniciales en la variable de avance:
$$\alpha_i(0) := \pi_i$$
(2) Inducción. Los siguientes pasos consistirán en ir actualizando la variable de avance a partir de los diferentes estados de la cadena:
$$\alpha_i(t+1) = \sum_{j=1}^N p(o_t|s_i) p(s_i|s_j)\alpha_j(t)$$
En este caso, utilizaremos operaciones entre matrices y vectores, pues es evidente que (tomando los elementos del modelo $A,B,\Pi$):
$$\alpha_i(t+1) = B_{t,\cdot} \odot (A_{i,\cdot} \alpha_{\cdot}(t))$$
donde $B_{t,\cdot}$ representa el vector renglón que corresponde a la observación en el estado $t$ y $\alpha_{\cdot}(t)$ es el vector que contiene todas las variables $\alpha_i(t)$. Aquí, $A_{i,\cdot}$ representa el i-ésimo vector renglón de $A$ y $\odot$ representa el producto de Hadamard. En general, podemos actualizar todo el vector $\alpha(t+1)$ (cuyas entradas son $\alpha_i(t+1)$) tomando la matriz completa $A$. De tal forma que:
$$\alpha(t+1) = B_{t,\cdot} \odot (A \alpha_{\cdot}(t)) $$
(3) Terminación. Para la terminación, tenemos que:
$$p(o^{(1)}...o^{(T)}) = \sum_{i=1}^N \alpha_i(T)$$
Pero ya que hemos guardado cada variable $\alpha_i$ como la entrada de un vector, basta sumar las entradas de este vector para obtener esta probabilidad:
alpha = forward(obs,(O,S,A,B,Pi))
pd.DataFrame(data=alpha, columns=S.values(), index=[0] + obs)
simbolo: la - probabilidad 0.20859788359788356 simbolo: la nina - probabilidad 0.03951877047115142 simbolo: la nina pasa - probabilidad 0.011220368570708704 simbolo: la nina pasa la - probabilidad 0.0020895512786453284 simbolo: la nina pasa la garza - probabilidad 0.00041117213569291135
| DA | N | V | |
|---|---|---|---|
| 0 | 0.666667 | 0.166667 | 0.166667 |
| la | 0.098413 | 0.072222 | 0.037963 |
| nina | 0.005147 | 0.024354 | 0.010018 |
| pasa | 0.000989 | 0.002713 | 0.007519 |
| la | 0.001212 | 0.000374 | 0.000504 |
| garza | 0.000052 | 0.000277 | 0.000082 |
print('Probabilidad de la cadena',obs, 'es:', alpha[-1].sum(0) )
Probabilidad de la cadena ['la', 'nina', 'pasa', 'la', 'garza'] es: 0.00041117213569291135
Backward Procedure¶
El procedimiento de retroceso o backward procedure es similar al procedimiento de avance, pero la cadena se recorre en sentido inverso, es decir, de derecha a izquierda. Aquí presentamos una implementación simple de este procedimiento utilizando las propiedades de operaciones entre arreglos en numpy.
Ahora es cuando iniciamos con el procedimiento de retroceso. El paso de inducción consiste en determinar las variables $\beta_i(T) = 1$. En este caso, utilizamos un vector, que podemos denotar como $\beta_\cdot(T)$, el cual contendrá en un primer momento sólo unos.
En el paso inductivo tenemos que obtener las variables:
$$\beta_j(t) = \sum_{i=1}^N p(o^{(t+1)} | s_i^{(t+1)}) p(s_i^{(t+1)}|s_{j}^{(t)}) \beta_i(t+1)$$
Ya que estamos trabajando con arreglos de numpy, podemos ver que la suma y el producto sobre los $i$-ésimos elementos determinan un producto punto. Así, $B_{(t+1),\cdot}$ representa el vector columna correspondiente al elemento de la cadena de observaciones en el estado $t+1$, $A_{\cdot, j}$ a la $j$-ésima columna de la matriz de transiciones $A$ y $\hat{\beta}(t+1)$ al vector de las variables $\beta$ en el estado subsecuente. Así, podemos expresar la ecuación anterior como:
$$\beta_j(t) = [\hat{\beta}(t+1) \odot B_{(t+1),\cdot}] \cdot A_{\cdot, j}$$
Donde $\cdot$ es el producto punto, y $\odot$ el producto de Hadarmard. En términos generales, podemos actualizar todas las variables $\beta$ expresándolo como un vector $\hat{\beta}(t)$ y de esta forma, tenemos que realizar la operación:
$$\hat{\beta(t)} = A^T [\hat{\beta}(t+1) \odot B_{(t+1),\cdot}]$$
Donde el producto punto se hace con toda la matriz de transiciones $A$.
Finalmente, el paso de terminación consiste en obtener la probabilidad de la cadena de la siguiente forma:
$$p(o^{(1)}...o^{(T)}) = \sum_{j=1}^N \pi_j \beta_j(0)$$
Donde $\pi_j$ es la $j$-ésima probabilidad inicial. De igual forma, podemos aprovechar la estructura de los arreglos y expresar la suma de los productos por medio de un producto punto con el vector de probabilidades iniciales $\Pi$ y el vector de variables $\beta(0)$. Esto es:
$$p(o^{(1)}...o^{(T)}) = \Pi \cdot \hat{\beta}(0)$$
beta = backward(obs,(O,S,A,B,Pi))
pd.DataFrame(data=beta, columns=S.values(), index=[0] + obs)
simbolo: garza - probabilidad: 0.20701058201058198 simbolo: la garza - probabilidad: 0.03951877047115142 simbolo: pasa la garza - probabilidad: 0.010412587433562487 simbolo: nina pasa la garza - probabilidad: 0.002226236896034973 simbolo: la nina pasa la garza - probabilidad: 0.0004111721356929113
| DA | N | V | |
|---|---|---|---|
| 0 | 0.000422 | 0.000390 | 0.000387 |
| la | 0.002688 | 0.001393 | 0.001213 |
| nina | 0.010966 | 0.011738 | 0.006875 |
| pasa | 0.041100 | 0.036659 | 0.036054 |
| la | 0.242063 | 0.154762 | 0.119048 |
| garza | 1.000000 | 1.000000 | 1.000000 |
print('La probabilidad de la observación es:', np.dot(Pi,beta[0]) )
La probabilidad de la observación es: 0.0004111721356929113
Forward-Backward Procedure¶
El procedimiento de avance-retroceso puede usarse para calcular la probabilidad de una cadena; en particular, tenemos que:
$$p(o^{(1)}...o^{(T)}) = \sum_{i=1}^N \alpha_i(t) \beta_i(t)$$
En este sentido para cada $t$ (de $1$ a $T$) podremos utilizar las variables en ese estado para calcular la probabilidad.
El avance-retroceso nos da la probabilidad conjunta de la cadena y el símbolo en un estado $t$:
$$p(o^{(1)}...o^{(T)}, s_i^{(t)}) = \alpha_i(t) \beta_i(t)$$
Más aún, también nos puede dar las probabilidades condicionales de la siguiente forma:
$$p(s_i^{(t)}|o^{(1)}...o^{(T)}) = \frac{\alpha_i(t) \beta_i(t)}{\sum_i \alpha_i(t) \beta_i(t)}$$
result = forward_backward(obs, (O,S,A,B,Pi), normalize=True)
#Probabilidad condicional
pd.DataFrame(data=result, columns=list(S.values()), index=['0']+obs)
simbolo: la - probabilidad 0.20859788359788356 simbolo: la nina - probabilidad 0.03951877047115142 simbolo: la nina pasa - probabilidad 0.011220368570708704 simbolo: la nina pasa la - probabilidad 0.0020895512786453284 simbolo: la nina pasa la garza - probabilidad 0.00041117213569291135 simbolo: garza - probabilidad: 0.20701058201058198 simbolo: la garza - probabilidad: 0.03951877047115142 simbolo: pasa la garza - probabilidad: 0.010412587433562487 simbolo: nina pasa la garza - probabilidad: 0.002226236896034973 simbolo: la nina pasa la garza - probabilidad: 0.0004111721356929113
| DA | N | V | |
|---|---|---|---|
| 0 | 0.684777 | 0.158185 | 0.157038 |
| la | 0.643323 | 0.244665 | 0.112012 |
| nina | 0.137279 | 0.695233 | 0.167488 |
| pasa | 0.098821 | 0.241848 | 0.659331 |
| la | 0.713422 | 0.140610 | 0.145968 |
| garza | 0.126837 | 0.674711 | 0.198452 |
Podemos entonces predecir la probabilidad de una etiqueta y obtener la etiqueta adecuada para cada estado de la cadena de observaciones. La predicción se hace como:
\begin{align} \hat{s}^{(t)} &= \arg\max_i p(s_i^{(t)}|o^{(1)}...o^{(T)}) = \arg\max_i \frac{\alpha_i(t) \beta_i(t)}{\sum_i \alpha_i(t) \beta_i(t)} \\ &= \arg\max_i p(o^{(1)}...o^{(T)}, s_i^{(t)}) = \arg\max_i \alpha_i(t) \beta_i(t) \end{align}
Como ya hemos guardado las variables $\alpha$ y $\beta$ basta multiplicarlas y obtener el argumento que maximiza el producto.
#Predicción en base a forward-backward
pred = [S[t] for t in np.argmax(result, axis=1)]
#Obtención de la secuencia etiquetada
tag_obs = list(zip(obs,pred[1:]))
print(tag_obs)
[('la', 'DA'), ('nina', 'N'), ('pasa', 'V'), ('la', 'DA'), ('garza', 'N')]