.. _2017session6rst:

====================================
Session 26/6/2017 - machine learning
====================================


.. only:: html

    **Links:** :download:`notebook <2017_session6.ipynb>`, :downloadlink:`html <2017_session62html.html>`, :download:`PDF <2017_session6.pdf>`, :download:`python <2017_session6.py>`, :downloadlink:`slides <2017_session6.slides.html>`, :githublink:`GitHub|_doc/notebooks/sessions/2017_session6.ipynb|*`


Découverte des trois problèmes de machine learning exposé dans l’article
`Machine Learning - session
6 <http://www.xavierdupre.fr/app/actuariat_python/helpsphinx/blog/2017/2017-06-25_plan_session6.html>`__.

.. code:: 

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





Problème 1 : comparaison random forest, linéaire
------------------------------------------------

C’est un problème de régression. On cherche à comparer une random forest
avec un modèle linéaire.


-  Comparaison des tests de coefficients pour un modèle linéaire
   `OLS <http://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.html%3E>`__
   et des `features
   importance <http://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_importances.html>`__
-  Résultat au niveau d’une observation
   `treeinterpreter <https://github.com/andosa/treeinterpreter%3E>`__
-  Données :
   `Housing <https://archive.ics.uci.edu/ml/datasets/housing%3E>`__,
   `Forest
   Fire <https://archive.ics.uci.edu/ml/datasets/Forest+Fires>`__

Données
~~~~~~~

.. code:: 

    import pandas
    df = pandas.read_csv("data/housing.data", delim_whitespace=True, header=None)
    df.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>0</th>
          <th>1</th>
          <th>2</th>
          <th>3</th>
          <th>4</th>
          <th>5</th>
          <th>6</th>
          <th>7</th>
          <th>8</th>
          <th>9</th>
          <th>10</th>
          <th>11</th>
          <th>12</th>
          <th>13</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.00632</td>
          <td>18.0</td>
          <td>2.31</td>
          <td>0</td>
          <td>0.538</td>
          <td>6.575</td>
          <td>65.2</td>
          <td>4.0900</td>
          <td>1</td>
          <td>296.0</td>
          <td>15.3</td>
          <td>396.90</td>
          <td>4.98</td>
          <td>24.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.02731</td>
          <td>0.0</td>
          <td>7.07</td>
          <td>0</td>
          <td>0.469</td>
          <td>6.421</td>
          <td>78.9</td>
          <td>4.9671</td>
          <td>2</td>
          <td>242.0</td>
          <td>17.8</td>
          <td>396.90</td>
          <td>9.14</td>
          <td>21.6</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.02729</td>
          <td>0.0</td>
          <td>7.07</td>
          <td>0</td>
          <td>0.469</td>
          <td>7.185</td>
          <td>61.1</td>
          <td>4.9671</td>
          <td>2</td>
          <td>242.0</td>
          <td>17.8</td>
          <td>392.83</td>
          <td>4.03</td>
          <td>34.7</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.03237</td>
          <td>0.0</td>
          <td>2.18</td>
          <td>0</td>
          <td>0.458</td>
          <td>6.998</td>
          <td>45.8</td>
          <td>6.0622</td>
          <td>3</td>
          <td>222.0</td>
          <td>18.7</td>
          <td>394.63</td>
          <td>2.94</td>
          <td>33.4</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.06905</td>
          <td>0.0</td>
          <td>2.18</td>
          <td>0</td>
          <td>0.458</td>
          <td>7.147</td>
          <td>54.2</td>
          <td>6.0622</td>
          <td>3</td>
          <td>222.0</td>
          <td>18.7</td>
          <td>396.90</td>
          <td>5.33</td>
          <td>36.2</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    cols = "CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT MEDV".split()
    df.columns = cols
    df.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>CRIM</th>
          <th>ZN</th>
          <th>INDUS</th>
          <th>CHAS</th>
          <th>NOX</th>
          <th>RM</th>
          <th>AGE</th>
          <th>DIS</th>
          <th>RAD</th>
          <th>TAX</th>
          <th>PTRATIO</th>
          <th>B</th>
          <th>LSTAT</th>
          <th>MEDV</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.00632</td>
          <td>18.0</td>
          <td>2.31</td>
          <td>0</td>
          <td>0.538</td>
          <td>6.575</td>
          <td>65.2</td>
          <td>4.0900</td>
          <td>1</td>
          <td>296.0</td>
          <td>15.3</td>
          <td>396.90</td>
          <td>4.98</td>
          <td>24.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.02731</td>
          <td>0.0</td>
          <td>7.07</td>
          <td>0</td>
          <td>0.469</td>
          <td>6.421</td>
          <td>78.9</td>
          <td>4.9671</td>
          <td>2</td>
          <td>242.0</td>
          <td>17.8</td>
          <td>396.90</td>
          <td>9.14</td>
          <td>21.6</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.02729</td>
          <td>0.0</td>
          <td>7.07</td>
          <td>0</td>
          <td>0.469</td>
          <td>7.185</td>
          <td>61.1</td>
          <td>4.9671</td>
          <td>2</td>
          <td>242.0</td>
          <td>17.8</td>
          <td>392.83</td>
          <td>4.03</td>
          <td>34.7</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.03237</td>
          <td>0.0</td>
          <td>2.18</td>
          <td>0</td>
          <td>0.458</td>
          <td>6.998</td>
          <td>45.8</td>
          <td>6.0622</td>
          <td>3</td>
          <td>222.0</td>
          <td>18.7</td>
          <td>394.63</td>
          <td>2.94</td>
          <td>33.4</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.06905</td>
          <td>0.0</td>
          <td>2.18</td>
          <td>0</td>
          <td>0.458</td>
          <td>7.147</td>
          <td>54.2</td>
          <td>6.0622</td>
          <td>3</td>
          <td>222.0</td>
          <td>18.7</td>
          <td>396.90</td>
          <td>5.33</td>
          <td>36.2</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    X = df.drop("MEDV", axis=1)
    y = df["MEDV"]
    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)

