Cours

TME 8: Méthodes discriminantes

Nous allons travailler dans cette séance sur un comparatif des modèles génératifs et discriminants pour la classification de chiffres manuscrits. Les résultats de référence sont ceux obtenus avec un modèle génératif issu de la séance 3 (cf sujet semaine 3).

Résultats de référence avec la méthode génératrice

Afin que tout le monde reparte de la même base, nous redonnons les données (au format pickle) et un code pour la résolution:

Données USPS (zipées): usps.pkl.zip

import numpy as np
import pickle as pkl

# fonction de suppression des 0 (certaines variances sont nulles car les pixels valent tous la même chose)
def woZeros(x):
    y = np.where(x==0., 1., x)
    return y

# Apprentissage d'un modèle naïf où chaque pixel est modélisé par une gaussienne (+hyp. d'indépendance des pixels)
# cette fonction donne 10 modèles (correspondant aux 10 classes de chiffres)
# USAGE: theta = learnGauss ( X,Y )
# theta[0] : modèle du premier chiffre,  theta[0][0] : vecteur des moyennes des pixels, theta[0][1] : vecteur des variances des pixels
def learnGauss (X,Y):
    theta = [(X[Y==y].mean(0),woZeros(X[Y==y].var(0))) for y in np.unique(Y)]
    return (np.array(theta))

# Application de TOUS les modèles sur TOUTES les images: résultat = matrice (nbClasses x nbImages)
def logpobs(X, theta):
    logp = [[-0.5*np.log(mod[1,:] * (2 * np.pi )).sum() + -0.5 * ( ( x - mod[0,:] )**2 / mod[1,:] ).sum () for x in X] for mod in theta ]
    return np.array(logp)

######################################################################
#########################     script      ############################


# Données au format pickle: le fichier contient X, XT, Y et YT
# X et Y sont les données d'apprentissage; les matrices en T sont les données de test
data = pkl.load(file('ressources/usps_small.pkl','rb'))

X = data['X']
Y = data['Y']
XT = data['XT']
YT = data['YT']

theta = learnGauss ( X,Y ) # apprentissage

logp  = logpobs(X, theta)  # application des modèles sur les données d'apprentissage
logpT = logpobs(XT, theta) # application des modèles sur les données de test

ypred  = logp.argmax(0)    # indice de la plus grande proba (colonne par colonne) = prédiction
ypredT = logpT.argmax(0)

print "Taux bonne classification en apprentissage : ",np.where(ypred != Y, 0.,1.).mean()
print "Taux bonne classification en test : ",np.where(ypredT != YT, 0.,1.).mean()

Validation (du copier-coller!)

  Taux bonne classification en apprentissage :  0.826498422713
  Taux bonne classification en test :  0.782760338814

Approche discriminante

Nous allons coder l'algorithme de la régression logistique, vu en TD. Cet algorithme itératif de classification binaire présente les caractéristiques suivantes:

{$$ f (\mathbf{x_i}) = p(y=1 | \mathbf{x_i} ) = \frac{1}{1 + \exp( -( \mathbf{x_i} \mathbf{w} + b))} $$}

Formule de la log-vraisemblance:

{$$ \mathcal L_{\log}= \sum_i y_i log (p(y=1 | \mathbf{x_i} )) + (1-y_i) ( \log(1- p(y=1 | \mathbf{x_i} ))) $$} {$$ \mathcal L_{\log}= \sum_i y_i log ( f (\mathbf{x_i}) ) + (1-y_i) (\log(1- f (\mathbf{x_i}))) $$}

Problème d'apprentissage:

{$$ (w^\star, b^\star) = \arg\max_{w,b} \mathcal L_{\log} $$}

Malheureusement, pas de solution analytique (cf TD)... Nous utiliserons donc une résolution itérative basée sur la dérivation de la vraisemblance (montée de gradient, détaillée ci-dessous).

Formule de dérivation de la vraisemblance:

{$$ \frac{\partial }{\partial w_j} \mathcal L_{\log} =\sum_{i=1}^N x_{ij}( y_i- f (\mathbf{x_i}))$$}

{$$ \frac{\partial }{\partial b} \mathcal L_{\log} = \sum_{i=1}^N (y_i - f (\mathbf{x_i}))$$}

Apprentissage d'un modèle binaire

L'apprentissage de la régression logistique est non-trivial:

  • les formules de gradient ne sont pas annulables analytiquement
  • nous allons chercher à augmenter progressivement la vraisemblance du modèle (cf figure ci-dessous)
    • Initialisation aléatoire des paramètres {$\lambda_0 = \{w_0, b_0\}$}
    • Tant que la convergence n'est pas atteinte:
      • Mise à jour:

{$$ \forall j,\ w_j = w_j + \epsilon \frac{\partial }{\partial w_j} \mathcal L_{\log}$$} {$$ b = b + \epsilon \frac{\partial }{\partial b} \mathcal L_{\log}$$}

Précisions diverses:

  • il faut donner {$\epsilon$} à l'algorithme: tester des valeurs entre {$10^{-3}$} et {$10^{-5}$}
  • pour l'initialisation de {$\lambda_0 = \{w_0, b_0\}$}, il faut choisir des valeurs petites (de l'ordre de {$\epsilon$} par exemple)
  • il faut que les valeurs de {$\lambda_0 $} soient positives et négatives (elles n'ont pas de raison d'être uniquement positives)
  • pour faciliter le codage initial, coder en dur le nombre d'itérations à réaliser (100 par exemple) au lieu de vérifier la convergence
  • pour vérifier que votre algorithme fonctionne bien:
    • calculer la vraisemblance à chaque itération et la stocker
    • tracer l'évolution de la vraisemblance à la fin des itérations

