.. _td2acorrectionsession4Arst:

===================================================
2A.ml - Machine Learning et Marketting - correction
===================================================


.. only:: html

    **Links:** :download:`notebook <td2a_correction_session_4A.ipynb>`, :downloadlink:`html <td2a_correction_session_4A2html.html>`, :download:`python <td2a_correction_session_4A.py>`, :downloadlink:`slides <td2a_correction_session_4A.slides.html>`, :githublink:`GitHub|_doc/notebooks/td2a_ml/td2a_correction_session_4A.ipynb|*`


Classification binaire, correction.

.. code:: ipython3

    %matplotlib inline


.. parsed-literal::
    Populating the interactive namespace from numpy and matplotlib


.. code:: ipython3

    import matplotlib.pyplot as plt

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





Données
~~~~~~~

Tout d’abord, on récupère la base de données : `Bank Marketing Data
Set <https://archive.ics.uci.edu/ml/datasets/Bank+Marketing>`__.

.. code:: ipython3

    url = "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/"
    file = "bank.zip"
    import pyensae.datasource
    data = pyensae.datasource.download_data(file, website=url)

.. code:: ipython3

    import pandas
    df = pandas.read_csv("bank.csv",sep=";")
    df.tail()






.. raw:: html

    <div>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>age</th>
          <th>job</th>
          <th>marital</th>
          <th>education</th>
          <th>default</th>
          <th>balance</th>
          <th>housing</th>
          <th>loan</th>
          <th>contact</th>
          <th>day</th>
          <th>month</th>
          <th>duration</th>
          <th>campaign</th>
          <th>pdays</th>
          <th>previous</th>
          <th>poutcome</th>
          <th>y</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>4516</th>
          <td>33</td>
          <td>services</td>
          <td>married</td>
          <td>secondary</td>
          <td>no</td>
          <td>-333</td>
          <td>yes</td>
          <td>no</td>
          <td>cellular</td>
          <td>30</td>
          <td>jul</td>
          <td>329</td>
          <td>5</td>
          <td>-1</td>
          <td>0</td>
          <td>unknown</td>
          <td>no</td>
        </tr>
        <tr>
          <th>4517</th>
          <td>57</td>
          <td>self-employed</td>
          <td>married</td>
          <td>tertiary</td>
          <td>yes</td>
          <td>-3313</td>
          <td>yes</td>
          <td>yes</td>
          <td>unknown</td>
          <td>9</td>
          <td>may</td>
          <td>153</td>
          <td>1</td>
          <td>-1</td>
          <td>0</td>
          <td>unknown</td>
          <td>no</td>
        </tr>
        <tr>
          <th>4518</th>
          <td>57</td>
          <td>technician</td>
          <td>married</td>
          <td>secondary</td>
          <td>no</td>
          <td>295</td>
          <td>no</td>
          <td>no</td>
          <td>cellular</td>
          <td>19</td>
          <td>aug</td>
          <td>151</td>
          <td>11</td>
          <td>-1</td>
          <td>0</td>
          <td>unknown</td>
          <td>no</td>
        </tr>
        <tr>
          <th>4519</th>
          <td>28</td>
          <td>blue-collar</td>
          <td>married</td>
          <td>secondary</td>
          <td>no</td>
          <td>1137</td>
          <td>no</td>
          <td>no</td>
          <td>cellular</td>
          <td>6</td>
          <td>feb</td>
          <td>129</td>
          <td>4</td>
          <td>211</td>
          <td>3</td>
          <td>other</td>
          <td>no</td>
        </tr>
        <tr>
          <th>4520</th>
          <td>44</td>
          <td>entrepreneur</td>
          <td>single</td>
          <td>tertiary</td>
          <td>no</td>
          <td>1136</td>
          <td>yes</td>
          <td>yes</td>
          <td>cellular</td>
          <td>3</td>
          <td>apr</td>
          <td>345</td>
          <td>2</td>
          <td>249</td>
          <td>7</td>
          <td>other</td>
          <td>no</td>
        </tr>
      </tbody>
    </table>
    </div>



