.. _artificieldurationpredictionrst:

======================
Prédiction d’une durée
======================


.. only:: html

    **Links:** :download:`notebook <artificiel_duration_prediction.ipynb>`, :downloadlink:`html <artificiel_duration_prediction2html.html>`, :download:`PDF <artificiel_duration_prediction.pdf>`, :download:`python <artificiel_duration_prediction.py>`, :downloadlink:`slides <artificiel_duration_prediction.slides.html>`, :githublink:`GitHub|_doc/notebooks/lectures/artificiel_duration_prediction.ipynb|*`


Ce notebook cherche à prédire la durée de stockage de paquets préparés
par un magasin. Chaque paquet met plus ou moins de temps à être préparé.
Si la commande arrive le soir et ne peut être finie avant la fin de la
journée, elle est reportée sur la journée suivante. C’est la
particularité de ce jeu de données.


-  L’heure de la commande est uniformément répartie sur la semaine et la
   journée.
-  La durée suit une loi :math:`\Gamma(3, 0.5)`.

.. code:: 

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





Distribution gamma
------------------

La loi `Gamma <https://fr.wikipedia.org/wiki/Loi_Gamma>`__ est souvent
utilisée pour modéliser une durée de vie.

.. code:: 

    %matplotlib inline

.. code:: 

    from numpy.random import gamma
    g = gamma(6, 0.25, 1000)
    
    import matplotlib.pyplot as plt
    import pandas
    fig, ax = plt.subplots(1, 1, figsize=(8, 3))
    df = pandas.DataFrame(dict(g=g))
    df.hist('g', bins=50, ax=ax)
    ax.set_title('Distribution gamma(6, 0.25)')
    ax.set_xlim([0, 4]);



.. image:: artificiel_duration_prediction_5_0.png


On simule un jeu de données avec la fonction
`duration_selling <http://www.xavierdupre.fr/app/papierstat/helpsphinx/papierstat/datasets/duration.html?papierstat.datasets.duration.duration_selling>`__.

.. code:: 

    from papierstat.datasets import duration_selling
    df = duration_selling()
    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>commande</th>
          <th>reception</th>
          <th>true_duration</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>2018-02-05 13:15:09.062249</td>
          <td>2018-02-05 14:18:58.540886</td>
          <td>1.063744</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2018-02-05 17:16:52.648499</td>
          <td>2018-02-05 18:21:28.592673</td>
          <td>1.076651</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2018-02-05 09:37:59.627662</td>
          <td>2018-02-05 10:58:42.387613</td>
          <td>1.345211</td>
        </tr>
        <tr>
          <th>3</th>
          <td>2018-02-05 09:02:20.357816</td>
          <td>2018-02-05 10:03:05.420104</td>
          <td>1.012517</td>
        </tr>
        <tr>
          <th>4</th>
          <td>2018-02-05 13:52:33.335866</td>
          <td>2018-02-05 15:20:27.454535</td>
          <td>1.465033</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    df.shape




.. parsed-literal::
    (2481, 3)



Quelques statistiques
---------------------

.. code:: 

    stat = df.copy()
    stat["wk"] = stat.commande.dt.weekday
    stat["wk2"] = stat.reception.dt.weekday
    stat["hr"] = stat.commande.dt.hour
    stat["hr2"] = stat.reception.dt.hour

.. code:: 

    wk = stat.groupby('wk').count()[["commande"]]
    wk2 = stat.groupby('wk2').count()[["reception"]]
    wks = wk.merge(wk2, left_on=wk.index, right_on=wk2.index)
    wks.columns = ['jour'] + list(wks.columns)[1:]
    wks






.. 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>jour</th>
          <th>commande</th>
          <th>reception</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0</td>
          <td>502</td>
          <td>498</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1</td>
          <td>531</td>
          <td>529</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2</td>
          <td>535</td>
          <td>544</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3</td>
          <td>443</td>
          <td>458</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4</td>
          <td>470</td>
          <td>452</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1, 1, figsize=(8, 3))
    wks.plot(x="jour", y=["commande", "reception"], kind="bar", ax=ax)
    ax.set_xlabel('jour de la semaine')
    ax.set_title("Distribution des commandes sur la semaine")
    ax.set_ylim([440,560]);



.. image:: artificiel_duration_prediction_12_0.png


.. code:: 

    hr = stat.groupby('hr').count()[["commande"]]
    hr2 = stat.groupby('hr2').count()[["reception"]]
    hrs = hr.merge(hr2, left_on=hr.index, right_on=hr2.index)
    hrs.columns = ['heure'] + list(wks.columns)[1:]
    hrs






