Algorithme Forward-Backward - Exemple 2

by Joseph Razik, on 2019-10-18
Forward-Backward_exemple_2

Basé sur les données de 'Projet de Mathématiques: Chaînes de Markov cachées, Algorithme de Beaum-Welsh', Agen et Michot, 2007

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

Description du système:

Nombre d'observations: 2 (1: Pile, 2: Face) (valeur des observations possibles)

Nombre d'états: 3

In [2]:
N = 3  # nombre d'états. On ne compte pas les deux états fictifs START et STOP
In [3]:
# Probabilités initiales de transition
a_ij = array(
    [[0.3, 0.5, 0.2],
     [0.0, 0.3, 0.7],
     [0.0, 0.0, 1.0],
     ])
In [4]:
# probabilités initiales de commencer dans un des états
pi_i = [0.6, 0.4, 0.0]
In [5]:
# probabilités initiales des observations selons les états
b_i_o = array([[1.0, 0.0],
               [0.5, 0.5],
               [0.0, 1.0]])
In [6]:
# une suite d'observations
obs = [1, 1, 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$$

In [7]:
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
In [8]:
alpha = calc_alpha()
In [9]:
alpha
Out[9]:
array([[ 0.6   ,  0.2   ,  0.    ],
       [ 0.18  ,  0.18  ,  0.    ],
       [ 0.    ,  0.072 ,  0.162 ],
       [ 0.    ,  0.0108,  0.2124]])

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$$

In [10]:
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'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, obs[t+1]-1] for j in range(N)]) for i in range(N)]
    return beta
In [11]:
beta = calc_beta()
In [12]:
beta
Out[12]:
array([[ 0.330625,  0.124125,  0.      ],
       [ 0.4125  ,  0.8275  ,  1.      ],
       [ 0.45    ,  0.85    ,  1.      ],
       [ 1.      ,  1.      ,  1.      ]])

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$$
In [13]:
p_o_lambda = sum(alpha*beta, axis=1)
print(p_o_lambda)

[ 0.2232  0.2232  0.2232  0.2232]
In [14]:
# 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)}$$
In [15]:
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
In [16]:
gamma = calc_gamma()
In [17]:
gamma
Out[17]:
array([[ 0.88877688,  0.11122312,  0.        ],
       [ 0.33266129,  0.66733871,  0.        ],
       [ 0.        ,  0.27419355,  0.72580645],
       [ 0.        ,  0.0483871 ,  0.9516129 ]])
In [18]:
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
In [19]:
xi = calc_xi()
In [20]:
xi
Out[20]:
array([[[ 0.33266129,  0.55611559,  0.        ],
        [ 0.        ,  0.11122312,  0.        ],
        [ 0.        ,  0.        ,  0.        ]],

       [[ 0.        ,  0.17137097,  0.16129032],
        [ 0.        ,  0.10282258,  0.56451613],
        [ 0.        ,  0.        ,  0.        ]],

       [[ 0.        ,  0.        ,  0.        ],
        [ 0.        ,  0.0483871 ,  0.22580645],
        [ 0.        ,  0.        ,  0.72580645]]])

Les nouvelles estimations des probabilités initiales $\pi_i$ sont données par la formule suivant:

$$\pi_i = \gamma_1(i)$$
In [21]:
# les nouvelles valeurs des p_i
new_pi = gamma[0, :]
In [22]:
new_pi
Out[22]:
array([ 0.88877688,  0.11122312,  0.        ])

La somme des $\pi_i$ doit faire 1

In [23]:
sum(new_pi)
Out[23]:
1.0

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$$
In [24]:
# 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)])
In [25]:
new_a_ij
Out[25]:
array([[ 0.27235213,  0.59559835,  0.13204952],
       [ 0.        ,  0.24928184,  0.75071816],
       [ 0.        ,  0.        ,  1.        ]])

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)}$$
In [26]:
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(2)]).T
In [27]:
new_b_i_o
Out[27]:
array([[ 1.        ,  0.        ],
       [ 0.70704913,  0.29295087],
       [ 0.        ,  1.        ]])

Ensuite il faut itérer ce processus d'évaluation du nouveau modèle $\lambda'$ jusqu'à convergence (avec la même séquence d'observations).