Exercice 1 : prédire y en fonction des attributs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Les données ne sont pas toutes au format numérique, il faut convertir
les variables catégorielles. Pour cela, on utilise la fonction
`DictVectorizer <http://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.DictVectorizer.html>`__.

.. code:: ipython3

    import numpy
    import numpy as np
    numerique = [ c for c,d in zip(df.columns,df.dtypes) if d == numpy.int64 ]
    categories = [ c for c in df.columns if c not in numerique and c not in ["y"] ]
    target = "y"
    print(numerique)
    print(categories)
    print(target)
    num = df[ numerique ]
    cat = df[ categories ]
    tar = df[ target ]


.. parsed-literal::
    ['age', 'balance', 'day', 'duration', 'campaign', 'pdays', 'previous']
    ['job', 'marital', 'education', 'default', 'housing', 'loan', 'contact', 'month', 'poutcome']
    y


On traite les variables catégorielles :

.. code:: ipython3

    from sklearn.feature_extraction import DictVectorizer
    prep = DictVectorizer()
    cat_as_dicts = [dict(r.iteritems()) for _, r in cat.iterrows()]
    temp = prep.fit_transform(cat_as_dicts)
    cat_exp = temp.toarray()
    prep.feature_names_




.. parsed-literal::
    ['contact=cellular',
     'contact=telephone',
     'contact=unknown',
     'default=no',
     'default=yes',
     'education=primary',
     'education=secondary',
     'education=tertiary',
     'education=unknown',
     'housing=no',
     'housing=yes',
     'job=admin.',
     'job=blue-collar',
     'job=entrepreneur',
     'job=housemaid',
     'job=management',
     'job=retired',
     'job=self-employed',
     'job=services',
     'job=student',
     'job=technician',
     'job=unemployed',
     'job=unknown',
     'loan=no',
     'loan=yes',
     'marital=divorced',
     'marital=married',
     'marital=single',
     'month=apr',
     'month=aug',
     'month=dec',
     'month=feb',
     'month=jan',
     'month=jul',
     'month=jun',
     'month=mar',
     'month=may',
     'month=nov',
     'month=oct',
     'month=sep',
     'poutcome=failure',
     'poutcome=other',
     'poutcome=success',
     'poutcome=unknown']



On construit les deux matrices :math:`(X,Y)` = (features, classe).

**Remarque :** certains modèles d’apprentissage n’acceptent pas les
corrélations. Lorsqu’on crée des variables catégorielles à choix unique,
les sommes des colonnes associées à une catégories fait nécessairement
un. Avec deux variables catégorielles, on introduit nécessairement des
corrélations. On pense à enlever les dernières catégories :
``'contact=unknown', 'default=yes', 'education=unknown',``
``'housing=yes',``
``'job=unknown', 'loan=yes', 'marital=single', 'month=sep',``
``'poutcome=unknown'``.

.. code:: ipython3

    cat_exp_df = pandas.DataFrame( cat_exp, columns = prep.feature_names_ )
    reject = ['contact=unknown', 'default=yes', 'education=unknown', 'housing=yes','job=unknown', 
              'loan=yes', 'marital=single', 'month=sep', 'poutcome=unknown']
    keep = [ c for c in cat_exp_df.columns if c not in reject ]
    cat_exp_df_nocor = cat_exp_df [ keep ]
    X = pandas.concat ( [ num, cat_exp_df_nocor ], axis= 1)
    Y = tar.apply( lambda r : (1.0 if r == "yes" else 0.0))
    X.shape, Y.shape




.. parsed-literal::
    ((4521, 42), (4521,))



Quelques corrélations sont très grandes malgré tout :

.. code:: ipython3

    import numpy
    numpy.corrcoef(X)




.. parsed-literal::
    array([[ 1.        ,  0.99733687,  0.96756321, ...,  0.89646012,
             0.9809934 ,  0.94772372],
           [ 0.99733687,  1.        ,  0.98202501, ...,  0.8917544 ,
             0.99153146,  0.96013301],
           [ 0.96756321,  0.98202501,  1.        , ...,  0.90132068,
             0.998155  ,  0.98697606],
           ..., 
           [ 0.89646012,  0.8917544 ,  0.90132068, ...,  1.        ,
             0.90491984,  0.94963638],
           [ 0.9809934 ,  0.99153146,  0.998155  , ...,  0.90491984,
             1.        ,  0.98331872],
           [ 0.94772372,  0.96013301,  0.98697606, ...,  0.94963638,
             0.98331872,  1.        ]])



