Cours

Partie obligatoire

TME 7: MMC et reconnaissance de lettres

La base de données est la même que dans la séance précédente:
TME6_lettres.pkl

Les mécanismes de discrétisation et de tracé sont donc les mêmes (cf semaine 6)

import numpy as np
import matplotlib.pyplot as plt
# truc pour un affichage plus convivial des matrices numpy
np.set_printoptions(precision=2, linewidth=320)
plt.close('all')

data = pkl.load(file("ressources/lettres.pkl","rb"))
X = np.array(data.get('letters'))
Y = np.array(data.get('labels'))

nCl = 26

Passage au MMC

Le but de la séance est de mettre en oeuvre un modèle de Markov caché (MMC). Pour cela, nous allons développer les algorithmes vus en TD.

Apprentissage d'un modèle connaissant les états

Hypothèse Gauche-Droite

On fait l'hypothèse que les états sont connus... Alors que ce n'est pas le cas. Mais il existe des stratégies simples (et parfois efficaces) pour attribuer arbitrairement des états sur les chaines. La plus classique est l'hypothèse gauche-droite, qui est bien adaptée aux signaux courts et non périodiques:

  • On définit le nombre d'états N
  • On découpe chaque série d'observations en N portions à peu près égales
  • On affecte l'état 0 au début, puis on incrémente jusqu'à l'état N pour la dernière portion de la chaine

Sur un exemple:

  X0 = [ 1.  9.  8.  8.  8.  8.  8.  9.  3.  4.  5.  6.  6.  6.  7.  7.  8.  9.  0.  0.  0.  1.  1.]
  S0 = [ 0.  0.  0.  0.  0.  0.  1.  1.  1.  1.  1.  1.  2.  2.  2.  2.  2.  3.  3.  3.  3.  3.  3.]

Au niveau de la mise en oeuvre, vous définirez la méthode def initGD(X,N): qui prend un ensemble de séquences d'observations et qui retourne l'ensemble des séquences d'états. Pour chaque séquence x, vous pouvez utiliser:

np.floor(np.linspace(0,N-.00000001,len(x)))

Apprentissage

Etant donné la structure d'un MMC:

  • les observations n'influencent pas les états: les matrices {$\Pi, A$} s'obtiennent comme dans un modèle de Markov simple (cf semaine 6)
  • chaque observation ne dépend que de l'état courant

La nature des données nous pousse à considérer des lois de probabilités discrètes quelconques pour les émissions. L'idée est donc de procéder par comptage en définissant la matrice {$B$} comme suit:

  • K colonnes (nombre d'observations), N lignes (nombre d'états)
  • Chaque ligne correspond à une loi d'émission pour un état (ie, chaque ligne somme à 1)

Ce qui donne l'algorithme:

  1. {$b_{ij}$} = comptage des émissions depuis l'état {$s_i$} vers l'observation {$x_j$}
  2. normalisation des lignes de {$B$}

Donner le code de la fonction def learnHMM(allX, allS, N, K): qui apprend un modèle à partir d'un ensemble de couples (seq. d'observations, seq. d'états)

Variante stabilisée: afin d'éviter les transitions à probabilité nulle et de régulariser l'ensemble du système, il est intéressant d'initialiser les matrices de transition à une valeur non nulle. Cette initialisation spéciale peut être faite de manière optionnelle en utilisant les arguments par défaut de python:

def learnHMM(allx, allq, N, K, initTo0=False):
    if initTo0:
        A = np.zeros((N,N))
        B = np.zeros((N,K))
        Pi = np.zeros(N)
    else:
        eps = 1e-8
        A = np.ones((N,N))*eps
        B = np.ones((N,K))*eps
        Pi = np.ones(N)*eps
...

