.. _chap-tpAfd:

#####################################################
Travaux pratiques - Analyse factorielle discriminante
#####################################################

(correspond à 1 séance de TP)

.. only:: html

    .. container:: notebook

        .. image:: cnam_theme/static/jupyter_logo.png
            :class: svg-inline

        `Cahier Jupyter <jupyter/tpAfd.ipynb>`_

Références externes utiles :

  * `Documentation NumPy <https://docs.scipy.org/doc/numpy/user/index.html>`_
  * `Documentation SciPy <https://docs.scipy.org/doc/scipy/reference/>`_
  * `Documentation MatPlotLib <http://matplotlib.org/>`_
  * `Site scikit-learn <http://scikit-learn.org/stable/index.html>`_
  * `Site langage python <https://www.python.org>`_

**L'objectif** de cette séance de TP est de présenter l'utilisation des fonctionnalités de Scikit-learn concernant l'AFD, ainsi que de contribuer à une meilleure compréhension de cette méthode d'analyse de données. Après une brève présentation de la classe correspondante de Scikit-learn, nous examinons d'abord l'application de l'AFD à des données générées de façon contrôlée, sans ou avec une réduction de dimension préalable par ACP. Deux exemples du cours, basés sur des données réelles, sont ensuite repris et analysés.


L'AFD dans Scikit-learn
***********************

Dans Scikit-learn, l'analyse factorielle discriminante (AFD) est mise en œuvre dans la classe `LinearDiscriminantAnalysis <http://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.LinearDiscriminantAnalysis.html>`_). Le `guide utilisateur <http://scikit-learn.org/stable/modules/lda_qda.html>`_ fournit quelques explications supplémentaires.

Pour définir l'analyse, on appelle ``LinearDiscriminantAnalysis(solver='svd', shrinkage=None, priors=None, n_components=None, store_covariance=False, tol=0.0001)``.

.. admonition:: Note

    Attention, la classe ``sklearn.lda.LDA`` est obsolète (*deprecated*). Il ne faut plus l'utiliser.
    
