%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
# une suite d'observations
obs_1 = [11, 9, 10, 18, 19, 18, 14, 15, 14]
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
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])
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
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
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$$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
def calc_xi(alpha, beta, p_o_lambda, a_ij, 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
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).
def iteration(nb_iter, obs):
global mu_i, sigma_i
a_ij = array(
[[0.5, 0.5, 0.0],
[0.0, 0.5, 0.5],
[0.0, 0.0, 1.0],
])
pi_i = [1.0, 0.0, 0.0]
mu_i = None # [0.0, 0.0, 0.0]
sigma2_i = None # [1.0, 1.0, 1.0]
N = 3
T = len(obs)
obs = array(obs)
for iterate in range(nb_iter):
alpha = calc_alpha(pi_i, a_ij, obs, N)
beta = calc_beta(a_ij, obs, N)
p_o_lambda = sum(alpha*beta, axis=1)
gamma = calc_gamma(alpha, beta, p_o_lambda, T, N)
xi = calc_xi(alpha, beta, p_o_lambda, a_ij, obs, N)
# nouvelles estimations
new_pi = gamma[0, :]
new_a_ij = array([[sum(xi[:, i, j])/sum(gamma[:T-1, i]) for j in range(N)] for i in range(N)])
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_pi: ', new_pi)
print('new_a_ij', new_a_ij)
print('new_b_i_o', new_mu_i, new_sigma2_i)
iteration(1, obs_1)
iteration(5, obs_1)
Avec deux autres séquences d'observations
- 10, 11, 19, 17, 19, 19, 13, 14
- 9, 17, 18, 17, 17, 13, 15
obs_2 = [10, 11, 19, 17, 19, 19, 13, 14]
obs_3 = [9, 17, 18, 17, 17, 13, 15]
iteration(1, obs_2)
iteration(1, obs_3)