On divise en base d’apprentissage et de test :

.. code:: ipython3

    from sklearn.model_selection import train_test_split
    X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.33)

Puis on cale un modèle d’apprentissage :

.. code:: ipython3

    from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
    type_classifier = GradientBoostingClassifier
    clf = type_classifier()
    clf = clf.fit(X_train, Y_train.ravel()) 

La méthode
`ravel <http://pandas.pydata.org/pandas-docs/stable/generated/pandas.Series.ravel.html#pandas.Series.ravel>`__
évite de prendre en compte l’index de *Y_train*. La méthode
`train_test_split <http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html>`__
conserve dans l’index les positions initiales des élèments. Mais l’index
fait que ``Y_train[0]`` ne désigne pas le premier élément de ``Y_train``
mais le premier élément du tableau initial. ``Y_train.ravel()[0]``
désigne bien le premier élément du tableau. On calcule ensuite la
matrice de confusion (`Confusion
matrix <http://scikit-learn.org/stable/auto_examples/plot_confusion_matrix.html>`__)
:

.. code:: ipython3

    from sklearn.metrics import confusion_matrix
    for x,y in [ (X_train, Y_train), (X_test, Y_test) ]:
        yp  = clf.predict(x)
        cm = confusion_matrix(y.ravel(), yp.ravel())
        print(cm)
        
    import matplotlib.pyplot as plt
    plt.matshow(cm)
    plt.title('Confusion matrix')
    plt.colorbar()
    plt.ylabel('True label')
    plt.xlabel('Predicted label') 


.. parsed-literal::
    [[2685   14]
     [ 161  169]]
    [[1261   40]
     [ 120   71]]




.. parsed-literal::

    <matplotlib.text.Text at 0x294ed68>




.. image:: td2a_correction_session_4A_21_2.png


Si le model choisi est un
`GradientBoostingClassifier <http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html>`__,
on peut regarder l’importance des variables dans la construction du
résultat. Le graphe suivant est inspiré de la page `Gradient Boosting
regression <http://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_regression.html>`__
même si ce n’est pas une régression qui a été utilisée ici.

.. code:: ipython3

    import numpy as np
    feature_name = X.columns
    
    limit = 20
    feature_importance = clf.feature_importances_[:20]
    feature_importance = 100.0 * (feature_importance / feature_importance.max())
    sorted_idx = np.argsort(feature_importance)
    pos = np.arange(sorted_idx.shape[0]) + .5
    plt.subplot(1, 2, 2)
    plt.barh(pos, feature_importance[sorted_idx], align='center')
    plt.yticks(pos, feature_name[sorted_idx])
    plt.xlabel('Relative Importance')
    plt.title('Variable Importance')




.. parsed-literal::
    <matplotlib.text.Text at 0x149d5a90>




.. image:: td2a_correction_session_4A_23_1.png


Il faut tout de même rester prudent quant à l’interprétation du graphe
précédent. La documentation au sujet de l\ `importance des
features <http://scikit-learn.org/stable/modules/ensemble.html#feature-importance-evaluation>`__
précise plus ou moins comment sont calculés ces chiffres. Toutefois,
lorsque des variables sont très corrélées, elles sont plus ou moins
interchangeables. Tout dépend alors comment l’algorithme d’apprentissage
choisit telle ou telle variables, toujours dans le même ordre ou dans un
ordre aléatoire.

variables
^^^^^^^^^

On utilise le code de la séance 3 `Analyse en Composantes
Principales <http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx/notebooks/td2a_cenonce_session_3A.html?highlight=acp#PCA>`__
pour observer les variables.

.. code:: ipython3

    from sklearn.decomposition import PCA
    pca = PCA(n_components=4)
    x_transpose = X.T
    pca.fit(x_transpose)




