.. _chap-tpEstimationDensiteNoyaux:

####################################################
Travaux pratiques - Estimation de densité par noyaux
####################################################

(1 séance de TP)

.. only:: html

    .. container:: notebook

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

        `Cahier Jupyter <jupyter/tpEstimationDensiteNoyaux.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'estimation de densité par noyaux, ainsi que de contribuer à une meilleure compréhension de cette méthode et de l'impact de la distribution des données sur les résultats. Pour cela, sont d'abord examinées des données générées de façon contrôlée et ensuite des données réelles vues en cours.


L'estimation de densités dans Scikit-learn
******************************************

Pour estimer des densités, Scikit-learn propose les histogrammes et les noyaux comme méthodes non paramétriques, et les modèles de mélange comme méthodes paramétriques. Manquent, en revanche, des outils de choix du paramètre unique (largeur de « fenêtre ») pour les méthodes non paramétriques et des outils de comparaison entre densités (même lorsque celles-ci sont issues de mélanges de lois normales).

Les méthodes non paramétriques proposées sont présentées dans `sklearn.density <http://scikit-learn.org/stable/modules/density.html>`_ et l'utilisation de modèles de mélange dans `sklearn.mixture <http://scikit-learn.org/stable/modules/mixture.html>`_.

Estimation par noyaux
*********************

Afin de permettre la visualisation des résultats et sachant que des outils de comparaison entre densités ne sont pas directement accessibles, nous examinerons uniquement des données unidimensionnelles et bidimensionnelles.

La description de l'implémentation de l'estimation de densités par noyaux se trouve dans `KernelDensity <http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html>`_.

Le type de noyau est spécifié dans le paramètre ``kernel`` est il est possible d'utiliser comme valeurs ``'gaussian'`` (choix par défaut), ``'tophat'``, ``'epanechnikov'``, ``'exponential'``, ``'linear'`` ou ``'cosine'``. Une fois le type de noyau choisi, tous les noyaux (ou « fenêtres ») employés présentent une même valeur pour l'unique paramètre de largeur de « fenêtre » appelé ici ``bandwidth``.

Il faut noter que le coût du calcul de la densité en un point donné peut être élevé si l'ensemble d'observations à partir duquel l'estimation de densité a été réalisée est très volumineux (en effet, cela impose de calculer la valeur du noyau entre le point donné et chacune des observations de départ !). Pour réduire le coût de calcul, on observe que les noyaux sont « localisés », c'est à dire leur valeur est numériquement non nulle seulement pour une valeur relativement faible de la norme de la différence entre les arguments ; la valeur du noyau calculé entre un point et une observation très éloignée est ainsi assimilable à 0 et il n'est donc pas nécessaire de faire le calcul. Il est alors possible d'utiliser un index spatial (ici de type *KD tree* ou *Ball tree*, ``algorithm ['kd_tree'|'ball_tree'|'auto']``) pour éviter de calculer le noyau entre données trop éloignées les unes des autres.

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

* ``fit(X, y=None)`` : calcul du modèle à partir des observations qui sont les lignes de X.
* ``sample([n_samples, random_state])`` : générer un échantillon (composé de ``n_samples`` données) à partir du modèle (pour le moment, accessible seulement pour les noyaux ``'gaussian'`` et ``'tophat'``).
* ``score(X, y=None)`` : retourne la log-vraisemblance totale des données de X par rapport au modèle.
* ``score_samples(X)`` : retourne le logarithme de la densité calculée pour chaque donnée de X.
* ``get_params([deep])`` : lire les valeurs des paramètres de l'estimateur employé.
* ``set_params(**params)`` : donner des valeurs aux paramètres de l'estimateur employé.

Il n'y a pas d'attributs publics, pourquoi ?


Estimation à partir de données générées
---------------------------------------

Cette partie de la séance a pour objectif de faire comprendre l'impact des différents aspects qui interviennent dans l'estimation par noyaux : le volume de données, le choix du noyau, le paramètre d'échelle (``bandwidth``, correspondant en général à la variance du noyau). Pour cela, nous générons des données à partir d'une densité connue, estimons la densité à partir de ces données et comparons le résultat de l'estimation à la densité qui a généré les données.

Une première utilisation sur des données unidimensionnelles, sur lesquelles il est possible d'afficher à la fois la vraie densité et la densité estimée :