Résultats envisageables en jouant sur {$\epsilon$} et le nombre d'itération (un seul est satisfaisant!):


Divergence ({$\epsilon$} trop grand)

Convergence incomplète: {$\epsilon$} petit ou nbIter faible

bonne convergence

Indication pour refaire les courbes

Les valeurs qui marchent bien sont: {$\epsilon=5e-5$} et nombre d'itérations = 120. Pour afficher ce qui se passe en cas de divergence, il faut éviter les valeurs trop faibles qui font des nan qui ne sont pas affichables... J'ai un peu triché dans le calcul de vraisemblance en introduisant le calcul robuste suivant:

(Y * np.log(np.maximum(f,1e-8)) + (1-Y)*np.log(1-np.minimum(f,1-1e-8))).sum()

Paradigme un-contre-tous pour le passage au multi-classe

Nous avons mis en place un classifieur binaire alors que le problème qui nous intéresse compte 10 classes. L'idée est:

  • d'apprendre 10 classifieurs : 1 classe contre tout le reste des données
# Pour transformer le vecteur Y afin de traiter la classe 0:
Yc0 = np.where(Y==0,1.,0.)
  • chaque classifieur {$f_c$} donne un score (la probabilité d'appartenir à la classe cible c)

{$$f_c(\mathbf{x_i}) = p(y=1 | \mathbf{x_i}) $$}

  • en stockant tous les scores dans une matrice, nous nous retrouvons dans le même cas de figure qu'avec les modèles génératifs
# si vos paramètres w et b, correspondant à chaque classe, sont stockés sur les lignes de thetaRL... Alors:

pRL  = np.array([[1./(1+np.exp(-x.dot(mod[0]) - mod[1])) for x in X] for mod in thetaRL ])
pRLT = np.array([[1./(1+np.exp(-x.dot(mod[0]) - mod[1])) for x in XT] for mod in thetaRL ])
ypred  = pRL.argmax(0)
ypredT = pRLT.argmax(0)
print "Taux bonne classification en apprentissage : ",np.where(ypred != Y, 0.,1.).mean()
print "Taux bonne classification en test : ",np.where(ypredT != YT, 0.,1.).mean()

Calculer la performance en apprentissage et en test de la regression logistique sur les données USPS et comparer avec les modèles génératifs.

Tracer des matrices de confusion à la manière de ce qui a été fait dans les dernières séances.

Note: le code est quasi-identique.

Analyse qualitative comparée des modèles (génératifs/discriminants)

Pour le modèle génératif codé en début de séance:

Ce modèle est bien un modèle génératif: en tirant des pixels aléatoirement selon la loi d'un chiffre {$\mathcal N(\mu_{chiffre}, \sigma^2_{chiffre})$}, il est possible de générer une infinité de chiffres, tous différents les uns des autres mais tous issus du même modèle!

Comparer plusieurs chiffres générer à partir du même modèle et afficher le vecteur moyen du modèle.

import matplotlib.pyplot as plt
import matplotlib.cm as cm

# pour afficher le vecteur des moyennes des pixels pour le modèle 0:
plt.figure()
plt.imshow(theta[0][0].reshape(16,16), cmap = cm.Greys_r, interpolation='nearest')
plt.show()

Note: pour tirer un nombre x selon une loi gaussienne de paramètre {$\mathcal N(\mu_{chiffre}, \sigma^2_{chiffre})$}:

x=np.random.randn(1) * np.sqrt(sigma2) + mu

Note 2: dans un premier temps, nous avions détecté les variances nulles pour les mettre à 1 afin de pouvoir calculer les {$\log p(x|\lambda)$} sans tomber toujours sur nan... Mais cette aménagement doit être supprimé lors de la génération de chiffre.

Note 3: les pixels générés peuvent être positifs... Ou négatifs. Alors que dans la base de données, tous les pixels sont positifs... Du coup, pour obtenir un résultat plus vraisemblable, il est intéressant de ne garder que les valeurs positives (et de mettre à 0 les valeurs négatives).

Pour le modèle discriminant:

Visualiser le vecteur {$w$} de la régression logistique pour chaque classe: cela vous indique les pixels qui ''comptent' dans la prise de décision. En étudiant {$f_2$}, par défaut, vous visualisez ce qui distingue un 2 de tous les autres chiffres. Il est possible de pousser plus loin l'analyse qualitative en apprenant une régression logistique spécifique aux 2 et 3 par exemple pour comprendre comment le modèle distingue ces deux classes.

Zones blanches = si un pixel est présent, le modèle augmente la probabilité de reconnaitre un 2.

Optimisation du modèle discriminant

Mettre en place un critère d'arrêt basé sur l'évolution de la vraisemblance dans la procédure évolutive. L'idée est la même que dans les séances passée (cf TME de la semaine 7).

Mettre en place une stratégie de rejet des échantillons ambigus: étudier l'amélioration des résultats en fonction du nombre d'échantillon rejetés. Etudier en particulier les cas suivants:

  • Parmi tous les {$p(y_c=1 | \mathbf{x_i})$}, aucun ne dépasse 0.5 (ie, pour le modèle, le chiffre présenté ne ressemble à rien)
  • Parmi tous les {$p(y_c=1 | \mathbf{x_i})$}, le max et le second plus grand sont très proches (seuil à définir)

En général, on présente les résultats sous la forme d'une courbe de pourcentage de reconnaissance (en ordonnées) qui évolue en fonction du nombre d'échantillon rejetés.