.. _chap-tpAcp:

######################################################
Travaux pratiques - Analyse en composantes principales
######################################################
 
(correspond à 1 séance de TP)

.. only:: html

    .. container:: notebook

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

        `Cahier Jupyter <jupyter/tpAcp.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'ACP, 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 des données générées de façon contrôlée afin de mieux comprendre l'impact de la distribution des données sur les résultats de l'ACP. Deux exemples du cours, basés sur des données réelles, sont ensuite repris et examinés.


L'ACP dans Scikit-learn
***********************

Dans Scikit-learn, l'analyse en composantes principales (ACP) est mise en œuvre dans la classe ``sklearn.decomposition.PCA`` (voir `http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html <http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html>`_). Les variantes et quelques exemples sont présentés sur `http://scikit-learn.org/stable/modules/decomposition.html <http://scikit-learn.org/stable/modules/decomposition.html>`_.

Pour définir l'analyse on appelle ``PCA(n_components=None, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', random_state=None)``. Les paramètres d'appel sont (pour des explications plus détaillées voir le lien ci-dessus) :

 * ``n_components`` : nombre de composantes à conserver (par défaut toutes, c'est à dire ``min(n_samples, n_features)`` où ``n_samples`` est le nombre de lignes de la matrice de données et ``n_features`` le nombre de colonnes).
 * ``copy`` : si ``False``, la matrice de données est écrasée par les données transformées (par défaut ``True``).
 * ``whiten`` : si ``True``, les données transformées sont modifiées pour que les projections sur les axes principaux présentent une variance unitaire (par défaut ``False``).
 * ``svd_solver`` : solveur à employer parmi ``'auto'``, ``'full'``, ``'arpack'``, ``'randomized'`` (par défaut ``'auto'``).
 * ``tol``, ``iterated_power``, ``random_state`` : paramètres optionnels pour solveurs spécifiques.

Les attributs accessibles sont les suivants :

 * ``explained_variance_`` : array ``[n_components]`` dans lequel on trouve les valeurs propres en ordre décroissant.
 * ``components_`` : array ``[n_components, n_features]`` dans lequel chaque ligne correspond à un vecteur propre ; l'ordre est celui des valeurs propres.
 * ``explained_variance_ratio_`` : array ``[n_components]`` avec les valeurs propres exprimées en pourcentage de la variance expliquée.
 * ``mean_`` : array ``[n_features]`` contenant les moyennes des variables calculées lors de l'ACP.
 * ``n_components_`` : nombre estimé de composantes.
 * ``noise_variance_`` : variance du bruit estimée suivant l'approche de l'ACP probabiliste (voir lien ci-dessus).

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

 * ``fit(X[, y])`` : trouver les axes principaux et les valeurs propres associés à partir de la matrice de données X.
 * ``fit_transform(X[, y])`` : trouver les axes principaux et les valeurs propres associés, appliquer la projection sur les axes principaux (avec stockage des résultats dans la matrice de données X si ``copy=False`` et dans une autre matrice si ``copy=True``).
 * ``get_covariance()`` : calculer les covariances à partir de la matrice de données.
 * ``get_params([deep])`` : retourner les paramètres de l'estimateur.
 * ``get_precision()`` : calculer la matrice de précision des données.
 * ``inverse_transform(X[, y])`` : re-projeter les données transformées dans l'espace initial.
 * ``score(X[, y])`` : valeur moyenne de la log-vraisemblance pour toutes les observations.
 * ``score_samples(X)`` : log-vraisemblance pour chaque observation.
 * ``set_params(**params)`` : donner des valeurs aux paramètres de l'estimateur.
 * ``transform(X[, y])`` : appliquer la projection sur les axes principaux (avec stockage des résultats dans la matrice de données X si ``copy=False`` et dans une autre matrice si ``copy=True``).


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

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

.. code-block:: python
    
    import numpy as np
    # Génération de données selon une loi normale tridimensionnelle
    rndn3d = np.random.randn(500,3)
    
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    # Affichage du nuage de points
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(rndn3d[:,0], rndn3d[:,1], rndn3d[:,2])
    plt.title("Données initiales")
    plt.show()

Appliquez ensuite l'ACP à ces données :

.. code-block:: python
    
    from sklearn.decomposition import PCA

    pca = PCA(n_components=3)
    pca.fit(rndn3d)

    print("Pourcentage de variance expliquée : ")
    print(pca.explained_variance_ratio_)
    print("Composantes principales : ")
    print(pca.components_)

Chaque ligne de ``pca.components_`` représente le vecteur unitaire qui donne la direction d'un axe factoriel. Le vecteur ``pca.components_[i,:]`` correspond à la valeur propre ``pca.explained_variance_ratio_[i]``.

.. admonition:: Question :

   Générez un **autre** tableau de données en utilisant de nouveau ``np.random.randn(500,3)`` et appliquez l'ACP à ces nouvelles données. Comparez les résultats avec ceux obtenus sur ``rndn3d``. Que constatez-vous ? Expliquez.

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

   .. admonition:: Correction :

        On constate d'abord que les trois valeurs propres sont chaque fois très proches entre elles, ce qui indique que le nuage d'observations a une forme relativement sphérique. Ensuite, en examinant les vecteurs propres correspondants, on constate que d'un ensemble de données générées à un autre les directions des vecteurs propres changent, donc ces vecteurs n'ont pas de pertinence particulière pour décrire les données. Ce qui est cohérent avec les faibles différences entre valeurs propres successives.

Appliquez maintenant une déformation et une rotation dans l'espace tridimensionnel au nuage des observations de ``rndn3d``, visualisez le résultat :

.. code-block:: python
    
    s1 = np.array([[3,0,0],[0,1,0],[0,0,0.2]])  # 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
    rndef = rndn3d.dot(s1).dot(r1)
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(rndef[:,0], rndef[:,1], rndef[:,2])
    plt.title("Données déformées")
    plt.show()

Appliquez une ACP à ces observations :

.. code-block:: python
    
    pca = PCA(n_components=3)
    pca.fit(rndef)
    
    print("Pourcentage de variance expliquée : ")
    print(pca.explained_variance_ratio_)
    print("Composantes principales : ")
    print(pca.components_)

.. admonition:: Question :

   Que constatez-vous par rapport à l'application directe de l'ACP sur les données de ``rndn3d`` ? Expliquez. Peut-on lier les rapports entre les valeurs propres obtenues aux rapports entre les éléments de la diagonale de la matrice de déformation ``s1`` ?

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

   .. admonition:: Correction :

        La décroissance des valeurs propres est ici très rapide, le premier axe factoriel explique 90% de la variance. Les rapports entre valeurs propres correspondent aux carrés des rapports entre les éléments de la diagonale de la matrice de déformation. En effet, les valeurs propres sont des variances alors que les éléments de la diagonale de la matrice de déformation ont ici une signification d'écart-types.

.. admonition:: Question :

   Générez un autre ensemble de données en utilisant les mêmes matrices de transformation et appliquez l'ACP à ces nouvelles données. Comparez les résultats avec ceux obtenus sur ``rndef``. Que constatez-vous ? Expliquez.

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

   .. admonition:: Correction :

        Les valeurs propres obtenues sur le second ensemble de données sont très proches de celles obtenues sur ``rndef``. De plus, entre les deux ensembles de données, les vecteurs propres correspondant aux valeurs propres de même ordre sont très similaires. Les valeurs et vecteurs propres caractérisent les (fortes) déformations appliquées aux nuages d'observations générés et sont donc relativement stables par rapport à un changement d'échantillon.


ACP des données sur le sommeil des mammifères
*********************************************

L'ACP sera d'abord appliquée aux données concernant le sommeil des mammifères, issues de [AC76]_ et déjà vues dans le cours. Pour cela, il est d'abord nécessaire d'enregistrer ces données dans le répertoire de travail. Dans une fenêtre terminal entrez (ou télécharger le fichier via `ce lien <http://cedric.cnam.fr/~crucianm/src/mammals.csv>`_:

.. only:: html

    .. code-block:: bash
    
        wget -nc http://cedric.cnam.fr/~crucianm/src/mammals.csv

.. only:: jupyter

    .. code-block:: bash
    
        !wget -nc http://cedric.cnam.fr/~crucianm/src/mammals.csv

Vous pouvez examiner le contenu de ce fichier format ``.csv``. Les données doivent être ensuite lues dans Python :

.. code-block:: python
    
    mammals = np.loadtxt('mammals.csv', delimiter=';', usecols=[1,2,3,4,5,6,7,8,9,10], skiprows=1)
    print(mammals[:2,:])
    noms = np.genfromtxt('mammals.csv', dtype='str', delimiter=';', usecols=[0], skip_header=1)
    print(noms)

.. admonition:: Question :

   Appliquez une ACP aux données de ``mammals`` et affichez les valeurs propres. Comparez les deux premières aux valeurs vues en cours. Que constatez-vous ? Pourquoi ?

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

    .. only:: jupyter

        .. code-block:: python

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

    .. admonition:: Correction :

        .. code-block:: python
            
            pca = PCA()
            pca.fit(mammals)

            print(pca.explained_variance_ratio_)
            # array([  9.61492649e-01,   3.35026815e-02,   4.90760735e-03,
            #          8.11547630e-05,   1.21589728e-05,   2.12907093e-06,
            #          8.08341160e-07,   4.40101019e-07,   3.18913491e-07,
            #          5.23943276e-08])
          
        La première valeur propre obtenue explique la quasi-totalité de la variance, alors que dans le cours elle expliquait moins de la moitié. L'ACP directe n'est pas une ACP normée, les variances des variables « poids » (du corps, du cerveau) sont comparativement très élevées par rapport aux variances des autres variables et dominent très largement dans la définition de la direction du premier axe factoriel. Pour avoir une analyse utile il faut centrer et réduire les variables.

.. admonition:: Question :

   Appliquez une transformation de centrage et réduction aux variables décrivant les données de ``mammals`` (regardez `le centrage et la réduction des variables (« standardization ») avec scikit-learn <http://scikit-learn.org/stable/modules/preprocessing.html#standardization-or-mean-removal-and-variance-scaling>`_) et conservez le résultat dans ``mammalsCR``. Appliquez de nouveau l'ACP, de préférence en utilisant une nouvelle instance ``pcaCR`` pour conserver ``pca``. Comment ont évolué les valeurs propres les plus grandes ? Affichez le graphique de décroissance des valeurs propres.

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

    .. only:: jupyter

        .. code-block:: python

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

   .. admonition:: Correction :

        .. code-block:: python
            
            from sklearn import preprocessing

            mammalsCR = preprocessing.scale(mammals)
            pcaCR = PCA()
            pcaCR.fit(mammalsCR)

            ratios = pcaCR.explained_variance_ratio_
            print(ratios)
            # array([ 0.48017524,  0.22078825,  0.12605711,  0.06488067,  0.04739057,
            #         0.02607844,  0.01797744,  0.00920976,  0.00501433,  0.00242818])
            plt.bar(range(len(ratios)), ratios)
            plt.xticks(range(len(ratios)))
            plt.xlabel("Composante principale")
            plt.ylabel("% de variance expliquée")
            plt.show()

        La plus grande valeur propre a maintenant une valeur qui correspond au pourcentage vu en cours.

Nous pouvons maintenant obtenir et afficher les projections des observations sur les deux premiers axes factoriels. Rappelons que ``mammalsCR`` contient les données transformées (variables centrées et réduites) et que ``pcaCR`` est l'instance de PCA utilisée pour analyser ``mammalsCR`` ; ``mammalsCR`` et ``pcaCR`` ont été obtenues en réponse à la question ci-dessus.

.. code-block:: python
    
    mNt = pcaCR.transform(mammalsCR)
    print(mNt[:2,:])
    
.. code-block:: python

    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111)
    for i in range(len(noms)):
            x,y = mNt[i,0], mNt[i,1]
            ax.scatter(x,y)
            ax.text(x,y,noms[i])
       
    plt.show()

.. admonition:: Question :

   Réalisez l'affichage des projections sur les **3** premiers axes principaux ; la possibilité de rotation du graphique facilite la lecture des étiquettes textuelles.

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

    .. only:: jupyter

        .. code-block:: python


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

    .. admonition:: Correction :

        .. only:: jupyter

            .. code-block:: python

                %matplotlib notebook
                # Active le mode interactif de Matplotlib pour les notebooks

        .. code-block:: python
            
            fig = plt.figure(figsize=(8,8))
            ax = fig.add_subplot(111, projection='3d')
            for i in range(len(noms)):
                    x,y,z = mNt[i,0],mNt[i,1],mNt[i,2]
                    ax.scatter(x,y,z)
                    ax.text(x,y,z,noms[i])
                
            plt.show()

.. admonition:: Question :

   Les résultats de l'analyse du nuage des variables ne sont pas directement disponibles. Comment pouvez-vous les obtenir ? Affichez la projection du nuage des variables sur les deux premiers axes factoriels.

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

    .. only:: jupyter

        .. code-block:: python

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

    .. admonition:: Correction :
            
        Il suffit pour cela d'utiliser les relations de transition entre les deux analyses (voir le cours).

        .. code-block:: python
            
            mammalsCRTranspose = mammalsCR.transpose()
            irlambdas = 1/np.sqrt(pcaCR.explained_variance_)
            mirlambdas = np.diagflat(irlambdas)
            projectionsVars = (mammalsCRTranspose.dot(mNt)).dot(mirlambdas)/62
            nomsVars = np.genfromtxt('mammals.csv', dtype='str', delimiter=';', skip_header=0, max_rows=1)
            
            fig = plt.figure(figsize=(8,8))
            ax = fig.add_subplot(111)
            ax.set_xlim([-1, 1])
            ax.set_ylim([-1, 1])
            cercle = plt.Circle((0, 0), 1, color='b', fill=False)
            ax.add_artist(cercle)
            for i in range(len(nomsVars)-1):
                    x,y = projectionsVars[i,0], projectionsVars[i,1]
                    ax.scatter(x,y)
                    ax.arrow(0, 0, x, y, shape='full', lw=1, length_includes_head=True)
                    ax.text(x, y, nomsVars[i+1])
                
            plt.show()

..           .. code-block:: python
..            
..            mammalsCRTranspose = mammalsCR.transpose()
..            mammalsCR.shape
..            (62, 10)
..            mammalsCRTranspose.shape
..            (10, 62)
..            nomsVars = np.genfromtxt('mammals.csv', dtype='str', delimiter=';', skip_header=0, max_rows=1)
..            nomsVars
..            array(['Species', 'BodyW', 'BrainW', 'SWS', 'PS', 'TS', 'LifeSpan', 'GT',
..                   'PI', 'SEI', 'ODI'], 
..                  dtype='<U8')
..            pcaVars = PCA()
..            pcaVars.fit(mammalsCRTranspose)
..            mNTt = pcaVars.transform(mammalsCRTranspose)
..            fig = plt.subplots()
..            ax = fig.add_subplot(111)
..            cercle = plt.Circle((0, 0), np.sqrt(62), color='b', fill=False)
..            ax.add_artist(cercle)
..            for i in range(len(nomsVars)-1):
..                 x,y = mNTt[i,0],mNTt[i,1]
..                 ax.scatter(x,y)
..                 ax.text(x,y,nomsVars[i+1])
..                
..            plt.show()


ACP sur les données « feuilles »
********************************

Nous appliquerons maintenant l'ACP 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>`_. Pour cela, il est d'abord nécessaire d'enregistrer ces données dans le répertoire de travail. Dans une fenêtre terminal entrez (ou téléchargez depuis `cette adresse <http://cedric.cnam.fr/~crucianm/src/leaf.csv>`_:

.. only:: html
   
    .. code-block:: bash
    
        wget -nc http://cedric.cnam.fr/~crucianm/src/leaf.csv

.. only:: jupyter
   
    .. code-block:: bash
    
        !wget -nc http://cedric.cnam.fr/~crucianm/src/leaf.csv

Examinez le fichier téléchargé et regardez sur le site indiqué les explications concernant le format des données.

Appliquez l'ACP aux données et affichez la décroissance des valeurs propres :

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

    from sklearn import preprocessing
    leafN = preprocessing.scale(leaf[:,2:])
    pca = PCA()
    pca.fit(leafN)

    ratios = pca.explained_variance_ratio_

    plt.bar(range(len(ratios)), pca.explained_variance_ratio_)
    plt.xticks(range(len(ratios)))
    plt.xlabel("Composante principale")
    plt.ylabel("% de variance expliquée")
    plt.show()

.. admonition:: Question :

   Affichez les projections des données sur les 3 premiers axes principaux.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

        .. code-block:: python
            
            leafNt = pca.transform(leafN)
            leafNt[:2,:]
            
            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter(leafNt[:,0],leafNt[:,1],leafNt[:,2])
            plt.show()

.. admonition:: Question :

   Affichez les projections des données sur les 3 premiers axes principaux en utilisant l'étiquette de classe pour donner une couleur aux marques. L'étiquette de classe se trouve (sous forme numérique) dans la première colonne de ``leaf``. Il est nécessaire d'ajouter dans l'appel ``ax.scatter`` un argument pour indiquer que la couleur est obtenue à partir de données, ``c = ...``.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

        .. code-block:: python
            
            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter(leafNt[:,0],leafNt[:,1],leafNt[:,2],c=leaf[:,0])
            plt.show()

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

.. [AC76] Allison T., D. Cicchetti. *Sleep in mammals: ecological and constitutional correlates*. *Science*, 194: 732–734 (1976).