.. parsed-literal::
    PCA(copy=True, n_components=4, whiten=False)



.. code:: ipython3

    plt.bar(numpy.arange(len(pca.explained_variance_ratio_))+0.5, pca.explained_variance_ratio_)
    plt.title("Variance expliquée")




.. parsed-literal::
    <matplotlib.text.Text at 0x161b9630>




.. image:: td2a_correction_session_4A_27_1.png


.. code:: ipython3

    import warnings
    warnings.filterwarnings('ignore')
    X_reduced = pca.transform(x_transpose)
    plt.figure(figsize=(18,6))
    plt.scatter(X_reduced[:, 0], X_reduced[:, 1])
    for label, x, y in zip(x_transpose.index, X_reduced[:, 0], X_reduced[:, 1]):
        plt.annotate(
            label,
            xy = (x, y), xytext = (-10, 10),
            textcoords = 'offset points', ha = 'right', va = 'bottom',
            bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),
            arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))



.. image:: td2a_correction_session_4A_28_0.png


Les variables les plus dissemblables sont celles qui contribuent le
plus. Toutefois, à la vue de ce graphique, il apparaît qu’il faut
normaliser les données avant d’interpréter l’ACP :

.. code:: ipython3

    from sklearn.preprocessing import normalize
    xnorm = normalize(x_transpose)
    pca = PCA(n_components=10)
    pca.fit(xnorm)
    plt.bar(numpy.arange(len(pca.explained_variance_ratio_))+0.5, pca.explained_variance_ratio_)
    plt.title("Variance expliquée")




.. parsed-literal::
    <matplotlib.text.Text at 0x1628bd68>




.. image:: td2a_correction_session_4A_30_1.png


.. code:: ipython3

    X_reduced = pca.transform(xnorm)
    plt.figure(figsize=(18,6))
    plt.scatter(X_reduced[:, 0], X_reduced[:, 1])
    for label, x, y in zip(x_transpose.index, X_reduced[:, 0], X_reduced[:, 1]):
        plt.annotate(
            label,
            xy = (x, y), xytext = (-10, 10),
            textcoords = 'offset points', ha = 'right', va = 'bottom',
            bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),
            arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))



.. image:: td2a_correction_session_4A_31_0.png


Nettement mieux. En règle générale, il est préférable de normaliser ses
données avant d’apprendre un modèle. Cela n’est pas toujours nécessaire
(comme pour les arbres de décision). Toutefois, numériquement, avoir des
données d’ordre de grandeur très différent introduit toujours plus
d’approximations.

.. code:: ipython3

    from sklearn.pipeline import Pipeline
    from sklearn.preprocessing import Normalizer
    from sklearn.ensemble import GradientBoostingClassifier
    
    clf = Pipeline([
          ('normalize', Normalizer()),
          ('classification', GradientBoostingClassifier())
        ])
    clf = clf.fit(X_train, Y_train.ravel()) 

.. code:: ipython3

    from sklearn.metrics import confusion_matrix
    x,y = X_test, Y_test
    yp  = clf.predict(x)
    cm2 = confusion_matrix(y, yp)
    print("non normalisé\n",cm)
    print("normalisé\n",cm2)


.. parsed-literal::
    non normalisé
     [[1261   40]
     [ 120   71]]
    normalisé
     [[1270   31]
     [ 127   64]]


C’est plus ou moins équivalent lorsque les variables sont normalisées
dans ce cas. Il faudrait vérifier sur la courbe ROC.

Exercice 2 : tracer la courbe ROC
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On utilise l’exemple `Receiver Operating Characteristic
(ROC) <http://scikit-learn.org/stable/auto_examples/plot_roc.html>`__
qu’il faut modifié car la réponse juste dans notre cas est le fait de
prédire la bonne classe. Cela veut dire qu’il y a deux cas pour lesquels
le modèle prédit le bon résultat : on choisit la classe qui la
probabilité la plus forte.

.. code:: ipython3

    from sklearn.metrics import roc_curve, auc
    probas = clf.predict_proba(X_test)
    probas[:5]




