%pylab inline
Description du système:¶
Nombre d'observations: 21 possibles (entre 0 et 20, valeurs des observations possibles)
Nombre d'états: 3
Séquence d'bservations:
- 11, 9, 10, 18, 19, 18, 14, 15, 14
N = 3 # nombre d'états. On ne compte pas les deux états fictifs START et STOP
# Probabilités initiales de transition, hmm gauche-droite strict
a_ij = array(
[[0.5, 0.5, 0.0],
[0.0, 0.5, 0.5],
[0.0, 0.0, 1.0],
])
# P(C|C), P(H|C)
# P(C|H), P(H|H)
# plus de P(STOP|C), P(STOP|H)
# probabilités initiales de commencer dans un des états P(C|START), P(H|START) (et P(STOP|START) mais égal à 0)
pi_i = [1.0, 0.0, 0.0]
La différence entre les HMM à observation discrète et les HMM continus, est le fait que l'ensemble des observations n'est pas forcément limité à un ensemble de valeurs fixe et connu. On rencontre habituellement une modélisation par Gaussienne dans le cas continu. Ainsi, la probabilité d'une observersation est dorénavant donnée par sa loi de densité de probabilité:
$$N(\mu, \sigma^2) = \frac{1}{\sqrt{2\pi|\sigma|}} \exp^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}} $$# probabilités initiales des observations selons les états
mu_i = None # [0.0, 0.0, 0.0]
sigma2_i = None # [1.0, 1.0, 1.0]
def b_i_o(i, t):
if mu_i is None:
return 1.0
else:
return (1/sqrt(2*pi*abs(sigma2_i[i]))) * exp(-0.5*(obs[t]-mu_i[i])**2/sigma2_i[i])
# une suite d'observations
obs = array([11, 9, 10, 18, 19, 18, 14, 15, 14])
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(pi_i, a_ij, obs, N):
""" 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
T = len(obs)
alpha = empty((T, N))
# initialisation de la récurrence
alpha[0, :] = [pi_i[i] * b_i_o(i, 0) for i in range(N)]
# on avance petit à petit l'évaluation
for t in range(1, T):
alpha[t, :] = [b_i_o(i, t) * sum([a_ij[j, i] * alpha[t-1, j] for j in range(N)]) for i in range(N)]
return alpha
alpha = calc_alpha(pi_i, a_ij, obs, N)
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(a_ij, obs, N):
""" 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
T = len(obs)
beta = empty((T, N))
beta[T-1, :] = 1.0 # valeur d'nitialisation arbitraire mais égale pour tous les états (cf rabiner89)
for t in range(T - 2, -1, -1):
beta[t, :] = [sum([beta[t+1, j] * a_ij[i, j] * b_i_o(j, t+1) for j in range(N)]) for i in range(N)]
return beta
beta = calc_beta(a_ij, obs, N)
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(alpha, beta, p_o_lambda, T, N):
""" 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(alpha, beta, p_o_lambda, T, N)
gamma
def calc_xi(alpha, beta, obs, N):
""" Fonction qui calcule les valeurs zeta à partir des alphas et des betas."""
T = len(obs)
xi = empty((T-1, N, N))
for t in range(T-1):
for i in range(N):
for j in range(N):
xi[t, i, j] = alpha[t, i] * a_ij[i, j] * beta[t+1, j] * b_i_o(j, t+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(alpha, beta, obs, N)
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 dépendent des paramètres définissant le modèle gaussien: les moyennes et l'écart-type qui peuvent être calculées par la formule suivante:
$$\mu_i = \frac{\sum_{t=1}^T \gamma_t(i) o_t}{\sum_{t=1}^T \gamma_t(i)}$$$$\sigma^2_i = \frac{\sum_{t=1}^T \gamma_t(i) (o_t - \mu_i) (o_t - \mu_i)'}{\sum_{t=1}^T \gamma_t(i)} $$Bien remarquer maintenant dans l'équation de ré-estimation de $\mu_i$ la multiplication par l'apparition de la valeur elle-même. Dans le cas discret, il y avait uniquement observation ou non observation (0 ou 1).
Dans l'équation de $\sigma^2$, le caractère "'" signifie "transposé" (dans le cas multidimensionnel).
# estimer b_i_o revient à estimer mu et sigma
new_mu_i = [ sum(gamma[:, i]*obs[:]) / sum(gamma[:, i]) for i in range(N)]
new_sigma2_i = [ sum(gamma[:, i] * (obs[:] - new_mu_i[i])**2) / sum(gamma[:, i]) for i in range(N)]
print(new_mu_i)
print(new_sigma2_i)
Ensuite il faut itérer ce processus d'évaluation du nouveau modèle $\lambda'$ jusqu'à convergence (avec la même séquence d'observations).