.. 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>heure</th>
          <th>commande</th>
          <th>reception</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>9</td>
          <td>260</td>
          <td>249</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10</td>
          <td>233</td>
          <td>262</td>
        </tr>
        <tr>
          <th>2</th>
          <td>11</td>
          <td>244</td>
          <td>264</td>
        </tr>
        <tr>
          <th>3</th>
          <td>12</td>
          <td>239</td>
          <td>223</td>
        </tr>
        <tr>
          <th>4</th>
          <td>13</td>
          <td>268</td>
          <td>237</td>
        </tr>
        <tr>
          <th>5</th>
          <td>14</td>
          <td>249</td>
          <td>266</td>
        </tr>
        <tr>
          <th>6</th>
          <td>15</td>
          <td>249</td>
          <td>241</td>
        </tr>
        <tr>
          <th>7</th>
          <td>16</td>
          <td>247</td>
          <td>256</td>
        </tr>
        <tr>
          <th>8</th>
          <td>17</td>
          <td>228</td>
          <td>262</td>
        </tr>
        <tr>
          <th>9</th>
          <td>18</td>
          <td>264</td>
          <td>221</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    fig, ax = plt.subplots(1, 1, figsize=(8, 3))
    hrs.plot(x="heure", y=["commande", "reception"], kind="bar", ax=ax)
    ax.set_xlabel('heure du jour')
    ax.set_title("Heures des commandes sur la semaine")
    ax.set_ylim([200,300]);



.. image:: artificiel_duration_prediction_14_0.png


.. code:: 

    fig, ax = plt.subplots(1, 3, figsize=(12, 3))
    stat["duree"] = stat.reception - stat.commande
    stat["hours"] = stat.duree.dt.total_seconds() / 3600
    stat['hours'].hist(bins=500, ax=ax[0])
    stat[stat['hours'] < 10]['hours'].hist(bins=50, ax=ax[1])
    stat['true_duration'].hist(bins=50, ax=ax[2])
    ax[0].set_title("Distribution de la durée")
    ax[1].set_title("Distribution de la durée\ntronquée à 10h")
    ax[2].set_title("Distribution de la durée\nayant servi à l'estimation")
    ax[0].set_xlabel("Heures")
    ax[1].set_xlabel("Heures")
    ax[2].set_xlabel("Heures");



.. image:: artificiel_duration_prediction_15_0.png


La distribution montre trois pics : le premier correspond aux paquets
réceptionnés dans la journée, le second le lendemain et le troisième le
week-end. La durée de livraison dépend uniquement de l’heure d’arrivée
de la commande. Si elle arrive assez tôt, la livraison est faite
aléatoirement dans la journée, sinon, on passe au jour suivant.

.. code:: 

    import seaborn
    fig, ax = plt.subplots(1, 2, figsize=(12, 4))
    seaborn.boxplot(x="hr", y="hours", data=stat, palette="muted", ax=ax[0])
    seaborn.boxplot(x="hr", y="hours", data=stat, palette="muted", ax=ax[1])
    ax[1].set_ylim([0, 15])
    ax[0].set_title("Distribution de la durée de livraison\nen fonction de l'heure d'arrivée")
    ax[1].set_title("Distribution tronquée de la durée de livraison\nen fonction de l'heure d'arrivée");



.. image:: artificiel_duration_prediction_17_0.png


Prédiction de la durée
----------------------

Il n’y a pas grand-chose à prédire puisque tout est aléatoire excepté le
fait que le fait qu’un paquet soit réceptionné le lendemain ou le lundi
suivant. On construit donc une variable correspondant à l’heure de
commande et le jour de la semaine. On compare un modèle linéaire et un
arbre de décision.