.. parsed-literal::
    array([[ 0.9792842 ,  0.0207158 ],
           [ 0.98228488,  0.01771512],
           [ 0.83708208,  0.16291792],
           [ 0.98675536,  0.01324464],
           [ 0.98818146,  0.01181854]])



On construit le vecteur des bonnes réponses :

.. code:: ipython3

    rep = [ ]
    yt = Y_test.ravel()
    for i in range(probas.shape[0]):
        p0,p1 = probas[i,:]
        exp = yt[i]
        if p0 > p1 :
            if exp == 0 :
                # bonne réponse
                rep.append ( (1, p0) )
            else :
                # mauvaise réponse
                rep.append( (0,p0) )
        else :
            if exp == 0 :
                # mauvaise réponse
                rep.append ( (0, p1) )
            else :
                # bonne réponse
                rep.append( (1,p1) )
    mat_rep = numpy.array(rep)
    mat_rep[:5]




.. parsed-literal::
    array([[ 1.        ,  0.9792842 ],
           [ 1.        ,  0.98228488],
           [ 1.        ,  0.83708208],
           [ 1.        ,  0.98675536],
           [ 1.        ,  0.98818146]])



.. code:: ipython3

    "taux de bonne réponse",sum(mat_rep[:,0]/len(mat_rep)) # voir matrice de confusion




.. parsed-literal::
    ('taux de bonne réponse', 0.89410187667560337)



.. code:: ipython3

    fpr, tpr, thresholds = roc_curve(mat_rep[:,0], mat_rep[:, 1])
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.0])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC')
    plt.legend(loc="lower right")




.. parsed-literal::
    <matplotlib.legend.Legend at 0x17e0df60>




.. image:: td2a_correction_session_4A_41_1.png


Ce score n’est pas si mal pour un premier essai. On n’a pas tenu compte
du fait que la classe 1 est sous-représentée (voir `Quelques astuces
pour faire du machine
learning <http://www.xavierdupre.fr/blog/2014-03-28_nojs.html>`__. A
priori, ce ne devrait pas être le cas du
`GradientBoostingClassifier <http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html>`__.
C’est une famille de modèles qui, lors de l’apprentissage, pondère
davantage les exemples où ils font des erreurs. L’algorithme de
`boosting <http://fr.wikipedia.org/wiki/Boosting>`__ le plus connu est
`AdaBoost <http://fr.wikipedia.org/wiki/AdaBoost>`__.

On tire maintenant deux échantillons aléatoires qu’on ajoute au
graphique précédent :

.. code:: ipython3

    import random
    Y1 = numpy.array([ random.randint(0,1) == 0 for i in range(0,mat_rep.shape[0]) ])
    Y2 = numpy.array([ random.randint(0,1) == 0 for i in range(0,mat_rep.shape[0]) ])
    
    fpr1, tpr1, thresholds1 = roc_curve(mat_rep[Y1,0], mat_rep[Y1, 1])
    roc_auc1 = auc(fpr1, tpr1)
    fpr2, tpr2, thresholds2 = roc_curve(mat_rep[Y2,0], mat_rep[Y2, 1])
    roc_auc2 = auc(fpr2, tpr2)
    print(fpr1.shape,tpr1.shape,fpr2.shape,tpr2.shape)
    
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    
    ax.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
    ax.plot([0, 1,2], [0, 1,2], 'k--')
    ax.set_xlim([0.0, 1.0])
    ax.set_ylim([0.0, 1.0])
    ax.set_xlabel('False Positive Rate')
    ax.set_ylabel('True Positive Rate')
    ax.set_title('ROC')
    ax.plot(fpr1, tpr1, label='ech 1, area=%0.2f' % roc_auc1)
    ax.plot(fpr2, tpr2, label='ech 2, area=%0.2f' % roc_auc2)
    ax.legend(loc="lower right")


.. parsed-literal::
    (697,) (697,) (699,) (699,)




.. parsed-literal::

    <matplotlib.legend.Legend at 0x17e82940>




.. image:: td2a_correction_session_4A_43_2.png