Random Forest
~~~~~~~~~~~~~

.. code:: 

    from sklearn.ensemble import RandomForestRegressor
    clr = RandomForestRegressor()
    clr.fit(X, y)




.. parsed-literal::
    RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
                          max_features='auto', max_leaf_nodes=None,
                          min_impurity_decrease=0.0, min_impurity_split=None,
                          min_samples_leaf=1, min_samples_split=2,
                          min_weight_fraction_leaf=0.0, n_estimators=100,
                          n_jobs=None, oob_score=False, random_state=None,
                          verbose=0, warm_start=False)



.. code:: 

    importances = clr.feature_importances_
    importances




.. parsed-literal::
    array([0.04377988, 0.00098713, 0.00669368, 0.00092023, 0.02074126,
           0.4076033 , 0.01283279, 0.06458251, 0.0036521 , 0.01366027,
           0.01576497, 0.01209152, 0.39669034])



On s’inspire de l’exemple `Feature importances with forests of
trees <http://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_importances.html>`__.

.. code:: 

    %matplotlib inline

.. code:: 

    import matplotlib.pyplot as plt
    import numpy as np
    plt.figure(figsize=(12,4))
    indices = np.argsort(importances)[::-1]
    std = np.std([tree.feature_importances_ for tree in clr.estimators_],
                 axis=0)
    plt.title("Feature importances")
    plt.bar(range(X.shape[1]), importances[indices],
            color="r", yerr=std[indices], align="center")
    xlabels = list(df.columns[:-1])
    xlabels = [xlabels[i] for i in indices]
    plt.xticks(range(X.shape[1]), xlabels)
    plt.xlim([-1, X.shape[1]])
    plt.show()



.. image:: 2017_session6_12_0.png


.. code:: 

    from sklearn.metrics import r2_score
    r2_score(y_train, clr.predict(X_train))




.. parsed-literal::
    0.979236135448983



.. code:: 

    r2_score(y_test, clr.predict(X_test))




.. parsed-literal::
    0.9843720628811461



Modèle linéaire
~~~~~~~~~~~~~~~

.. code:: 

    import statsmodels.api as sm

.. code:: 

    model = sm.OLS(y_train, X_train)

.. code:: 

    results = model.fit()

.. code:: 

    results.params




.. parsed-literal::
    CRIM      -0.109867
    ZN         0.045855
    INDUS     -0.054801
    CHAS       3.758792
    NOX       -5.285538
    RM         6.267853
    AGE       -0.022554
    DIS       -1.169496
    RAD        0.221785
    TAX       -0.010410
    PTRATIO   -0.447963
    B          0.016474
    LSTAT     -0.259210
    dtype: float64