Comparaison d’un modèle linéaire et d’un arbre de précision
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: 

    data = df.copy()
    data["heure"] = data.commande.dt.hour + data.commande.dt.minute / 60
    data["wk"] = data.commande.dt.weekday
    data['duree'] = (data.reception - data.commande).dt.total_seconds() / 3600
    data.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>commande</th>
          <th>reception</th>
          <th>true_duration</th>
          <th>heure</th>
          <th>wk</th>
          <th>duree</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>2018-02-05 13:15:09.062249</td>
          <td>2018-02-05 14:18:58.540886</td>
          <td>1.063744</td>
          <td>13.250000</td>
          <td>0</td>
          <td>1.063744</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2018-02-05 17:16:52.648499</td>
          <td>2018-02-05 18:21:28.592673</td>
          <td>1.076651</td>
          <td>17.266667</td>
          <td>0</td>
          <td>1.076651</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2018-02-05 09:37:59.627662</td>
          <td>2018-02-05 10:58:42.387613</td>
          <td>1.345211</td>
          <td>9.616667</td>
          <td>0</td>
          <td>1.345211</td>
        </tr>
        <tr>
          <th>3</th>
          <td>2018-02-05 09:02:20.357816</td>
          <td>2018-02-05 10:03:05.420104</td>
          <td>1.012517</td>
          <td>9.033333</td>
          <td>0</td>
          <td>1.012517</td>
        </tr>
        <tr>
          <th>4</th>
          <td>2018-02-05 13:52:33.335866</td>
          <td>2018-02-05 15:20:27.454535</td>
          <td>1.465033</td>
          <td>13.866667</td>
          <td>0</td>
          <td>1.465033</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: 

    X = data[["heure", 'wk']]
    y = data["duree"]

.. code:: 

    from sklearn.model_selection import train_test_split
    X_train, X_test, y_train, y_test = train_test_split(X, y)

.. code:: 

    from sklearn.linear_model import LinearRegression
    clr = LinearRegression()
    clr.fit(X_train, y_train)
    clr.score(X_test, y_test)




.. parsed-literal::
    0.20497865247039226



.. code:: 

    from sklearn.tree import DecisionTreeRegressor
    clr2 = DecisionTreeRegressor(max_depth=3)
    clr2.fit(X_train, y_train)
    clr2.score(X_test, y_test)




.. parsed-literal::
    0.6591302839171569



.. code:: 

    fig, ax = plt.subplots(1, 1, figsize=(12,4))
    ax.plot(y_test, clr.predict(X_test), '.', label="linéaire")
    ax.plot(y_test, clr2.predict(X_test), '.', label="arbre")
    ax.plot([0,70], [0,70], '--', label="prédiction parfaite")
    ax.set_xlabel("valeur attendue")
    ax.set_ylabel("valeur prédite")
    ax.legend();



.. image:: artificiel_duration_prediction_25_0.png


L’arbre de décision est nettement meilleur simplement parce que la
relation entre la valeur à prédire et les variables n’est pas linéaire.

.. code:: 

    from seaborn import pairplot
    xy = X.copy()
    xy['duree'] = y
    pairplot(xy);



.. image:: artificiel_duration_prediction_27_0.png


.. code:: 

    from sklearn.tree import export_graphviz
    dot = export_graphviz(clr2, out_file=None, feature_names=X.columns)
    from jyquickhelper import RenderJsDot
    RenderJsDot(dot)






.. raw:: html

    <div id="Md046ff54fc7a47b99ea561b80ea9a284-cont"><div id="Md046ff54fc7a47b99ea561b80ea9a284" style="width:100%;height:100%;"></div></div>
    <script>

    require(['http://www.xavierdupre.fr/js/vizjs/viz.js'], function() { var svgGraph = Viz("digraph Tree {\nnode [shape=box] ;\n0 [label=\"heure <= 17.642\nmse = 145.685\nsamples = 1860\nvalue = 5.38\"] ;\n1 [label=\"heure <= 16.908\nmse = 30.643\nsamples = 1597\nvalue = 2.28\"] ;\n0 -> 1 [labeldistance=2.5, labelangle=45, headlabel=\"True\"] ;\n2 [label=\"heure <= 16.758\nmse = 3.876\nsamples = 1473\nvalue = 1.579\"] ;\n1 -> 2 ;\n3 [label=\"mse = 1.024\nsamples = 1450\nvalue = 1.529\"] ;\n2 -> 3 ;\n4 [label=\"mse = 173.978\nsamples = 23\nvalue = 4.672\"] ;\n2 -> 4 ;\n5 [label=\"wk <= 3.5\nmse = 273.284\nsamples = 124\nvalue = 10.615\"] ;\n1 -> 5 ;\n6 [label=\"mse = 55.104\nsamples = 100\nvalue = 7.227\"] ;\n5 -> 6 ;\n7 [label=\"mse = 935.174\nsamples = 24\nvalue = 24.734\"] ;\n5 -> 7 ;\n8 [label=\"wk <= 3.5\nmse = 431.543\nsamples = 263\nvalue = 24.205\"] ;\n0 -> 8 [labeldistance=2.5, labelangle=-45, headlabel=\"False\"] ;\n9 [label=\"heure <= 17.792\nmse = 20.839\nsamples = 203\nvalue = 14.081\"] ;\n8 -> 9 ;\n10 [label=\"mse = 56.652\nsamples = 21\nvalue = 9.53\"] ;\n9 -> 10 ;\n11 [label=\"mse = 14.041\nsamples = 182\nvalue = 14.606\"] ;\n9 -> 11 ;\n12 [label=\"heure <= 17.725\nmse = 301.289\nsamples = 60\nvalue = 58.455\"] ;\n8 -> 12 ;\n13 [label=\"mse = 976.78\nsamples = 7\nvalue = 28.082\"] ;\n12 -> 13 ;\n14 [label=\"mse = 74.138\nsamples = 53\nvalue = 62.466\"] ;\n12 -> 14 ;\n}");
    document.getElementById('Md046ff54fc7a47b99ea561b80ea9a284').innerHTML = svgGraph; });

    </script>