.. code-block:: python
    
    import numpy as np
    import matplotlib.pyplot as plt

    from scipy.stats import norm
    from sklearn.neighbors import KernelDensity
    
    # générer l'échantillon à partir de deux lois normales
    N = 100
    X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),
                        np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
    
    # préparer les points où on calculera la densité
    X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
    
    # préparation de l'affichage de la vraie densité, qui est celle à partir
    #  de laquelle les données ont été générées (voir plus haut)
    # la pondération des lois dans la somme est la pondération des lois
    #  dans l'échantillon généré (voir plus haut)
    true_density = (0.3 * norm(0,1).pdf(X_plot[:,0]) + 0.7 * norm(5,1).pdf(X_plot[:,0]))

La variable ``true_density`` correspond à la distribution réelle de nos données. En pratique, celle-ci est inconnue et nous cherchons à l'estimer.
Pour ce faire, nous allons utiliser la méthode des noyaux.

.. code-block:: python
    
    # estimation de densité par noyaux gaussiens
    kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(X)

Nous pouvons ensuite calculer, pour chaque valeur de ``X_plot``, la probabilité d'occurrence de cette valeur. Cela permet de tracer la densité estimée.

.. code-block:: python

    # calcul de la densité pour les données de X_plot
    density = np.exp(kde.score_samples(X_plot))
    
    # affichage : vraie densité et estimation
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.fill(X_plot[:,0], true_density, fc='b', alpha=0.2, label='Vraie densité')
    ax.plot(X_plot[:,0], density, '-', label="Estimation")
    ax.plot(X[:, 0], -0.005 - 0.01 * np.random.random(X.shape[0]), '+k')
    ax.legend(loc='upper left')
    plt.show()


.. admonition:: Question :

   Faites varier la taille de l'échantillon et examinez les résultats.

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

   .. admonition:: Correction :

        Pour des échantillons plus restreints (``N`` plus petit) la densité estimée s'écarte en général plus de la densité réélle, pour des échantillons plus larges (``N`` plus grand) l'écart entre les deux densités est en général plus réduit.

.. admonition:: Question :

   Faites varier le nombre de lois normales lors de la génération de l'échantillon et examinez les résultats.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

        Il suffit d'utiliser plus de lois lors de la génération de l'échantillon, en prenant soin de les « localiser » à des abscisses différentes (utiliser des valeurs différentes pour les espérances). Par ex. 3 lois :
           
        .. code-block:: python
            
            X = np.concatenate((np.random.normal(0, 1, 0.3 * N),
                                np.random.normal(5, 1, 0.4 * N),
                                np.random.normal(8, 1, 0.3 * N)))[:, np.newaxis]
           
        Plus le nombre de lois utilisé pour générer les données est élevé, plus la taille de l'échantillon doit être élevée pour que l'estimation reste de bonne qualité.

.. admonition:: Question :

   Faites varier la valeur de ``bandwidth`` et examinez les résultats.

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

   .. admonition:: Correction :

         Une valeur de ``bandwidth`` trop faible risque d'augmenter le nombre de « modes » (pics différents) dans l'estimation, surtout si l'échantillon est restreint. Une valeur trop élevée pour ``bandwidth`` produit un « lissage » excessiv, par ex. fusionne les « modes » estimés.

.. admonition:: Question :

   Faites varier le type de noyau et examinez les résultats.

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

   .. admonition:: Correction :

         Le type de noyau a un impact sur la forme de la densité estimée, mais l'impact sur la qualité d'estimation est plus faible que celui de la valeur de ``bandwidth``.

Il est possible de générer de nouvelles données suivant cette densité estimée et de visualiser le résultat d'une nouvelle estimation faite à partir de ces dernières données :

.. code-block:: python
    
    Xg = kde.sample(N)
    kde2 = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(Xg)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.fill(X_plot[:,0], true_density, fc='b', alpha=0.2, label='Vraie densité')
    ax.plot(X_plot[:,0], density, '-', label="Estimation")
    ax.plot(X_plot[:,0], np.exp(kde2.score_samples(X_plot)), 'r-', label="Estimation2")
    ax.legend(loc='upper left')
    plt.show()


Une deuxième utilisation de l'estimation de densité par noyaux sur des données bidimensionnelles, sur lesquelles il n'est malheureusement plus possible de bien visualiser à la fois la vraie densité et la densité estimée :