Les paramètres d'appel sont (pour des explications plus détaillées voir le lien ci-dessus) :

 * ``solver`` : solveur employé pour obtenir les axes ; ``'svd'`` traite directement la matrice de données, sans calculer au préalable la matrice des covariances empiriques ; ``'lsrq'`` est la solution des moindres carrés ; ``'eigen'`` calcule d'abord la matrice des covariances empiriques (par défaut, ``solver = 'svd'``).
 * ``shrinkage`` : si absent, la matrice des covariances empiriques est directement utilisée ; si ``'auto'``, choix automatique du paramètre suivant le lemme Ledoit-Wolf [LWH04]_.
 * ``priors`` : optionnel, array ``n_classes`` avec les *a priori* des différentes classes.
 * ``n_components`` : nombre de composantes à conserver (par défaut toutes, c'est à dire ``< n_classes - 1)``).
 * ``store_covariance`` : optionnel, si ``True`` demande le calcul de la maytrice des covariances empiriques.
 * ``tol`` : optionnel, seuil à utiliser pour l'estimation du rang dans le solveur SVD.

Les attributs accessibles sont les suivants :

 * ``coef_`` : array, unidimensionnel ``n_features`` ou bidimensionnel ``n_classes x n_features``, poids des variables, éventuellement par classe.
 * ``intercept_`` : array ``n_features``, intercept.
 * ``covariance_`` : array ``n_features x n_features``, matrice des covariances empiriques (partagée par toutes les classes).
 * ``explained_variance_ratio_`` : array ``n_components``, pourcentage de variance expliquée par chaque composante (ou facteur discriminant).
 * ``means_`` : array ``n_classes x n_features``, centres de gravité des classes.
 * ``priors_`` : array ``n_classes`` avec les *a priori* des différentes classes (leur somme est égale à 1).
 * ``scalings_`` : array ``rank x n_classes - 1``, transformation des variables dans l'espace des centres de gravité des classes.
 * ``xbar_`` : array ``n_features``, centre de gravité global.
 * ``classes_`` : array ``n_classes``, étiquettes des classes.

Les méthodes qui peuvent être employées :

* ``decision_function(X)`` : scores de confiance pour les données.
* ``fit(X, y[, store_covariance, tol])`` : construction du modèle à partir des données X et des étiquettes de classe y.
* ``fit_transform(X[, y])`` : construction du modèle et ensuite projection des données sur l'espace discriminant.
* ``get_params([deep])`` : lire les valeurs des paramètres de l'estimateur employé.
* ``predict(X)`` : prédire les étiquettes de classe pour les données de la matrice X (par ex. données de test).
* ``predict_log_proba(X)`` : estimer les log-probabilités.
* ``predict_proba(X)`` : estimer les probabilités.
* ``score(X, y[, sample_weight])`` : taux de bon classement des données X, avec les étiquettes désirées dans y.
* ``set_params(**params)`` : donner des valeurs aux paramètres de l'estimateur employé.
* ``transform(X)`` : projeter les données sur l'espace discriminant.


AFD de données générées
***********************

Démarrez par la génération de deux ensembles de 500 vecteurs chacun dans l'espace 3D, suivant une loi normale (de moyenne nulle et de variance unitaire), ensuite transformez les points générés :

.. code-block:: python
    
    import numpy as np    # si pas encore fait
    s1 = np.array([[3,0,0],[0,1,0],[0,0,0.01]])  # matrice de déformation
    r1 = np.array([[0.36,0.48,-0.8],[-0.8,0.6,0],[0.48,0.64,0.6]])  # matrice de rotation

    # On génère deux nuages de points déformés et tournés
    rndn3d1 = np.random.randn(500,3)
    rndef1 = rndn3d1.dot(s1).dot(r1)
    rndn3d2 = np.random.randn(500,3)
    # Le deuxième nuage de point est translaté selon l'axe z
    rndef2 = rndn3d2.dot(s1).dot(r1) + [0, 0, 1]
    rndef = np.concatenate((rndef1, rndef2))
    print(rndef.shape)

.. admonition:: Question :

   Expliquez l'addition qui permet de créer ``rndef2`` : à quoi sert cette addition, quel est le résultat de l'addition entre deux ``array`` de dimension différente ?

.. ifconfig:: tpml in ('public')

   .. admonition:: Correction :

           Cette addition sert à décaler un des deux nuages (représenté par les données de ``rndef2``) sur l'axe ``z`` et ainsi à séparer les données des deux classes. Pour l'addition, le vecteur ``[0, 0, 1]`` est (implicitement) utilisé pour produire une matrice de 1000 lignes (comme ``rndn3d2.dot(s1).dot(r1)``) dont chaque ligne est égale à ce vecteur.

Maintenant, générez des étiquettes de classe correspondantes : *1* pour les données de ``rndef1``, *2* pour celles de ``rndef2`` et ensuite visualisez les données :

.. code-block:: python
    
    lcls1 = np.ones(500)
    lcls2 = 2 * np.ones(500)
    lcls = np.concatenate((lcls1, lcls2))
    print(lcls)

Cela indique que les 500 premiers points sont de classe 1 et les 500 points suivants sont de classe 2.

Affichons maintenant le nuage de points avec matplotlib :
    
.. code-block:: python

    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(rndef[:,0], rndef[:,1], rndef[:,2], c=lcls)
    plt.show()

Appliquez d'abord une analyse discriminante à ces données, projetez-les sur l'espace discriminant (de quelle dimension est-il ?) et ensuite visualisez les résultats :

.. code-block:: python
    
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

    lda = LinearDiscriminantAnalysis()
    lda.fit(rndef,lcls)

L'AFD a ainsi été entraînée sur les données. Nous pouvons maintenant obtenir les données projetées grâce à la méthode ``transform`` :

.. code-block:: python

    rndt = lda.transform(rndef)
    print(rndt.shape)

Ce qui produit un vecteur de dimension (1000, 1), c'est-à-dire que l'AFD a projeté nos observations sur un unique axe discriminant.

Affichons désormais le résultat de cette projection :

.. code-block:: python

    plt.plot(rndt, lcls, 'r+')
    plt.xlabel("Axe discriminant")
    plt.ylabel("Classe")
    plt.yticks([1, 2])
    plt.show()

Essayez de réduire la dimension des données avec une ACP en conservant les projections sur les 2 premiers axes principaux, appliquez ensuite l'AFD et visualisez les résultats :

.. code-block:: python
    
    from sklearn.decomposition import PCA
    pca = PCA(n_components=2)
    pca.fit(rndef)

Projection des observations puis affichage du graphe :

.. code-block:: python

    rndp = pca.transform(rndef)
    print(rndp.shape)
    # (1000, 2)
    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(rndp[:,0], rndp[:,1], c=lcls)
    plt.show()

On applique ensuite une AFD par-dessus le résultat de cette projection :
    
.. code-block:: python

    lda = LinearDiscriminantAnalysis()
    lda.fit(rndp,lcls)
    rndpt = lda.transform(rndp)

    plt.plot(rndpt, lcls, 'r+')
    plt.xlabel("Axe discriminant")
    plt.ylabel("Classe")
    plt.yticks([1, 2])
    plt.show()

.. admonition:: Question :

   Expliquez ce que vous constatez.

.. ifconfig:: tpml in ('public')

   .. admonition:: Correction :

           L'ACP préalable a supprimé la direction qui permettait de séparer les deux classes car dans cette direction la variance du nuage global était la plus faible. En conséquence, les classes ne sont plus séparables pour les projections des observations sur ce premier plan factoriel.


ACP et AFD des données « textures »
***********************************

Ces données correspondent à 5500 observations qui sont des pixels extraits d'images de textures et décrits par 40 variables (descripteurs visuels). Chaque observation appartient à une des 11 classes de textures ; chaque classe est représentée par 500 pixels. Les données sont issues de `https://www.elen.ucl.ac.be/neural-nets/Research/Projects/ELENA/databases/REAL/texture/ <https://www.elen.ucl.ac.be/neural-nets/Research/Projects/ELENA/databases/REAL/texture/>`_

.. only:: jupyter

    .. code-block:: bash

        !wget -nc http://cedric.cnam.fr/~crucianm/src/texture.dat

.. only:: not jupyter

    .. code-block:: bash

        wget -nc http://cedric.cnam.fr/~crucianm/src/texture.dat

Examinez le fichier téléchargé et regardez sur le site indiqué ci-dessus les explications concernant le format des données. Le jeu de données comporte 5500 observations comportant 41 colonnes : 40 attributs visuels décrivant le type de texture et 1 attribut concernant à la catégorie de la texture.

.. admonition:: Question :

   Lisez ces données dans un ``array`` en utilisant ``np.loadtxt``. Appliquez une ACP et visualisez les projections sur les 3 premiers axes principaux. Appliquez également une AFD (sur les données initiales, non sur les projections issues de l'ACP) et visualisez les projections sur les 3 premiers axes discriminants.

.. ifconfig:: tpml in ('public')

   .. admonition:: Correction :

        .. code-block:: python
        
            textures = np.loadtxt('texture.dat')

            pca = PCA()
            pca.fit(textures[:,:40])
            texturesp = pca.transform(textures[:,:40])

            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter(texturesp[:,0], texturesp[:,1], texturesp[:,2], c=textures[:,40])
            plt.title("Projection après ACP")
            plt.show()

            lda.fit(textures[:,:40],textures[:,40])
            texturest = lda.transform(textures[:,:40])

            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter(texturest[:,0], texturest[:,1], texturest[:,2], c=textures[:,40])
            plt.title("Projection après AFD")
            plt.show()


.. admonition:: Question :

    Faites l'extraction d'une partie des observations avant d'appliquer l'AFD. Réalisez ensuite la prédiction (étape décisionnelle de l'AFD) pour ces observations et affichez l'erreur de prédiction. Pour l'extraction des observations vous pouvez employer (entre autres) :

    .. code-block:: python
    
        from sklearn.model_selection import train_test_split
        train, test = train_test_split(textures, test_size = 0.2)

.. ifconfig:: tpml in ('private')

    .. only:: jupyter

        .. code-block:: python


.. ifconfig:: tpml in ('public')

    .. admonition:: Correction :

        .. code-block:: python
            
            from sklearn.model_selection import train_test_split
            train, test = train_test_split(textures, test_size = 0.2)
            print(train.shape)
            # (4400, 41)
            print(test.shape)
            # (1100, 41)

            lda.fit(train[:,:40],train[:,40])
            lda.score(test[:,:40],test[:,40])
            # 0.99818181818181817


AFD des données « feuilles »
********************************

Nous appliquerons enfin l'AFD sur les données concernant des feuilles d'arbres, issues de `https://archive.ics.uci.edu/ml/datasets/Leaf <https://archive.ics.uci.edu/ml/datasets/Leaf>`_. Si le données ne se trouvent pas déjà dans le répertoire de travail, entrez dans une fenêtre terminal :

.. only:: jupyter

    .. code-block:: bash

        !wget -nc http://cedric.cnam.fr/~crucianm/src/leaf.csv

.. only:: not jupyter

    .. code-block:: bash

        wget -nc http://cedric.cnam.fr/~crucianm/src/leaf.csv

Examinez une nouvelle fois le fichier téléchargé et regardez sur le site indiqué ci-dessus les explications concernant le format des données.

.. code-block:: python
    
    leaf = np.loadtxt('leaf.csv', delimiter=',')
    print(leaf[:2,:])

    lda = LinearDiscriminantAnalysis()
    lda.fit(leaf[:,2:],leaf[:,0])
    
    leafl = lda.transform(leaf[:,2:])
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(leafl[:,0],leafl[:,1],leafl[:,2],c=leaf[:,0])
    plt.show()


.. admonition:: Question :

   Comparez les résultats avec ceux de l'ACP.

.. ifconfig:: tpml in ('public')

   .. admonition:: Correction :

          Pour comparer il est nécessaire de reprendre l'affichage des projections des observations sur les trois axes principaux en utilisant pour chaque point la couleur correspondant à la classe, comme nous l'avons fait ci-dessus pour les projections obtenues par AFD. En comparant les deux graphiques on constate que l'AFD sépare un peu mieux les classes, mais celles-ci restent difficiles à déceler car leur nombre est très important (40 !) et donc les nuances de couleur correspondant à des classes différentes sont parfois très proches. L'utilisation de la couleur pour distinguer les classes a atteint ses limites.

.. admonition:: Question :

   Faites l'extraction d'une partie des observations avant d'appliquer l'AFD. Réalisez ensuite la prédiction (étape décisionnelle de l'AFD) pour ces observations et affichez l'erreur de prédiction.

.. ifconfig:: tpml in ('private')

    .. only:: jupyter

        .. code-block:: python


.. ifconfig:: tpml in ('public')

   .. admonition:: Correction :

           .. code-block:: python
            
            from sklearn.model_selection import train_test_split
            train, test = train_test_split(leaf, test_size = 0.2)
            lda.fit(train[:,2:],train[:,0])
            lda.score(test[:,2:],test[:,0])
            # 0.8970588235294118

|

--------------------------------

.. [LWH04] Ledoit O., Wolf M. *Honey, I Shrunk the Sample Covariance Matrix*. *Journal of Portfolio Management* 30(4), 110-119, 2004.