.. code:: 

    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>          <td>MEDV</td>       <th>  R-squared:         </th> <td>   0.958</td> 
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.957</td> 
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   575.3</td> 
    </tr>
    <tr>
      <th>Date:</th>             <td>Fri, 07 Jun 2019</td> <th>  Prob (F-statistic):</th> <td>7.10e-216</td>
    </tr>
    <tr>
      <th>Time:</th>                 <td>16:38:39</td>     <th>  Log-Likelihood:    </th> <td> -1025.5</td> 
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>   339</td>      <th>  AIC:               </th> <td>   2077.</td> 
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   326</td>      <th>  BIC:               </th> <td>   2127.</td> 
    </tr>
    <tr>
      <th>Df Model:</th>              <td>    13</td>      <th>                     </th>     <td> </td>    
    </tr>
    <tr>
      <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    
    </tr>
    </table>
    <table class="simpletable">
    <tr>
         <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
    </tr>
    <tr>
      <th>CRIM</th>    <td>   -0.1099</td> <td>    0.039</td> <td>   -2.845</td> <td> 0.005</td> <td>   -0.186</td> <td>   -0.034</td>
    </tr>
    <tr>
      <th>ZN</th>      <td>    0.0459</td> <td>    0.018</td> <td>    2.543</td> <td> 0.011</td> <td>    0.010</td> <td>    0.081</td>
    </tr>
    <tr>
      <th>INDUS</th>   <td>   -0.0548</td> <td>    0.079</td> <td>   -0.697</td> <td> 0.486</td> <td>   -0.209</td> <td>    0.100</td>
    </tr>
    <tr>
      <th>CHAS</th>    <td>    3.7588</td> <td>    1.116</td> <td>    3.370</td> <td> 0.001</td> <td>    1.564</td> <td>    5.953</td>
    </tr>
    <tr>
      <th>NOX</th>     <td>   -5.2855</td> <td>    4.398</td> <td>   -1.202</td> <td> 0.230</td> <td>  -13.938</td> <td>    3.367</td>
    </tr>
    <tr>
      <th>RM</th>      <td>    6.2679</td> <td>    0.405</td> <td>   15.461</td> <td> 0.000</td> <td>    5.470</td> <td>    7.065</td>
    </tr>
    <tr>
      <th>AGE</th>     <td>   -0.0226</td> <td>    0.017</td> <td>   -1.325</td> <td> 0.186</td> <td>   -0.056</td> <td>    0.011</td>
    </tr>
    <tr>
      <th>DIS</th>     <td>   -1.1695</td> <td>    0.249</td> <td>   -4.692</td> <td> 0.000</td> <td>   -1.660</td> <td>   -0.679</td>
    </tr>
    <tr>
      <th>RAD</th>     <td>    0.2218</td> <td>    0.083</td> <td>    2.672</td> <td> 0.008</td> <td>    0.058</td> <td>    0.385</td>
    </tr>
    <tr>
      <th>TAX</th>     <td>   -0.0104</td> <td>    0.005</td> <td>   -2.155</td> <td> 0.032</td> <td>   -0.020</td> <td>   -0.001</td>
    </tr>
    <tr>
      <th>PTRATIO</th> <td>   -0.4480</td> <td>    0.137</td> <td>   -3.261</td> <td> 0.001</td> <td>   -0.718</td> <td>   -0.178</td>
    </tr>
    <tr>
      <th>B</th>       <td>    0.0165</td> <td>    0.003</td> <td>    4.767</td> <td> 0.000</td> <td>    0.010</td> <td>    0.023</td>
    </tr>
    <tr>
      <th>LSTAT</th>   <td>   -0.2592</td> <td>    0.065</td> <td>   -3.959</td> <td> 0.000</td> <td>   -0.388</td> <td>   -0.130</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td>170.276</td> <th>  Durbin-Watson:     </th> <td>   1.981</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1726.398</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td> 1.839</td>  <th>  Prob(JB):          </th> <td>    0.00</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td>13.426</td>  <th>  Cond. No.          </th> <td>8.91e+03</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 8.91e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems.