Validation sur les séquences de la classe A:

  K = 10 # discrétisation (=10 observations possibles)
  N = 5  # 5 états possibles (de 0 à 4 en python) 
  # Xd = angles observés discrétisés

  Pi, A, B = learnHMM(Xd[Y=='a'],q[Y=='a'],N,K)

  Pi=[ 1.  0.  0.  0.  0.]

  A=[[ 0.79  0.21  0.    0.    0.  ]
     [ 0.    0.76  0.24  0.    0.  ]
     [ 0.    0.    0.77  0.23  0.  ]
     [ 0.    0.    0.    0.76  0.24]
     [ 0.    0.    0.    0.    1.  ]]

  B=[[ 0.06  0.02  0.    0.    0.    0.    0.    0.04  0.49  0.4 ]
     [ 0.    0.04  0.    0.13  0.09  0.13  0.02  0.09  0.41  0.09]
     [ 0.    0.    0.    0.02  0.12  0.5   0.31  0.04  0.    0.  ]
     [ 0.07  0.    0.    0.    0.    0.    0.26  0.33  0.2   0.15]
     [ 0.73  0.12  0.    0.    0.    0.    0.    0.02  0.02  0.12]]

Viterbi (en log)

Rappels:

  • Viterbi sert à estimer la séquence d'états la plus probable étant donnés les observations et le modèle.
  • Viterbi peut servir à approximer la probabilité de la séquence d'observation étant donné le modèle.

1. Initialisation (avec les indices à 0 en python): {$$\begin{array}{ccccccccc} \delta_{0} (i) &=& \log \pi_{i} +\log b_{i} (x_{1}) \\ \Psi_{0}(i) &=& -1 \mbox{ Note: -1 car non utilisé normalement} \end{array}$$} 2. Récursion: {$$ \begin{array}{ccccccccc} \delta_{t} (j) &=&\displaystyle \left[\max_{i} \delta_{t-1}(i) + \log a_{ij}\right] + \log b_{j}(x_{t}) \\ \Psi_{t}(j) &=&\displaystyle \arg\max_{i\in [1,\ N]} \delta_{t-1} (i) + \log a_{ij} \end{array}$$} 3. Terminaison (indices à {$T-1$} en python) {$$ S^{\star} = \max_{i} \delta_{T-1}(i)$$} 4. Chemin {$$\begin{array}{ccccccccc} s_{T-1}^{\star} & = &\displaystyle \arg\max_{i} \delta_{T-1}(i) \\ s_{t}^{\star} & = & \displaystyle \Psi_{t+1}(s_{t+1}^{\star}) \end{array}$$}

L'estimation de {$\log p(x_0^{T-1} | \lambda)$} est obtenue en cherchant la plus grande probabilité dans la dernière colonne de {$\delta$}.
Donner le code de la méthode def viterbi(x,Pi,A,B):

validation en utilisant le modèle précédent {$\lambda_A$} (mêmes valeurs de N et K):

  s_est, p_est = viterbi(Xd[0], Pi, A, B)
  s_est = [ 0.  0.  0.  0.  0.  0.  0.  0.  1.  2.  2.  2.  2.  2.  3.  3.  3.  3.  4.  4.  4.  4.  4.]
  p_est = -38.0935655456

[OPT] Probabilité d'une séquence d'observation

En fonction de votre vitesse d'avancement, coder la version {$\alpha$} de l'estimation (exacte) des probabilités d'une séquence d'observations. Comparer avec l'approximation estimée précédemment (le résultat attendu est-il plus ou moins élevé avec la procédure exacte?).

validation

  p =  calc_log_pobs_v2(Xd[0],Pi, A, B)
  p = -34.8950422805

Apprentissage complet (Baum-Welch simplifié)

En utilisant la procédure de Baum-Welch simplifiée vue en TD et rappelée ci-dessous, proposer un code pour apprendre les modèles correspondant aux 26 lettres de l'alphabet.

Baum-Welch simplifié:

  1. Initialisation des états cachés arbitraire (eg méthode gauche-droite)
  2. Tant que critère de convergence non atteint
    1. Apprentissage des modèles {$\lambda_{lettre}=\{\Pi, A, B\}$}
    2. Estimation des états cachés par Viterbi

Le critère de convergence sera la vraisemblance.

  • A chaque itération {$k$} et pour toutes les lettres {$lettre$}, calculer pour l'ensemble des séquences d'observation :

