From 'An interactive spreadsheet for teaching the forward-backward algorithm', Jason Eisner 2002¶
%pylab inline
Description du système:¶
Nombre de glaces mangées: entre 1 et 3 (valeur des observations possibles)
Nombre d'états: 2
N = 2 # nombre d'états C: Cold et H: Hot. On ne compte pas les deux états fictifs START et STOP
# Probabilités initiales de transition
a_ij = array(
[[0.8, 0.1],
[0.1, 0.8],
])
# C->C, C->H,
# H->H, H->C,
# probabilités initiales de commencer dans un des états P(C|START), P(H|START)
pi_i = [0.5, 0.5]
# probabilités initiales des observations selons les états
b_i_o = array([[0.7, 0.2, 0.1],
[0.1, 0.2, 0.7],
])
# P(1|C) P(2|C) P(3|C)
# P(1|H) P(2|H) P(3|H)
# une suite d'observations
obs = [2, 3, 3, 2, 3, 2, 3, 2, 2, 3, 1, 3, 3, 1, 1, 1, 2, 1, 1, 1, 3, 1, 2, 1, 1, 1, 2, 3, 3, 2, 3, 2, 2]
T = len(obs) # le nombre d'observations
Soit $\lambda$ le système initial défini par les probabilités d'observations $b_j(o_t)$, les probabilités de transition $a_{i,j}$ et les probabilités initiales $\pi_i$, $\alpha_t(i) = P(o_1o_2\cdots o_t, q_t=i|\lambda)$ et peut être calculé par récurrence tel que: $$\alpha_1(i) = \pi_i b_i(o_1)$$ $$\alpha_t(j) = \sum_{i=1}^N \alpha_{t-1}(i) a_{i,j} b_j(o_t)\, pour\, tout \, 1 < t \leq T$$
def calc_alpha():
""" Fonction qui calcul les alphas, sens forward, par récurence.
(probabilité de la suite d'observations partielle se terminant à l'instant t et à l'état i)
"""
# obs[t]-1 pour le décalage par rapport à la matrice
alpha = empty((T, N))
# initialisation de la récurrence
alpha[0, :] = [pi_i[i] * b_i_o[i, obs[0]-1] for i in range(N)]
# on avance petit à petit l'évaluation
for t in range(1, T):
alpha[t, :] = [b_i_o[i, obs[t]-1] * sum([a_ij[j, i] * alpha[t-1, j] for j in range(N)]) for i in range(N)]
return alpha
alpha = calc_alpha()
alpha
Soit $\lambda$ le système initial défini par les probabilités d'observations $b_j(o_t)$, les probabilités de transition $a_{i,j}$ et les probabilités initiales $\pi_i$, $\beta_t(i) = P(o_{t+1}o_{t+2}\cdots o_T, q_t=i|\lambda)$ et peut être calculé par récurrence tel que: $$ \beta_T(i) = 1$$ $$ \beta_t(i) = \sum_{j=1}^N \beta_{t+1}(j) a_{i,j} b_j(o_{t+1}) \, pour \, tout \, 0 \leq t < T$$
def calc_beta():
""" Fonction qui calcul les betas, sens backward, par récurence.
(probabilité d'observer la suite partielle ot+1ot+2...o_T, de t+1 à T est dans l'état i à l'instant t)
"""
# obs[t]-1 pour le décalage par rapport à la matrice
beta = empty((T, N))
beta[T-1, :] = 1.0 # valeur d'initialisation arbitraire, égale pour tous les états (cf rabiner89) (0.1 dans iesner)
for t in range(T - 2, -1, -1):
beta[t, :] = [sum([beta[t+1, j] * a_ij[i, j] * b_i_o[j, obs[t+1]-1] for j in range(N)]) for i in range(N)]
return beta
beta = calc_beta()
beta
Pour verifier que tout s'est bien passé on peut utiliser le fait que
$$P(O/\lambda) = \sum_{i=1}^N \alpha_T(i) = \sum_{i=1}^N \beta_1(i) = \sum_{i=1}^N \alpha_t(i) \beta_t(i) \,\, \forall t$$p_o_lambda = sum(alpha*beta, axis=1)
print(p_o_lambda)
# tout va bien ?
np.testing.assert_almost_equal(p_o_lambda, [p_o_lambda[0]]*T)
Une fois les coefficients $\alpha_t(i)$ et $\beta_t(i)$ calculés, il est possible d'estimer les coefficients $\gamma_t(i)$, la probabilité d'être dans l'état $i$ à l'instant $t$ sachant les observations et le modèle $\lambda$.
$$\gamma_t(i) = \frac{\alpha_t(i) \beta_t(i)}{P(O/\lambda)}$$def calc_gamma():
""" Fonction qui calcule les valeurs gamma en fonction des alphas et betas. """
gamma = empty((T, N))
for t in range(T):
gamma[t, :] = (alpha[t, :]*beta[t, :]) / p_o_lambda[t]
return gamma
gamma = calc_gamma()
gamma
def calc_xi():
""" Fonction qui calcule les valeurs zeta à partir des alphas et des betas."""
xi = empty((T-1, N, N))
for t in range(T-1):
for i in range(N):
xi[t, i, :] = alpha[t, i] * a_ij[i, :N] * beta[t+1, :] * b_i_o[:, obs[t+1]-1] / p_o_lambda[t]
# -1 pour faire correspondre avec la matrice d'observation, les valeurs commençant à 1 et non à 0 (nb glaces)
return xi
xi = calc_xi()
xi
Les nouvelles estimations des probabilités initiales $\pi_i$ sont données par la formule suivant:
$$\pi_i = \gamma_1(i)$$# les nouvelles valeurs des p_i
new_pi = gamma[0, :]
new_pi
La somme des $\pi_i$ doit faire 1
sum(new_pi)
Les nouvelles estimations des probabilités des transitions sont données par la formule suivante:
$$a_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)} \, pour\, i \ge 1 \, et \, j \leq N$$# les nouvelles probabilités de transition
new_a_ij = array([[sum(xi[:, i, j])/sum(gamma[:T-1, i]) for j in range(N)] for i in range(N)])
new_a_ij
Les nouvelles probabilités d'observations peuvent être calculées par la formule suivante:
$$b_i(k) = \frac{\sum_{t=1\,et\, o_t=k}^T \gamma_t(i)}{\sum_{t=1}^T \gamma_t(i)}$$new_b_i_o = array([
[sum([gamma[t, i] for t in range(T) if (obs[t]-1) == k])/sum(gamma[:, i]) for i in range(N)]
for k in range(3)]).T
new_b_i_o
Ensuite il faut itérer ce processus d'évaluation du nouveau modèle $\lambda'$ jusqu'à convergence (avec la même séquence d'observations).