.. code:: 

    model = sm.OLS(y,X.drop("LSTAT", axis=1))
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>          <td>MEDV</td>       <th>  R-squared:         </th> <td>   0.954</td> 
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.953</td> 
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   846.6</td> 
    </tr>
    <tr>
      <th>Date:</th>             <td>Fri, 07 Jun 2019</td> <th>  Prob (F-statistic):</th> <td>2.38e-320</td>
    </tr>
    <tr>
      <th>Time:</th>                 <td>16:38:39</td>     <th>  Log-Likelihood:    </th> <td> -1556.1</td> 
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>   506</td>      <th>  AIC:               </th> <td>   3136.</td> 
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   494</td>      <th>  BIC:               </th> <td>   3187.</td> 
    </tr>
    <tr>
      <th>Df Model:</th>              <td>    12</td>      <th>                     </th>     <td> </td>    
    </tr>
    <tr>
      <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    
    </tr>
    </table>
    <table class="simpletable">
    <tr>
         <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
    </tr>
    <tr>
      <th>CRIM</th>    <td>   -0.1439</td> <td>    0.036</td> <td>   -3.990</td> <td> 0.000</td> <td>   -0.215</td> <td>   -0.073</td>
    </tr>
    <tr>
      <th>ZN</th>      <td>    0.0413</td> <td>    0.015</td> <td>    2.696</td> <td> 0.007</td> <td>    0.011</td> <td>    0.071</td>
    </tr>
    <tr>
      <th>INDUS</th>   <td>   -0.0370</td> <td>    0.068</td> <td>   -0.540</td> <td> 0.589</td> <td>   -0.172</td> <td>    0.098</td>
    </tr>
    <tr>
      <th>CHAS</th>    <td>    3.2525</td> <td>    0.961</td> <td>    3.384</td> <td> 0.001</td> <td>    1.364</td> <td>    5.141</td>
    </tr>
    <tr>
      <th>NOX</th>     <td>  -10.8653</td> <td>    3.422</td> <td>   -3.175</td> <td> 0.002</td> <td>  -17.590</td> <td>   -4.141</td>
    </tr>
    <tr>
      <th>RM</th>      <td>    7.1436</td> <td>    0.289</td> <td>   24.734</td> <td> 0.000</td> <td>    6.576</td> <td>    7.711</td>
    </tr>
    <tr>
      <th>AGE</th>     <td>   -0.0449</td> <td>    0.014</td> <td>   -3.235</td> <td> 0.001</td> <td>   -0.072</td> <td>   -0.018</td>
    </tr>
    <tr>
      <th>DIS</th>     <td>   -1.2292</td> <td>    0.206</td> <td>   -5.980</td> <td> 0.000</td> <td>   -1.633</td> <td>   -0.825</td>
    </tr>
    <tr>
      <th>RAD</th>     <td>    0.2008</td> <td>    0.071</td> <td>    2.829</td> <td> 0.005</td> <td>    0.061</td> <td>    0.340</td>
    </tr>
    <tr>
      <th>TAX</th>     <td>   -0.0100</td> <td>    0.004</td> <td>   -2.391</td> <td> 0.017</td> <td>   -0.018</td> <td>   -0.002</td>
    </tr>
    <tr>
      <th>PTRATIO</th> <td>   -0.6575</td> <td>    0.112</td> <td>   -5.881</td> <td> 0.000</td> <td>   -0.877</td> <td>   -0.438</td>
    </tr>
    <tr>
      <th>B</th>       <td>    0.0165</td> <td>    0.003</td> <td>    5.779</td> <td> 0.000</td> <td>    0.011</td> <td>    0.022</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td>277.013</td> <th>  Durbin-Watson:     </th> <td>   0.927</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>3084.310</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td> 2.148</td>  <th>  Prob(JB):          </th> <td>    0.00</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td>14.307</td>  <th>  Cond. No.          </th> <td>8.13e+03</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 8.13e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems.



TPOT
~~~~

`TPOT <https://github.com/rhiever/tpot>`__ est un module d’apprentissage
automatique.

.. code:: 

    try:
        from tpot import TPOTRegressor
    except ImportError:
        # for sklearn 0.22
        import sklearn.preprocessing    
        from sklearn.impute import SimpleImputer
        sklearn.preprocessing.Imputer = SimpleImputer
        from tpot import TPOTRegressor
    
    tpot = TPOTRegressor(generations=2, population_size=50, verbosity=2)
    tpot.fit(X_train, y_train)
    print(tpot.score(X_test, y_test))
    tpot.export('tpot_boston_pipeline.py')