{$$\log\mathcal L^k = \sum_{lettre}\sum_i \log p(\mathbf x_i^{lettre} | \lambda_{lettre}^k)$$}

  • Lorsque la vraisemblance n'évolue plus (ie {$\frac{\log\mathcal L^k - \log\mathcal L^{k+1}}{\log\mathcal L^k} < 1e-4$}), sortir de la boucle de mise à jour.
  • Donner l'implémentation de la méthode d'apprentissage
  • Tracer la courbe de l'évolution de la vraisemblance au cours des itérations

Evaluation des performances

Comme dans la séance précédente, il faut trouver un moyen d'évaluer les performances sans biaiser les résultats... Nous utilisons la même procédure.

  • Récupérer la méthode separeTrainTest(y, pc) documentée lors de la séance précédente.
  • Apprendre les modèles sur les échantillons d'apprentissage.
  • Evaluer les performances sur les données de test.

Le code de cette partie est grossièrement identique à la séance passée. Pour plus de détails : semaine 6 Comme lors de la semaine précédente, dessiner la matrice de confusion pour détecter les lettres faciles/difficiles à distinguer.

Partie optionnelle

Génération de lettres

Comme dans le TME précédent, proposer une procédure de génération de lettres.

  • Donner le code de generateHMM qui génère une séquence d'observations (et une séquence d'états) à partir des sommes cumulées des modèles de lettres (cf usage ci-dessous)
  • Faire tourner le code suivant qui réalise la génération de n échantillon pour nClred classes de lettres

# affichage d'une lettre (= vérification bon chargement)
def tracerLettre(let):
    a = -let*np.pi/180;
    coord = np.array([[0, 0]]);
    for i in range(len(a)):
        x = np.array([[1, 0]]);
        rot = np.array([[np.cos(a[i]), -np.sin(a[i])],[ np.sin(a[i]),np.cos(a[i])]])
        xr = x.dot(rot) # application de la rotation
        coord = np.vstack((coord,xr+coord[-1,:]))
    plt.plot(coord[:,0],coord[:,1])
    return

#Trois lettres générées pour 5 classes (A -> E)
n = 3          # nb d'échantillon par classe
nClred = 5   # nb de classes à considérer
fig = plt.figure()
for cl in xrange(nClred):
    Pic = models[cl][0].cumsum() # calcul des sommes cumulées pour gagner du temps
    Ac = models[cl][1].cumsum(1)
    Bc = models[cl][2].cumsum(1)
    long = np.floor(np.array([len(x) for x in Xd[itrain[cl]]]).mean()) # longueur de seq. à générer = moyenne des observations
    for im in range(n):
        s,x = generateHMM(Pic, Ac, Bc, int(long))
        intervalle = 360./d  # pour passer des états => angles
        newa_continu = np.array([i*intervalle for i in x]) # conv int => double
        sfig = plt.subplot(nClred,n,im+n*cl+1)
        sfig.axes.get_xaxis().set_visible(False)
        sfig.axes.get_yaxis().set_visible(False)
        tracerLettre(newa_continu)
plt.savefig("lettres_hmm.png")

Dilemme Biais-Variance, reflexion sur la complexité des modèles

Durant ces deux semaines, on a étudié des modèles de complexités croissantes: architecture plus complexe, nombre de paramètres en hausse. Cela nous amène à discuter le compromis biais-variance. On considère en général qu'un modèle pauvre:

  • est en général plus loin de la solution optimale {$h^\star$} (grand biais),
  • mais comme l'espace des hypothèses {$\mathcal H$} est peu étendu, les résultats sont relativement stable (faible variance).

A l'inverse, un modèle plus riche:

  • peut s'approcher plus près de la solution optimale {$h^\star$} (biais faible),
  • mais présente un problème de stabilité du fait de l'étendue de l'espace d'hypothèses {$\mathcal H$} à explorer... (grande variance)

Les solutions de régularisation (comme celle proposée dans l'algorithme de comptage) cherche à optimiser le compromis, c'est à dire à isoler une partie de l'espace des hypothèses {$\mathcal H$} particulièrement intéressante:

En faisant tourner plusieurs fois les algorithmes de chaines de Markov simples et plusieurs fois les algorithmes de MMC avec différentes valeurs de N et K:

  • calculer les moyennes et variances des performances obtenues,
  • analyser les résultats par rapport au compromis biais-variance