On voit bien que l’arbre de décision sépare le début et la fin de la
journée. Après 17h30, la durée est augmentée de 8h en moyenne, si c’est
vendredi, la durée est augmentée d’une vingtaine d’heure. Un modèle
linéaire ne peut pas prendre cela en compte.

Remarque sur le score R2
~~~~~~~~~~~~~~~~~~~~~~~~

On voit que le score de prédiction :math:`R^2` est nettement meilleur
pour l’arbre de décision mais c’est un petit peu en trompe l’oeil. La
durée de la période d’inactivité joue un rôle non négligeable. Plus elle
est longue, plus la durée entre la commande et la réception est grande
en moyenne et c’est cette donnée que le modèle linéaire n’arrive pas à
prendre en compte.

.. code:: 

    from sklearn.model_selection import cross_val_score
    from tqdm import tqdm
    import numpy
    res = []
    for end in tqdm(range(1, 25)):
        scores_tree = []
        scores_lin = []
        for i in range(0, 5):
            df = duration_selling(hour_end=end, hour_begin=0, week_pattern=[1] * 7)
            data = df.copy()
            data["heure"] = data.commande.dt.hour + data.commande.dt.minute / 60
            data["wk"] = data.commande.dt.weekday
            data['duree'] = (data.reception - data.commande).dt.total_seconds() / 3600    
            X = data[["heure", 'wk']]
            y = data["duree"]
            clr = LinearRegression()
            clr2 = DecisionTreeRegressor(max_depth=3)
            score_tree = cross_val_score(clr2, X, y, cv=5).mean()
            score_lin = cross_val_score(clr, X, y, cv=5).mean()
            scores_tree.append(score_tree)
            scores_lin.append(score_lin)
        res.append(dict(end=end, score_tree=numpy.array(scores_tree).mean(),
                       score_lin=numpy.array(scores_lin).mean()))
    res = pandas.DataFrame(res)


.. parsed-literal::
    100%|██████████████████████████████████████████████████████████████████████████████████| 24/24 [00:09<00:00,  2.83it/s]


.. code:: 

    ax = res.plot(x='end', y=['score_tree', 'score_lin'], figsize=(5,4))
    ax.set_title("Evolution du R2 en fonction\ndu nombre d'heures d'ouvertures")
    ax.set_xlabel("Heures d'ouverture")
    ax.set_ylabel("R2");



.. image:: artificiel_duration_prediction_32_0.png


On voit que le :math:`R^2` pour l’arbre de décision et le modèle tendent
vers 0 lorsque le nombre d’heures d’ouverture tend vers 24, cela
signifie qu’il n’y a plus d’effets de seuil : la durée entre la
réception et la commande suit une loi gamma et est complètement
aléatoire. C’est attendu. A l’opposé, lorsque le nombre d’heures
d’ouverture est très faible, la durée entre la commande et la réception
est quasiment un mulitple de 24h puisque le magasin n’est quasiment
jamais ouvert : l’arbre de décision est nettement meilleur dans ce cas.
Le pic de :math:`R^2` est obtenu vers 4-5h d’ouvertures, cela correspond
au moment où il y a au pire une seule nuit entre la réception et la
commande.

Pour ce problème, le coefficient :math:`R^2` ne paraît pas être la bonne
solution puisque de toute façon, le modèle prédit simplement si le
paquet sera réceptionné le jour même ou le suivant. Ce serait plutôt un
problème de classification binaire.