.. code-block:: python
    
    # générer l'échantillon
    N = 20
    kd = np.random.rand(N, 2)
    
    # définir la grille pour la visualisation
    grid_size = 100
    Gx = np.arange(0, 1, 1/grid_size)
    Gy = np.arange(0, 1, 1/grid_size)
    Gx, Gy = np.meshgrid(Gx, Gy)
    
    # définir la largeur de bande pour le noyau
    bw = 0.02
    
    # estimation, puis calcul densité sur la grille
    kde3 = KernelDensity(kernel='gaussian', bandwidth=bw).fit(kd)
    Z = np.exp(kde3.score_samples(np.hstack(((Gx.reshape(grid_size*grid_size))[:,np.newaxis],
            (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))
    
    # affichage
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(Gx, Gy, Z.reshape(grid_size,grid_size), rstride=1, 
                        cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
    ax.scatter(kd[:,0], kd[:,1], -10)
    plt.show()

.. admonition:: Question :

   Faites varier la valeur de ``bandwidth`` et examinez les résultats.

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

   .. admonition:: Correction :

         Nous pouvons tirer les mêmes conclusions que pour le cas unidimensionnel.

.. admonition:: Question :

   Générez l'échantillon à partir d'une somme de 2 lois normales différentes (avec ``N`` plus grand) et examinez les résultats.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

         Les données avaient été générées suivant une loi uniforme dans l'intervalle :math:`[0, 1)^2` et la grille de visualisation définie en fonction de ce choix. Pour générer et visualiser les données à partir d'une somme de 2 lois normales il faut non seulement modifier la génération des données, mais aussi la définition de la grille de visualisation.
           
           .. code-block:: python
            
            # générer l'échantillon
            N = 100
            kd1 = np.random.randn(N/2, 2) + (1.0, 1.0)
            kd2 = np.random.randn(N/2, 2)
            kd = np.concatenate((kd1, kd2))
            
            # définir la grille pour la visualisation
            grid_size = 100
            Gx = np.arange(-1, 3, 4/grid_size)
            Gy = np.arange(-1, 3, 4/grid_size)
            Gx, Gy = np.meshgrid(Gx, Gy)
            
            # définir la largeur de bande pour le noyau
            bw = 0.2
            
            # estimation, puis calcul densité sur la grille
            kde4 = KernelDensity(kernel='gaussian', bandwidth=bw).fit(kd)
            Z = np.exp(kde4.score_samples(np.hstack(((Gx.reshape(grid_size * 
                                                                     grid_size))[:,np.newaxis],
                           (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))
            
            # affichage
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            ax.plot_surface(Gx, Gy, Z.reshape(grid_size,grid_size), rstride=1, 
                                cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
            ax.scatter(kd[:,0], kd[:,1], -10)
            plt.show()


Pour mieux comprendre ce qu'est une estimation de densité par noyaux, il peut être utile de construire le modèle et visualiser les résultats directement, sans passer par la classe spécifique de scikit-learn. En effet, dans l'estimation de densité par noyaux il n'y a pas de « modèle synthétique » obtenu à partir des N données (ou observations) de départ, le modèle de densité est constitué par les N noyaux centrés sur ces N données.

Pour la visualisation directe nous nous servirons de la fonction de calcul d'une loi normale bidimensionelle en ``(x,y)``, ``bivariate_normal(x,y,sigmax,sigmay,mux,muy,sigmaxy)``, où ``sigmax`` est la variance de X, ``sigmay`` la variance de Y, ``sigmaxy`` la covariance entre X et Y, ``mux`` la moyenne de X et ``muy`` la moyenne de Y.

.. code-block:: python
    
    # générer l'échantillon
    N = 10
    kd = np.random.rand(N, 2)
    
    # définir la grille pour la visualisation
    grid_size = 100
    Gx = np.arange(0, 1, 1/grid_size)
    Gy = np.arange(0, 1, 1/grid_size)
    Gx, Gy = np.meshgrid(Gx, Gy)
    bw = 0.01
    
    # calcul direct de la densité sur la grille
    from scipy.stats import multivariate_normal
    Z = np.zeros((grid_size,grid_size))
    for i in range(N):
        # génération d'une loi normale bidimensionnelle centrée sur le point kd[i] et d'écart-type bw
        Z = Z + multivariate_normal.pdf(np.dstack((Gx, Gy)),
                                        mean=[kd[i,0], kd[i,1]], cov=[[bw**2, 0], [0, bw**2]])
    
    # affichage
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(Gx, Gy, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                        linewidth=0, antialiased=True)
    ax.scatter(kd[:,0], kd[:,1], -500)
    plt.show()

.. admonition:: Question :

   Faites varier la taille de l'échantillon et examinez les résultats. Quelles sont les conséquences d'une augmentation de ``N`` ?

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

   .. admonition:: Correction :

         L'augmentation de la taille de l'échantillon avec une valeur inchangée pour ``bw`` produit une estimation plus lisse.


Estimation à partir des données « textures »
--------------------------------------------

Pour rappel, ces données correspondent à 5500 observations décrites par 40 variables. Chaque observation appartient à une des 11 classes de textures ; chaque classe est représentée par 500 observations. 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/>`_. Nous appliquerons l'estimation de densité par noyaux à ces données.

Si les données ne sont pas dans le répertoire de travail il est nécessaire de les chercher d'abord :

.. code-block:: bash
    
    %%bash
    wget -nc http://cedric.cnam.fr/~crucianm/src/texture.dat

Pour visualiser facilement les résultats, il faudra travailler sur des projections unidimensionnelles ou bidimensionnelles de ces données.

.. admonition:: Question :

   Appliquez l'analyse en composantes principales à ces données et estimez la densité des données projetées sur le premier axe principal (ou les deux premiers axes principaux). Faites varier la largeur de « fenêtre » (``bandwidth``). Visualisez les résultats.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

        Le programme suivant fait l'estimation sur les 2 premiers axes principaux :
           
        .. code-block:: python
            
            # lecture des données
            textures = np.loadtxt('texture.dat')

            from sklearn.decomposition import PCA
            pca = PCA(n_components=2)
            pca.fit(textures[:,:40])
            texturesp = pca.transform(textures[:,:40])
            
            # définir la grille pour la visualisation
            grid_size = 100
            mx = min(texturesp[:,0])
            Mx = max(texturesp[:,0])
            my = min(texturesp[:,1])
            My = max(texturesp[:,1])
            xstep = (Mx - mx) / grid_size
            ystep = (My - my) / grid_size
            Gx = np.arange(mx, Mx, xstep)
            Gy = np.arange(my, My, ystep)
            Gx, Gy = np.meshgrid(Gx, Gy)
            
            # définir la largeur de bande pour le noyau
            bw = 0.06
            
            # estimation, puis calcul densité sur la grille
            kdet = KernelDensity(kernel='gaussian', bandwidth=bw).fit(texturesp)
            Z = np.exp(kdet.score_samples(np.hstack((
                                (Gx.reshape(grid_size*grid_size))[:,np.newaxis],
                                (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))
            
            # affichage
            from mpl_toolkits.mplot3d import Axes3D
            from matplotlib import cm
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            ax.plot_surface(Gx, Gy, Z.reshape(grid_size,grid_size), rstride=1, 
                                cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
            ax.scatter(texturesp[:,0], texturesp[:,1], -0.5)
            plt.show()

.. admonition:: Question :

   Appliquez l'analyse discriminante à ces données et estimez la densité des données projetées sur le premier axe discriminant (ou les deux premiers axes discriminants). Faites varier la largeur de « fenêtre » (``bandwidth``). Visualisez les résultats.

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

    .. only:: jupyter

        .. code-block:: python


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

   .. admonition:: Correction :

        Le programme suivant fait l'estimation sur les 2 premiers axes discriminants :
           
        .. code-block:: python
            
            # lecture des données
            textures = np.loadtxt('texture.dat')

            from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
            lda = LinearDiscriminantAnalysis(n_components=2)
            lda.fit(textures[:,:40],textures[:,40])
            texturest = lda.transform(textures[:,:40])
            
            # définir la grille pour la visualisation
            grid_size = 100
            mx = min(texturest[:,0])
            Mx = max(texturest[:,0])
            my = min(texturest[:,1])
            My = max(texturest[:,1])
            xstep = (Mx - mx) / grid_size
            ystep = (My - my) / grid_size
            Gx = np.arange(mx, Mx, xstep)
            Gy = np.arange(my, My, ystep)
            Gx, Gy = np.meshgrid(Gx, Gy)
            
            # définir la largeur de bande pour le noyau
            bw = 0.25
            
            # estimation, puis calcul densité sur la grille
            kdet = KernelDensity(kernel='gaussian', bandwidth=bw).fit(texturest)
            Z = np.exp(kdet.score_samples(np.hstack((
                                (Gx.reshape(grid_size*grid_size))[:,np.newaxis],
                                (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))
            
            # affichage
            from mpl_toolkits.mplot3d import Axes3D
            from matplotlib import cm
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            ax.plot_surface(Gx, Gy, Z.reshape(grid_size,grid_size), rstride=1, 
                                cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
            ax.scatter(texturest[:,0], texturest[:,1], -0.05)
            plt.show()