.. parsed-literal::
    HBox(children=(IntProgress(value=0, description='Optimization Progress', max=150, style=ProgressStyle(descript…


.. parsed-literal::

    Generation 1 - Current best internal CV score: -10.19955600172555
    Generation 2 - Current best internal CV score: -10.19955600172555
    
    Best pipeline: XGBRegressor(StandardScaler(input_matrix), learning_rate=0.1, max_depth=9, min_child_weight=2, n_estimators=100, nthread=1, subsample=0.8)
    -9.01100125963658


Le module optimise les hyperparamètres, parfois un peu trop à en juger
la mauvaise performance obtenue sur la base de test.

.. code:: 

    r2_score(y_train, tpot.predict(X_train))




.. parsed-literal::
    0.9994364901579464



.. code:: 

    r2_score(y_test, tpot.predict(X_test))




.. parsed-literal::
    0.8978805496355459



Feature importance pour une observations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On reprend la première random forest et on utilise le module
`treeinterpreter <https://github.com/andosa/treeinterpreter>`__.

.. code:: 

    clr = RandomForestRegressor()
    clr.fit(X, y)




.. parsed-literal::
    RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
                          max_features='auto', max_leaf_nodes=None,
                          min_impurity_decrease=0.0, min_impurity_split=None,
                          min_samples_leaf=1, min_samples_split=2,
                          min_weight_fraction_leaf=0.0, n_estimators=100,
                          n_jobs=None, oob_score=False, random_state=None,
                          verbose=0, warm_start=False)



.. code:: 

    from treeinterpreter import treeinterpreter as ti
    prediction, bias, contributions = ti.predict(clr, X_test)

.. code:: 

    for i in range(min(2, X_train.shape[0])):
        print("Instance", i)
        print("Bias (trainset mean)", bias[i])
        print("Feature contributions:")
        for c, feature in sorted(zip(contributions[i], df.columns), 
                                 key=lambda x: -abs(x[0])):
            print(feature, round(c, 2))
        print( "-"*20)


.. parsed-literal::
    Instance 0
    Bias (trainset mean) 22.527195652173905
    Feature contributions:
    RM 7.53
    LSTAT 1.58
    TAX -0.64
    PTRATIO 0.6
    DIS -0.4
    INDUS 0.19
    RAD -0.18
    AGE 0.13
    ZN 0.12
    CRIM 0.07
    NOX 0.06
    B -0.01
    CHAS -0.01
    --------------------
    Instance 1
    Bias (trainset mean) 22.527195652173905
    Feature contributions:
    LSTAT -5.04
    RM -1.64
    B -0.9
    NOX -0.89
    CRIM -0.81
    DIS 0.64
    AGE 0.17
    PTRATIO -0.07
    TAX -0.04
    CHAS -0.02
    INDUS -0.0
    ZN 0.0
    RAD 0.0
    --------------------


Problème 2 : série temporelle
-----------------------------

On prend une série sur `Google
Trends <https://trends.google.fr/trends/explore?q=tennis%20live>`__,
dans notre cas, c’est la requête *tennis live*. On compare une approche
linéaire et une approche non linéaire.

Approche linéaire
~~~~~~~~~~~~~~~~~

.. code:: 

    import pandas
    df = pandas.read_csv("data/multiTimeline.csv", skiprows=1)

.. code:: 

    df.columns= ["Semaine", "compte"]

.. code:: 

    df["SemaineDt"] = pandas.to_datetime(df.Semaine)

.. code:: 

    df=df.set_index("SemaineDt")

.. code:: 

    df["compte"] = df["compte"].astype(float)

.. code:: 

    df.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Semaine</th>
          <th>compte</th>
        </tr>
        <tr>
          <th>SemaineDt</th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2012-07-01</th>
          <td>2012-07-01</td>
          <td>70.0</td>
        </tr>
        <tr>
          <th>2012-07-08</th>
          <td>2012-07-08</td>
          <td>49.0</td>
        </tr>
        <tr>
          <th>2012-07-15</th>
          <td>2012-07-15</td>
          <td>18.0</td>
        </tr>
        <tr>
          <th>2012-07-22</th>
          <td>2012-07-22</td>
          <td>22.0</td>
        </tr>
        <tr>
          <th>2012-07-29</th>
          <td>2012-07-29</td>
          <td>88.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    %matplotlib inline

.. code:: 

    df.plot()




.. parsed-literal::
    <matplotlib.axes._subplots.AxesSubplot at 0x19c15163358>




.. image:: 2017_session6_40_1.png


.. code:: 

    from statsmodels.tsa.arima_model import ARIMA
    arma_mod = ARIMA(df["compte"].values, order=(6 ,1, 1))
    res = arma_mod.fit()
    res.params


.. parsed-literal::

    




.. parsed-literal::

    array([ 0.00418581,  0.59035757, -0.32540695,  0.23286807, -0.03300838,
            0.06434307, -0.07204017, -0.99999983])



.. code:: 

    res.summary()






.. raw:: html

    <table class="simpletable">
    <caption>ARIMA Model Results</caption>
    <tr>
      <th>Dep. Variable:</th>        <td>D.y</td>       <th>  No. Observations:  </th>    <td>259</td>   
    </tr>
    <tr>
      <th>Model:</th>          <td>ARIMA(6, 1, 1)</td>  <th>  Log Likelihood     </th> <td>-1055.581</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>css-mle</td>     <th>  S.D. of innovations</th>  <td>14.116</td>  
    </tr>
    <tr>
      <th>Date:</th>          <td>Fri, 07 Jun 2019</td> <th>  AIC                </th> <td>2129.161</td> 
    </tr>
    <tr>
      <th>Time:</th>              <td>16:45:54</td>     <th>  BIC                </th> <td>2161.173</td> 
    </tr>
    <tr>
      <th>Sample:</th>                <td>1</td>        <th>  HQIC               </th> <td>2142.032</td> 
    </tr>
    <tr>
      <th></th>                       <td> </td>        <th>                     </th>     <td> </td>    
    </tr>
    </table>
    <table class="simpletable">
    <tr>
          <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
    </tr>
    <tr>
      <th>const</th>     <td>    0.0042</td> <td>    0.021</td> <td>    0.196</td> <td> 0.845</td> <td>   -0.038</td> <td>    0.046</td>
    </tr>
    <tr>
      <th>ar.L1.D.y</th> <td>    0.5904</td> <td>    0.063</td> <td>    9.431</td> <td> 0.000</td> <td>    0.468</td> <td>    0.713</td>
    </tr>
    <tr>
      <th>ar.L2.D.y</th> <td>   -0.3254</td> <td>    0.072</td> <td>   -4.507</td> <td> 0.000</td> <td>   -0.467</td> <td>   -0.184</td>
    </tr>
    <tr>
      <th>ar.L3.D.y</th> <td>    0.2329</td> <td>    0.075</td> <td>    3.097</td> <td> 0.002</td> <td>    0.085</td> <td>    0.380</td>
    </tr>
    <tr>
      <th>ar.L4.D.y</th> <td>   -0.0330</td> <td>    0.076</td> <td>   -0.433</td> <td> 0.665</td> <td>   -0.182</td> <td>    0.116</td>
    </tr>
    <tr>
      <th>ar.L5.D.y</th> <td>    0.0643</td> <td>    0.076</td> <td>    0.842</td> <td> 0.400</td> <td>   -0.085</td> <td>    0.214</td>
    </tr>
    <tr>
      <th>ar.L6.D.y</th> <td>   -0.0720</td> <td>    0.066</td> <td>   -1.096</td> <td> 0.274</td> <td>   -0.201</td> <td>    0.057</td>
    </tr>
    <tr>
      <th>ma.L1.D.y</th> <td>   -1.0000</td> <td>    0.010</td> <td>  -96.075</td> <td> 0.000</td> <td>   -1.020</td> <td>   -0.980</td>
    </tr>
    </table>
    <table class="simpletable">
    <caption>Roots</caption>
    <tr>
        <td></td>   <th>            Real</th>  <th>         Imaginary</th> <th>         Modulus</th>  <th>        Frequency</th>
    </tr>
    <tr>
      <th>AR.1</th> <td>          -1.2011</td> <td>          -1.2144j</td> <td>           1.7080</td> <td>          -0.3741</td>
    </tr>
    <tr>
      <th>AR.2</th> <td>          -1.2011</td> <td>          +1.2144j</td> <td>           1.7080</td> <td>           0.3741</td>
    </tr>
    <tr>
      <th>AR.3</th> <td>           0.1840</td> <td>          -1.4018j</td> <td>           1.4138</td> <td>          -0.2292</td>
    </tr>
    <tr>
      <th>AR.4</th> <td>           0.1840</td> <td>          +1.4018j</td> <td>           1.4138</td> <td>           0.2292</td>
    </tr>
    <tr>
      <th>AR.5</th> <td>           1.4636</td> <td>          -0.4882j</td> <td>           1.5429</td> <td>          -0.0512</td>
    </tr>
    <tr>
      <th>AR.6</th> <td>           1.4636</td> <td>          +0.4882j</td> <td>           1.5429</td> <td>           0.0512</td>
    </tr>
    <tr>
      <th>MA.1</th> <td>           1.0000</td> <td>          +0.0000j</td> <td>           1.0000</td> <td>           0.0000</td>
    </tr>
    </table>



Méthode non linéaire
~~~~~~~~~~~~~~~~~~~~

On construire la matrice des séries décalées. Cette méthode permet de
sortir du cadre linéaire et d’ajouter d’autres variables.

.. code:: 

    from statsmodels.tsa.tsatools import lagmat
    lag = 8
    X = lagmat(df["compte"], lag)
    lagged = df.copy()
    for c in range(1,lag+1):
        lagged["lag%d" % c] = X[:, c-1]
    pandas.concat([lagged.head(), lagged.tail()])






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Semaine</th>
          <th>compte</th>
          <th>lag1</th>
          <th>lag2</th>
          <th>lag3</th>
          <th>lag4</th>
          <th>lag5</th>
          <th>lag6</th>
          <th>lag7</th>
          <th>lag8</th>
        </tr>
        <tr>
          <th>SemaineDt</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2012-07-01</th>
          <td>2012-07-01</td>
          <td>70.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2012-07-08</th>
          <td>2012-07-08</td>
          <td>49.0</td>
          <td>70.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2012-07-15</th>
          <td>2012-07-15</td>
          <td>18.0</td>
          <td>49.0</td>
          <td>70.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2012-07-22</th>
          <td>2012-07-22</td>
          <td>22.0</td>
          <td>18.0</td>
          <td>49.0</td>
          <td>70.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2012-07-29</th>
          <td>2012-07-29</td>
          <td>88.0</td>
          <td>22.0</td>
          <td>18.0</td>
          <td>49.0</td>
          <td>70.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
          <td>0.0</td>
        </tr>
        <tr>
          <th>2017-05-21</th>
          <td>2017-05-21</td>
          <td>23.0</td>
          <td>40.0</td>
          <td>35.0</td>
          <td>21.0</td>
          <td>29.0</td>
          <td>27.0</td>
          <td>14.0</td>
          <td>23.0</td>
          <td>40.0</td>
        </tr>
        <tr>
          <th>2017-05-28</th>
          <td>2017-05-28</td>
          <td>44.0</td>
          <td>23.0</td>
          <td>40.0</td>
          <td>35.0</td>
          <td>21.0</td>
          <td>29.0</td>
          <td>27.0</td>
          <td>14.0</td>
          <td>23.0</td>
        </tr>
        <tr>
          <th>2017-06-04</th>
          <td>2017-06-04</td>
          <td>55.0</td>
          <td>44.0</td>
          <td>23.0</td>
          <td>40.0</td>
          <td>35.0</td>
          <td>21.0</td>
          <td>29.0</td>
          <td>27.0</td>
          <td>14.0</td>
        </tr>
        <tr>
          <th>2017-06-11</th>
          <td>2017-06-11</td>
          <td>28.0</td>
          <td>55.0</td>
          <td>44.0</td>
          <td>23.0</td>
          <td>40.0</td>
          <td>35.0</td>
          <td>21.0</td>
          <td>29.0</td>
          <td>27.0</td>
        </tr>
        <tr>
          <th>2017-06-18</th>
          <td>2017-06-18</td>
          <td>28.0</td>
          <td>28.0</td>
          <td>55.0</td>
          <td>44.0</td>
          <td>23.0</td>
          <td>40.0</td>
          <td>35.0</td>
          <td>21.0</td>
          <td>29.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    xc = ["lag%d" % i for i in range(1,lag+1)]
    split = 0.66
    isplit = int(len(lagged) * split)
    xt = lagged[10:][xc]
    yt = lagged[10:]["compte"]
    X_train, y_train, X_test, y_test = xt[:isplit], yt[:isplit], xt[isplit:], yt[isplit:]

.. code:: 

    from sklearn.ensemble import RandomForestRegressor
    clr = RandomForestRegressor()
    clr.fit(X_train, y_train)




.. parsed-literal::
    RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
                          max_features='auto', max_leaf_nodes=None,
                          min_impurity_decrease=0.0, min_impurity_split=None,
                          min_samples_leaf=1, min_samples_split=2,
                          min_weight_fraction_leaf=0.0, n_estimators=100,
                          n_jobs=None, oob_score=False, random_state=None,
                          verbose=0, warm_start=False)



.. code:: 

    from sklearn.metrics import r2_score
    r2 = r2_score(y_test.values, clr.predict(X_test))
    r2


.. parsed-literal::

    




.. parsed-literal::

    0.4878325848907805



.. code:: 

    plt.scatter(y_test.values, clr.predict(X_test));


.. parsed-literal::

    




.. parsed-literal::

    <matplotlib.collections.PathCollection at 0x19c26b9a2e8>




.. image:: 2017_session6_48_2.png


Texte
-----

On cherche à comparer une
`LDA <http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.LatentDirichletAllocation.html>`__
avec `word2vec <https://radimrehurek.com/gensim/models/word2vec.html>`__
et
`kmeans <http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html>`__
et les données qui sont sur
`ensae_teaching_cs/src/ensae_teaching_cs/data/data_web/ <https://github.com/sdpython/ensae_teaching_cs/tree/master/src/ensae_teaching_cs/data/data_web>`__.

.. code:: 

    from ensae_teaching_cs.data import twitter_zip
    df = twitter_zip(as_df=True)
    df.head(n=2).T






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>0</th>
          <th>1</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>index</th>
          <td>776066992054861825</td>
          <td>776067660979245056</td>
        </tr>
        <tr>
          <th>nb_user_mentions</th>
          <td>0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>nb_extended_entities</th>
          <td>0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>nb_hashtags</th>
          <td>1</td>
          <td>1</td>
        </tr>
        <tr>
          <th>geo</th>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>text_hashtags</th>
          <td>, SiJétaisPrésident</td>
          <td>, SiJétaisPrésident</td>
        </tr>
        <tr>
          <th>annee</th>
          <td>2016</td>
          <td>2016</td>
        </tr>
        <tr>
          <th>delimit_mention</th>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>lang</th>
          <td>fr</td>
          <td>fr</td>
        </tr>
        <tr>
          <th>id_str</th>
          <td>7.76067e+17</td>
          <td>7.76068e+17</td>
        </tr>
        <tr>
          <th>text_mention</th>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>retweet_count</th>
          <td>4</td>
          <td>5</td>
        </tr>
        <tr>
          <th>favorite_count</th>
          <td>3</td>
          <td>8</td>
        </tr>
        <tr>
          <th>type_extended_entities</th>
          <td>[]</td>
          <td>[]</td>
        </tr>
        <tr>
          <th>text</th>
          <td>#SiJétaisPrésident se serait la fin du monde.....</td>
          <td>#SiJétaisPrésident je donnerai plus de vacance...</td>
        </tr>
        <tr>
          <th>nb_user_photos</th>
          <td>0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>nb_urls</th>
          <td>0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>nb_symbols</th>
          <td>0</td>
          <td>0</td>
        </tr>
        <tr>
          <th>created_at</th>
          <td>Wed Sep 14 14:36:04 +0000 2016</td>
          <td>Wed Sep 14 14:38:43 +0000 2016</td>
        </tr>
        <tr>
          <th>delimit_hash</th>
          <td>, 0, 18</td>
          <td>, 0, 18</td>
        </tr>
      </tbody>
    </table>
    </div>



Des mots aux coordonnées - tf-idf
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: 

    keep = df.text.dropna().index

.. code:: 

    dfnonan = df.iloc[keep, :]
    dfnonan.shape




.. parsed-literal::
    (5087, 20)



.. code:: 

    from sklearn.feature_extraction.text import TfidfVectorizer
    tfidf_vectorizer = TfidfVectorizer(max_df=0.95, min_df=2, max_features=1000)
    tfidf = tfidf_vectorizer.fit_transform(dfnonan["text"])

.. code:: 

    tfidf[:2, :]




.. parsed-literal::
    <2x1000 sparse matrix of type '<class 'numpy.float64'>'
    	with 21 stored elements in Compressed Sparse Row format>



.. code:: 

    tfidf[:2, :].todense()




.. parsed-literal::
    matrix([[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]])



LDA
~~~

.. code:: 

    from sklearn.decomposition import LatentDirichletAllocation
    lda = LatentDirichletAllocation(n_components=10, max_iter=5,
                                    learning_method='online',
                                    learning_offset=50.,
                                    random_state=0)

.. code:: 

    lda.fit(tfidf)




.. parsed-literal::
    LatentDirichletAllocation(batch_size=128, doc_topic_prior=None,
                              evaluate_every=-1, learning_decay=0.7,
                              learning_method='online', learning_offset=50.0,
                              max_doc_update_iter=100, max_iter=5,
                              mean_change_tol=0.001, n_components=10, n_jobs=None,
                              perp_tol=0.1, random_state=0, topic_word_prior=None,
                              total_samples=1000000.0, verbose=0)



.. code:: 

    tf_feature_names = tfidf_vectorizer.get_feature_names()
    tf_feature_names[100:103]




.. parsed-literal::
    ['avoir', 'bac', 'bah']



.. code:: 

    lda.components_.shape




.. parsed-literal::
    (10, 1000)



On obtient dix vecteurs qui représentent les dix vecteurs associés aux
dix clusters. Chaque dimension relié au fait que le mot appartient ou
non au cluster.

.. code:: 

    def print_top_words(model, feature_names, n_top_words):
        for topic_idx, topic in enumerate(model.components_):
            print("Topic #%d:" % topic_idx)
            print(" ".join([feature_names[i]
                            for i in topic.argsort()[:-n_top_words - 1:-1]]))
        print()

.. code:: 

    print_top_words(lda, tf_feature_names, 10)


.. parsed-literal::
    Topic #0:
    gratuit mcdo supprimerai école soir kebab macdo kfc domicile cc
    Topic #1:
    macron co https de la est le il et hollande
    Topic #2:
    sijetaispresident je les de la et le des en pour
    Topic #3:
    notaires eu organiserais mets carte nouveaux journées installation cache créer
    Topic #4:
    sijetaispresident interdirais les je ballerines la serait serais bah de
    Topic #5:
    ministre de sijetaispresident la je premier mort et nommerais président
    Topic #6:
    cours le supprimerais jour sijetaispresident lundi samedi semaine je vendredi
    Topic #7:
    port interdirait démissionnerais promesses heure rendrai ballerine mes changement christineboutin
    Topic #8:
    seraient sijetaispresident gratuits aux les nos putain éducation nationale bonne
    Topic #9:
    bordel seront légaliserai putes gratuites pizza mot virerais vitesse dutreil
    


Clustering
~~~~~~~~~~

.. code:: 

    from sklearn.cluster import KMeans
    km = KMeans(n_clusters=10)
    km.fit(tfidf)




.. parsed-literal::
    KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
           n_clusters=10, n_init=10, n_jobs=None, precompute_distances='auto',
           random_state=None, tol=0.0001, verbose=0)



.. code:: 

    km.cluster_centers_.shape




.. parsed-literal::
    (10, 1000)



.. code:: 

    def print_top_words(model, feature_names, n_top_words):
        for topic_idx, topic in enumerate(model.cluster_centers_):
            print("Topic #%d:" % topic_idx)
            print(" ".join([feature_names[i]
                            for i in topic.argsort()[:-n_top_words - 1:-1]]))
        print()

.. code:: 

    print_top_words(km, tf_feature_names, 10)


.. parsed-literal::
    Topic #0:
    ferais je sijetaispresident les en sorte que de et un
    Topic #1:
    le sijetaispresident je monde et pour de tout des les
    Topic #2:
    de sijetaispresident la en et des une un plus aurait
    Topic #3:
    https co macron de la le les via et sijetaispresident
    Topic #4:
    les sijetaispresident je de et tous des seraient pour ballerines
    Topic #5:
    serait sijetaispresident la le merde de ça et on dans
    Topic #6:
    je sijetaispresident la des me et ferai pas en au
    Topic #7:
    macron est il de pas la hollande le gauche que
    Topic #8:
    ministre premier nommerais sijetaispresident de je la mickey serait nommerai
    Topic #9:
    serais je sijetaispresident un pas ne président que la de
    


Ah les stop words….