.. _regressionlineairerst:

===================
Régression linéaire
===================


.. only:: html

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


Ce notebook s’intéresse à la façon d’interpréter les résultats d’une
régression linéaire lorsque les variables sont corrélées puis il explore
une façon d’associer arbre de décision et régression linéaire pour
construire une régression linéaire par morceaux.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





.. code:: ipython3

    %matplotlib inline

Un cas simple
-------------

Une façon d’interpréter des résultats statistiques est de les calculer
dans un cas où la réponse cherchée est connue. On simule un modèle
simple :math:`Y=\alpha X_1 + 0.X_2 + \epsilon` et on cale une régression
linéaire. On suppose que :math:`X_1, X_2, \epsilon` sont des variables
aléatoires gaussiennes de même variance et moyenne.

.. code:: ipython3

    import numpy.random as npr
    eps = npr.normal(1000)
    X = npr.normal(size=(1000, 3))
    alpha = 2
    Y = alpha * X[:,0] + X[:, 2]
    X.shape, Y.shape




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



.. code:: ipython3

    from numpy import corrcoef
    corrcoef(X.T)




.. parsed-literal::
    array([[ 1.        , -0.0312982 ,  0.05188551],
           [-0.0312982 ,  1.        , -0.00356494],
           [ 0.05188551, -0.00356494,  1.        ]])



.. code:: ipython3

    from statsmodels.regression.linear_model import OLS

.. code:: ipython3

    model = OLS(Y,X[:, :2])
    results = model.fit()
    su = results.summary()
    su






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.815</td>
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.815</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   2204.</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> 
    </tr>
    <tr>
      <th>Time:</th>                 <td>10:34:12</td>     <th>  Log-Likelihood:    </th> <td> -1385.2</td>
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   2774.</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   998</td>      <th>  BIC:               </th> <td>   2784.</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     2</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>x1</th> <td>    2.0519</td> <td>    0.031</td> <td>   66.347</td> <td> 0.000</td> <td>    1.991</td> <td>    2.113</td>
    </tr>
    <tr>
      <th>x2</th> <td>   -0.0032</td> <td>    0.033</td> <td>   -0.097</td> <td> 0.922</td> <td>   -0.067</td> <td>    0.061</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td> 0.709</td> <th>  Durbin-Watson:     </th> <td>   1.990</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.701</td> <th>  Jarque-Bera (JB):  </th> <td>   0.674</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td> 0.063</td> <th>  Prob(JB):          </th> <td>   0.714</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 3.010</td> <th>  Cond. No.          </th> <td>    1.07</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



.. code:: ipython3

    results.rsquared, results.rsquared_adj




.. parsed-literal::
    (0.8153831029946165, 0.8150131292531227)



On vérifie que le coefficient devant :math:`X_1` est non nul (P-value
nulle, 0 n’est pas l’intervalle de confiance). Le coefficient devant
:math:`X_2` n’est pas nul mais presque, la P-value est élevée, le
coefficient :math:`R^2` est élevé. Dessinons.

.. code:: ipython3

    import matplotlib.pyplot as plt
    import seaborn
    fig, ax = plt.subplots(1, 2, figsize=(10,4))
    ax[0].plot(X[:, 0], Y, '.')
    seaborn.kdeplot(X[:, 0], Y, cmap="Reds", shade=True, shade_lowest=False, ax=ax[1])
    ax[0].set_title("nuage de points")
    ax[1].set_title("estimation de la densité");


.. parsed-literal::
    c:\python370_x64\lib\site-packages\scipy\stats\stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
      return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval



.. image:: regression_lineaire_10_1.png


Evolution de R2
---------------

Dans la régression précédente, le coefficient :math:`R^2` transcrit en
quelque sorte la part du bruit :math:`\epsilon` par rapport au terme
:math:`\alpha X_1`. Faisons varier :math:`\alpha`.

.. code:: ipython3

    alphas = []
    r2s = []
    for a in [0.1 * i for i in range(0, 50)]:
        Y = a*X[:,0] + X[:, 2]
        model = OLS(Y,X[:, :2])
        results = model.fit()
        alphas.append(a)
        r2s.append(results.rsquared)

.. code:: ipython3

    fig, ax = plt.subplots(1, 1)
    ax.plot(alphas, r2s, '.', label="observed")
    ax.plot(alphas, [a**2/(1+a**2) for a in alphas], label='theoretical')
    ax.set_xlabel("alpha")
    ax.set_ylabel("r2")
    ax.legend();



.. image:: regression_lineaire_13_0.png


Dans ce cas de régression simple, la valeur à prédire est :math:`y_i`,
la valeur prédite est :math:`\hat{y_i}=\alpha X_{1i}` et la moyenne
:math:`\bar{y} = \alpha \bar{X_1} + \bar{\epsilon} = 0`.

.. math:: R^2 = 1 - \frac{\sum_{i=1}^n (\hat{y_i}-\bar{y})^2}{\sum_{i=1}^n (y_i - \bar{y})^2}=1-\frac{\mathbb{V}\epsilon}{\alpha^2\mathbb{V}X_1+\mathbb{V}\epsilon} = 1 - \frac{1}{1+\alpha^2}=\frac{\alpha^2}{1+\alpha^2}

Deux variables corrélées
------------------------

On ne change pas le modèle mais on fait en sorte que :math:`X_2=X_1`.
Les deux variables sont corrélées.

.. code:: ipython3

    X[:, 1] = X[:, 0]
    Y = 2*X[:,0] + X[:, 2]
    model = OLS(Y,X[:, :2])
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.810</td>
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.810</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   4271.</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 27 Nov 2017</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> 
    </tr>
    <tr>
      <th>Time:</th>                 <td>12:06:03</td>     <th>  Log-Likelihood:    </th> <td> -1411.2</td>
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   2824.</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   999</td>      <th>  BIC:               </th> <td>   2829.</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     1</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>x1</th> <td>    1.0288</td> <td>    0.016</td> <td>   65.349</td> <td> 0.000</td> <td>    0.998</td> <td>    1.060</td>
    </tr>
    <tr>
      <th>x2</th> <td>    1.0288</td> <td>    0.016</td> <td>   65.349</td> <td> 0.000</td> <td>    0.998</td> <td>    1.060</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td> 8.165</td> <th>  Durbin-Watson:     </th> <td>   1.944</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.017</td> <th>  Jarque-Bera (JB):  </th> <td>   6.024</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td>-0.064</td> <th>  Prob(JB):          </th> <td>  0.0492</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 2.642</td> <th>  Cond. No.          </th> <td>1.61e+16</td>
    </tr>
    </table>



.. code:: ipython3

    model.rank




.. parsed-literal::
    1



Les variables corrélées n’ont pas l’air de déranger l’algorithme de
résolution car il utilise la méthode
`SVD <https://en.wikipedia.org/wiki/Singular-value_decomposition>`__
pour résoudre le même problème dans un espace de moindre dimension. Le
problème survient que les deux variables ne sont pas complétement
corrélées. On étudie le modèle :math:`Y \sim X_1 + X'_2` avec
:math:`X'_2 = \alpha X_1 + (1-\alpha) X_2` et on réduit la variance du
bruit pour en diminuer les effets.

.. code:: ipython3

    X_ = npr.normal(size=(1000, 3))

.. code:: ipython3

    alphas = [0.9 + i * 0.01 for i in range(0,11)]
    res = []
    for a in alphas:
        X = X_.copy()
        X[:, 1] = a * X[:, 0] + (1-a) * X[:, 1]
        Y = X[:, 0] + X[:, 1] + 0.1 * X[:, 2]
        model = OLS(Y,X[:, :2])
        results = model.fit()
        res.append(dict(alpha=a, r2=results.rsquared, rank=model.rank, c1=results.params[0], c2=results.params[1]))
        
    import pandas
    df = pandas.DataFrame(res)
    df = df.set_index('alpha')
    df






.. 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>c1</th>
          <th>c2</th>
          <th>r2</th>
          <th>rank</th>
        </tr>
        <tr>
          <th>alpha</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.90</th>
          <td>0.976370</td>
          <td>1.028982</td>
          <td>0.997391</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.91</th>
          <td>0.973150</td>
          <td>1.032202</td>
          <td>0.997416</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.92</th>
          <td>0.969125</td>
          <td>1.036227</td>
          <td>0.997440</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.93</th>
          <td>0.963950</td>
          <td>1.041402</td>
          <td>0.997464</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.94</th>
          <td>0.957049</td>
          <td>1.048303</td>
          <td>0.997489</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.95</th>
          <td>0.947389</td>
          <td>1.057963</td>
          <td>0.997513</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.96</th>
          <td>0.932898</td>
          <td>1.072454</td>
          <td>0.997536</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.97</th>
          <td>0.908747</td>
          <td>1.096605</td>
          <td>0.997560</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.98</th>
          <td>0.860444</td>
          <td>1.144908</td>
          <td>0.997583</td>
          <td>2</td>
        </tr>
        <tr>
          <th>0.99</th>
          <td>0.715536</td>
          <td>1.289816</td>
          <td>0.997606</td>
          <td>2</td>
        </tr>
        <tr>
          <th>1.00</th>
          <td>1.001279</td>
          <td>1.001279</td>
          <td>0.997627</td>
          <td>1</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1,2, figsize=(10,4))
    df[["r2"]].plot(ax=ax[0])
    df[["c1", "c2"]].plot(ax=ax[1])
    ax[0].set_title("R2")
    ax[1].set_title("coefficients");



.. image:: regression_lineaire_21_0.png


Le :math:`r^2` augmente quand la corrélation augmente mais les
coefficients sont moins fiables. Les résultats devraient être
sensiblement identiques en théorie mais en pratique, plus le déterminant
devient proche de zéro, plus l’ordinateur est limité par sa précision
numérique. Pour en savoir plus, vous pouvez lire un examen écrit que
j’ai rédigé, en python bien sûr : `Examen Programmation ENSAE première
année
2006 <http://www.xavierdupre.fr/site2013/enseignements/tdnote/ecrit_2006.pdf>`__.
Cette précision est aux alentours de :math:`10^{-15}` ce qui correspond
à la précision numérique des
`double <https://en.wikipedia.org/wiki/Double-precision_floating-point_format>`__.

.. code:: ipython3

    alphas = [1 - 10**(-i) for i in range(10,18)]
    res = []
    for a in alphas:
        X = X_.copy()
        X[:, 1] = a * X[:, 0] + (1-a) * X[:, 1]
        Y = X[:, 0] + X[:, 1] + X[:, 2]
        model = OLS(Y,X[:, :2])
        results = model.fit()
        res.append(dict(alpha_1=a-1, r2=results.rsquared, rank=model.rank, c1=results.params[0], c2=results.params[1]))
        
    import pandas
    df = pandas.DataFrame(res)
    df = df.set_index('alpha_1')
    df






.. 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>c1</th>
          <th>c2</th>
          <th>r2</th>
          <th>rank</th>
        </tr>
        <tr>
          <th>alpha_1</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>-1.000000e-10</th>
          <td>-2.898180e+08</td>
          <td>2.898180e+08</td>
          <td>0.811519</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-1.000000e-11</th>
          <td>-2.898201e+09</td>
          <td>2.898201e+09</td>
          <td>0.811519</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-9.999779e-13</th>
          <td>-2.898941e+10</td>
          <td>2.898941e+10</td>
          <td>0.811519</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-1.000311e-13</th>
          <td>-2.891422e+11</td>
          <td>2.891422e+11</td>
          <td>0.811518</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-9.992007e-15</th>
          <td>-2.915101e+12</td>
          <td>2.915101e+12</td>
          <td>0.811508</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-9.992007e-16</th>
          <td>1.012789e+00</td>
          <td>1.012789e+00</td>
          <td>0.811359</td>
          <td>2</td>
        </tr>
        <tr>
          <th>-1.110223e-16</th>
          <td>1.012789e+00</td>
          <td>1.012789e+00</td>
          <td>0.811359</td>
          <td>1</td>
        </tr>
        <tr>
          <th>0.000000e+00</th>
          <td>1.012789e+00</td>
          <td>1.012789e+00</td>
          <td>0.811359</td>
          <td>1</td>
        </tr>
      </tbody>
    </table>
    </div>



On fait un dernier test avec
`scikit-learn <http://scikit-learn.org/stable/>`__ pour vérifier que
l’algorithme de résolution donne des résultats similaires pour un cas où
le déterminant est quasi-nul.

.. code:: ipython3

    from sklearn.linear_model import LinearRegression
    from sklearn.metrics import r2_score
    
    alphas = [0.9 + i * 0.01 for i in range(0,11)]
    res = []
    for a in alphas:
        X = X_.copy()
        X[:, 1] = a * X[:, 0] + (1-a) * X[:, 1]
        Y = X[:, 0] + X[:, 1] + X[:, 2]
        model = LinearRegression()
        model.fit(X[:, :2], Y)
        r2 = r2_score(Y, model.predict(X[:, :2]))
        res.append(dict(alpha=a, c1=model.coef_[0], c2=model.coef_[1], r2=r2))
        
    import pandas
    df = pandas.DataFrame(res)
    df = df.set_index('alpha')
    df






.. 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>c1</th>
          <th>c2</th>
          <th>r2</th>
        </tr>
        <tr>
          <th>alpha</th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.90</th>
          <td>7.601931e-01</td>
          <td>1.293903e+00</td>
          <td>0.796916</td>
        </tr>
        <tr>
          <th>0.91</th>
          <td>7.275372e-01</td>
          <td>1.326559e+00</td>
          <td>0.798417</td>
        </tr>
        <tr>
          <th>0.92</th>
          <td>6.867173e-01</td>
          <td>1.367379e+00</td>
          <td>0.799911</td>
        </tr>
        <tr>
          <th>0.93</th>
          <td>6.342346e-01</td>
          <td>1.419862e+00</td>
          <td>0.801399</td>
        </tr>
        <tr>
          <th>0.94</th>
          <td>5.642576e-01</td>
          <td>1.489839e+00</td>
          <td>0.802880</td>
        </tr>
        <tr>
          <th>0.95</th>
          <td>4.662898e-01</td>
          <td>1.587807e+00</td>
          <td>0.804355</td>
        </tr>
        <tr>
          <th>0.96</th>
          <td>3.193382e-01</td>
          <td>1.734758e+00</td>
          <td>0.805823</td>
        </tr>
        <tr>
          <th>0.97</th>
          <td>7.441878e-02</td>
          <td>1.979678e+00</td>
          <td>0.807283</td>
        </tr>
        <tr>
          <th>0.98</th>
          <td>-4.154200e-01</td>
          <td>2.469516e+00</td>
          <td>0.808736</td>
        </tr>
        <tr>
          <th>0.99</th>
          <td>-1.884936e+00</td>
          <td>3.939033e+00</td>
          <td>0.810182</td>
        </tr>
        <tr>
          <th>1.00</th>
          <td>8.512221e+13</td>
          <td>-8.512221e+13</td>
          <td>0.811404</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1,3, figsize=(12,4))
    df[["c1", "c2"]].plot(ax=ax[1])
    df[["c1", "c2"]].plot(ax=ax[2])
    df[["r2"]].plot(ax=ax[0])
    ax[0].set_title("R2")
    ax[1].set_title("coefficients")
    ax[2].set_ylim([-5, 5])
    ax[2].set_title("coefficients, échelle tronquée");



.. image:: regression_lineaire_26_0.png


Le second graphe est trompeur mais il ne faut pas oublier de regarder
l’échelle de l’axe des ordonnées.

Indicatrices
------------

:math:`X_1` est une variable aléatoire gaussienne. On teste maintenant
un modèle :math:`Y = X'_1 + X'_2 + \epsilon` avec
:math:`X'_1 = X_1 \mathbf{1\!\!1}_{X_1 < 0}` et
:math:`X'_2 = X_1 \mathbf{1\!\!1}_{X_1 \geqslant 0}`.

.. code:: ipython3

    X = npr.normal(size=(1000, 3))
    X[:, 1] = X[:, 0]
    X[X[:, 0] >= 0, 0] = 0
    X[X[:, 1] < 0, 1] = 0
    Y = X[:, 0] + X[:, 1] + X[:, 2]
    corrcoef(X.T)




.. parsed-literal::
    array([[ 1.        ,  0.47358312, -0.03083914],
           [ 0.47358312,  1.        , -0.01293737],
           [-0.03083914, -0.01293737,  1.        ]])



.. code:: ipython3

    from pandas import DataFrame
    names = ["X%d" % i for i in range(X.shape[1]-1)]
    ax = DataFrame(X[:50,:2], columns=names).sort_values(names).reset_index(drop=True).plot()
    ax.set_title("Représentation des features tronquées");



.. image:: regression_lineaire_30_0.png


.. code:: ipython3

    model = OLS(Y,X[:, :3])
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th>  <td>   1.000</td> 
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   1.000</td> 
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>2.212e+33</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th>   <td>  0.00</td>  
    </tr>
    <tr>
      <th>Time:</th>                 <td>10:56:29</td>     <th>  Log-Likelihood:    </th>  <td>  33713.</td> 
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>-6.742e+04</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   997</td>      <th>  BIC:               </th> <td>-6.740e+04</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     3</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>x1</th> <td>    1.0000</td> <td> 2.42e-17</td> <td> 4.14e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    <tr>
      <th>x2</th> <td>    1.0000</td> <td> 2.39e-17</td> <td> 4.18e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    <tr>
      <th>x3</th> <td>    1.0000</td> <td> 1.73e-17</td> <td> 5.78e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td> 4.249</td> <th>  Durbin-Watson:     </th> <td>   2.031</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.119</td> <th>  Jarque-Bera (JB):  </th> <td>   4.338</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td>-0.107</td> <th>  Prob(JB):          </th> <td>   0.114</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 3.242</td> <th>  Cond. No.          </th> <td>    1.40</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



On découpe en trois.

.. code:: ipython3

    import numpy
    X = npr.normal(size=(1000, 4))
    for i in range(0, 3):
        X[:, i] = X_[:, 0]
    X[:, 3] = X_[:, 2]
    X[X_[:, 0] > -1, 0] = 0 
    X[(X_[:, 0] < -1) | (X_[:, 0] > 1), 1] = 0 
    X[X_[:, 0] < 1, 2] = 0 
    Y = X[:, 0] + X[:, 1] + X[:, 2] + X[:, 3]
    corrcoef(X.T)




.. parsed-literal::
    array([[ 1.        , -0.00347584,  0.16846101,  0.06722762],
           [-0.00347584,  1.        ,  0.00326437, -0.04707208],
           [ 0.16846101,  0.00326437,  1.        ,  0.08754832],
           [ 0.06722762, -0.04707208,  0.08754832,  1.        ]])



.. code:: ipython3

    from pandas import DataFrame
    names = ["X%d" % i for i in range(X.shape[1]-1)]
    ax = DataFrame(X[:50,:3], columns=names).sort_values(names).reset_index(drop=True).plot()
    ax.set_title("Représentation des features tronquées");



.. image:: regression_lineaire_34_0.png


.. code:: ipython3

    model = OLS(Y,X[:, :4])
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th>  <td>   1.000</td> 
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   1.000</td> 
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>1.910e+32</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th>   <td>  0.00</td>  
    </tr>
    <tr>
      <th>Time:</th>                 <td>10:57:27</td>     <th>  Log-Likelihood:    </th>  <td>  32608.</td> 
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>-6.521e+04</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   996</td>      <th>  BIC:               </th> <td>-6.519e+04</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     4</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>x1</th> <td>    1.0000</td> <td> 8.75e-17</td> <td> 1.14e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    <tr>
      <th>x2</th> <td>    1.0000</td> <td> 1.22e-16</td> <td> 8.23e+15</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    <tr>
      <th>x3</th> <td>    1.0000</td> <td> 8.33e-17</td> <td>  1.2e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    <tr>
      <th>x4</th> <td>    1.0000</td> <td> 5.23e-17</td> <td> 1.91e+16</td> <td> 0.000</td> <td>    1.000</td> <td>    1.000</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td>457.967</td> <th>  Durbin-Watson:     </th> <td>   1.816</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1967.636</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td>-2.198</td>  <th>  Prob(JB):          </th> <td>    0.00</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 8.282</td>  <th>  Cond. No.          </th> <td>    2.35</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



Régression linéaire par morceaux
--------------------------------

On se place dans un cas particulier où le problème est linéaire par
morceaux :

.. math:: Y = -2 X_1 \mathbf{1\!\!1}_{X_1 + \epsilon_1 <0} + 4 X_1 \mathbf{1\!\!1}_{X + \epsilon_1 > 0} + \epsilon_2

La régression donne de très mauvais résultat sur ce type de problèmes
mais on cherche une façon systématique de découper le problème en
segments linéaires.

.. code:: ipython3

    X = npr.normal(size=(1000,4))
    alpha = [4, -2]
    t = (X[:, 0] + X[:, 3] * 0.5) > 0
    switch = numpy.zeros(X.shape[0])
    switch[t] = 1
    Y = alpha[0] * X[:, 0] * t + alpha[1] * X[:, 0] * (1-t) + X[:, 2]

.. code:: ipython3

    fig, ax = plt.subplots(1, 1)
    ax.plot(X[:, 0], Y, ".")
    ax.set_title("Nuage de points linéaire par morceaux");



.. image:: regression_lineaire_38_0.png


.. code:: ipython3

    model = OLS(Y,X[:, :1])
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.094</td>
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.093</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   104.0</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th> <td>2.69e-23</td>
    </tr>
    <tr>
      <th>Time:</th>                 <td>10:59:28</td>     <th>  Log-Likelihood:    </th> <td> -2594.9</td>
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   5192.</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   999</td>      <th>  BIC:               </th> <td>   5197.</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     1</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>x1</th> <td>    1.0252</td> <td>    0.101</td> <td>   10.197</td> <td> 0.000</td> <td>    0.828</td> <td>    1.222</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td> 3.882</td> <th>  Durbin-Watson:     </th> <td>   1.092</td>
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.144</td> <th>  Jarque-Bera (JB):  </th> <td>   3.834</td>
    </tr>
    <tr>
      <th>Skew:</th>          <td> 0.151</td> <th>  Prob(JB):          </th> <td>   0.147</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 3.015</td> <th>  Cond. No.          </th> <td>    1.00</td>
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



.. code:: ipython3

    yp = results.predict(X[:, :1])

.. code:: ipython3

    fig, ax = plt.subplots(1, 1)
    ax.plot(X[:, 0], Y, ".", label="expected")
    ax.plot(X[:, 0], yp, ".", label="predicted")
    ax.legend()
    ax.set_title("Régression linéaire sur un nuage linéaire par morceaux");



.. image:: regression_lineaire_41_0.png


Passons à un arbre de décision qui n’est pas le meilleur modèle mais on
va détourner ses résultats pour revenir à un problème de régression par
morceaux.

.. code:: ipython3

    from sklearn.tree import DecisionTreeRegressor
    model = DecisionTreeRegressor(min_samples_leaf=10, max_depth=3)
    model.fit(X[:, :1], Y)
    yp = model.predict(X[:, :1])

.. code:: ipython3

    fig, ax = plt.subplots(1, 1)
    ax.plot(X[:, 0], Y, ".", label="expected")
    ax.plot(X[:, 0], yp, ".", label="predicted")
    ax.legend()
    r2 = r2_score(Y, model.predict(X[:, :1]))
    ax.set_title("Arbre de décision sur un nuage linéaire par morceaux\nR2=%f" % r2);



.. image:: regression_lineaire_44_0.png


.. code:: ipython3

    from sklearn.tree import export_graphviz
    export_graphviz(model, out_file="arbre.dot")
    from pyensae.graphhelper import run_dot
    run_dot("arbre.dot", "arbre.png")
    from IPython.display import Image
    Image("arbre.png")




.. image:: regression_lineaire_45_0.png



On extrait tous les seuils de l’arbre et on ajoute les milieux de
segments.

.. code:: ipython3

    th = list(sorted(set(model.tree_.threshold)))
    th += [(th[i] + th[i-1])/2 for i in range(1,len(th))]
    th = list(sorted(th))
    th




.. parsed-literal::
    [-2.0,
     -1.8018612563610077,
     -1.6037225127220154,
     -1.323736995458603,
     -1.0437514781951904,
     -0.3109976723790169,
     0.4217561334371567,
     0.678125374019146,
     0.9344946146011353,
     1.0011553764343262,
     1.067816138267517,
     1.2776717841625214,
     1.4875274300575256,
     1.7147845923900604,
     1.9420417547225952]



On fait une régression sur les variables
:math:`W_{i>0} = X_1 \mathbf{1\!\!1}_{X_1 > t_i}`, :math:`W_0 = X_1` où les
:math:`(t_i)` sont les seuils.

.. code:: ipython3

    W = numpy.zeros((X.shape[0], len(th)+1))
    x = X[:, 0]
    W[:, 0] = x
    for i in range(len(th)):
        W[x > th[i], i+1] = x[x > th[i]]

.. code:: ipython3

    model = OLS(Y,W)
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.849</td>
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.847</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   346.9</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> 
    </tr>
    <tr>
      <th>Time:</th>                 <td>11:07:03</td>     <th>  Log-Likelihood:    </th> <td> -1697.7</td>
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   3427.</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   984</td>      <th>  BIC:               </th> <td>   3506.</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>    16</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>x1</th>  <td>   -1.8316</td> <td>    0.124</td> <td>  -14.806</td> <td> 0.000</td> <td>   -2.074</td> <td>   -1.589</td>
    </tr>
    <tr>
      <th>x2</th>  <td>   -0.0573</td> <td>    0.210</td> <td>   -0.273</td> <td> 0.785</td> <td>   -0.469</td> <td>    0.354</td>
    </tr>
    <tr>
      <th>x3</th>  <td>   -0.2440</td> <td>    0.235</td> <td>   -1.038</td> <td> 0.300</td> <td>   -0.705</td> <td>    0.217</td>
    </tr>
    <tr>
      <th>x4</th>  <td>    0.3291</td> <td>    0.218</td> <td>    1.508</td> <td> 0.132</td> <td>   -0.099</td> <td>    0.757</td>
    </tr>
    <tr>
      <th>x5</th>  <td>   -0.3610</td> <td>    0.206</td> <td>   -1.756</td> <td> 0.079</td> <td>   -0.764</td> <td>    0.042</td>
    </tr>
    <tr>
      <th>x6</th>  <td>    0.5528</td> <td>    0.197</td> <td>    2.811</td> <td> 0.005</td> <td>    0.167</td> <td>    0.939</td>
    </tr>
    <tr>
      <th>x7</th>  <td>    3.2628</td> <td>    0.395</td> <td>    8.253</td> <td> 0.000</td> <td>    2.487</td> <td>    4.039</td>
    </tr>
    <tr>
      <th>x8</th>  <td>    1.4566</td> <td>    0.449</td> <td>    3.244</td> <td> 0.001</td> <td>    0.576</td> <td>    2.338</td>
    </tr>
    <tr>
      <th>x9</th>  <td>    0.2701</td> <td>    0.311</td> <td>    0.869</td> <td> 0.385</td> <td>   -0.340</td> <td>    0.880</td>
    </tr>
    <tr>
      <th>x10</th> <td>    0.7213</td> <td>    0.374</td> <td>    1.928</td> <td> 0.054</td> <td>   -0.013</td> <td>    1.456</td>
    </tr>
    <tr>
      <th>x11</th> <td>   -0.4599</td> <td>    0.457</td> <td>   -1.006</td> <td> 0.315</td> <td>   -1.357</td> <td>    0.437</td>
    </tr>
    <tr>
      <th>x12</th> <td>    0.4177</td> <td>    0.378</td> <td>    1.105</td> <td> 0.269</td> <td>   -0.324</td> <td>    1.159</td>
    </tr>
    <tr>
      <th>x13</th> <td>   -0.2703</td> <td>    0.253</td> <td>   -1.069</td> <td> 0.285</td> <td>   -0.766</td> <td>    0.226</td>
    </tr>
    <tr>
      <th>x14</th> <td>    0.5325</td> <td>    0.226</td> <td>    2.360</td> <td> 0.018</td> <td>    0.090</td> <td>    0.975</td>
    </tr>
    <tr>
      <th>x15</th> <td>   -0.3703</td> <td>    0.229</td> <td>   -1.618</td> <td> 0.106</td> <td>   -0.819</td> <td>    0.079</td>
    </tr>
    <tr>
      <th>x16</th> <td>    0.1996</td> <td>    0.194</td> <td>    1.026</td> <td> 0.305</td> <td>   -0.182</td> <td>    0.581</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td>228.022</td> <th>  Durbin-Watson:     </th> <td>   1.978</td> 
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td> 716.200</td> 
    </tr>
    <tr>
      <th>Skew:</th>          <td>-1.110</td>  <th>  Prob(JB):          </th> <td>3.01e-156</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 6.502</td>  <th>  Cond. No.          </th> <td>    36.3</td> 
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



Dessinons les résultats de la prédictions.

.. code:: ipython3

    yp = results.predict(W)
    fig, ax = plt.subplots(1, 1)
    ax.plot(X[:, 0], Y, ".", label="expected")
    ax.plot(X[:, 0], yp, ".", label="predicted")
    ax.legend()
    ax.set_title("Régression linéaire par morceaux\nsur un nuage linéaire par morceaux\nR2=%f" % results.rsquared);



.. image:: regression_lineaire_52_0.png


Le modèle nous suggère de ne garder que quelques seuils. En s’appuyant
sur les p-values :

.. code:: ipython3

    keep = numpy.arange(len(results.pvalues))[results.pvalues < 0.05]
    keep




.. parsed-literal::
    array([ 0,  5,  6,  7, 13])



.. code:: ipython3

    W2 = W[:, keep]

.. code:: ipython3

    model = OLS(Y,W2)
    results = model.fit()
    results.summary()






.. raw:: html

    <table class="simpletable">
    <caption>OLS Regression Results</caption>
    <tr>
      <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.846</td>
    </tr>
    <tr>
      <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.845</td>
    </tr>
    <tr>
      <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1094.</td>
    </tr>
    <tr>
      <th>Date:</th>             <td>Mon, 15 Oct 2018</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> 
    </tr>
    <tr>
      <th>Time:</th>                 <td>11:07:38</td>     <th>  Log-Likelihood:    </th> <td> -1708.6</td>
    </tr>
    <tr>
      <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   3427.</td>
    </tr>
    <tr>
      <th>Df Residuals:</th>          <td>   995</td>      <th>  BIC:               </th> <td>   3452.</td>
    </tr>
    <tr>
      <th>Df Model:</th>              <td>     5</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>x1</th> <td>   -1.9504</td> <td>    0.066</td> <td>  -29.574</td> <td> 0.000</td> <td>   -2.080</td> <td>   -1.821</td>
    </tr>
    <tr>
      <th>x2</th> <td>    0.3384</td> <td>    0.148</td> <td>    2.287</td> <td> 0.022</td> <td>    0.048</td> <td>    0.629</td>
    </tr>
    <tr>
      <th>x3</th> <td>    3.2628</td> <td>    0.397</td> <td>    8.209</td> <td> 0.000</td> <td>    2.483</td> <td>    4.043</td>
    </tr>
    <tr>
      <th>x4</th> <td>    2.0247</td> <td>    0.385</td> <td>    5.260</td> <td> 0.000</td> <td>    1.269</td> <td>    2.780</td>
    </tr>
    <tr>
      <th>x5</th> <td>    0.4635</td> <td>    0.119</td> <td>    3.901</td> <td> 0.000</td> <td>    0.230</td> <td>    0.697</td>
    </tr>
    </table>
    <table class="simpletable">
    <tr>
      <th>Omnibus:</th>       <td>248.807</td> <th>  Durbin-Watson:     </th> <td>   1.984</td> 
    </tr>
    <tr>
      <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td> 829.417</td> 
    </tr>
    <tr>
      <th>Skew:</th>          <td>-1.190</td>  <th>  Prob(JB):          </th> <td>7.84e-181</td>
    </tr>
    <tr>
      <th>Kurtosis:</th>      <td> 6.774</td>  <th>  Cond. No.          </th> <td>    20.1</td> 
    </tr>
    </table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



.. code:: ipython3

    yp = results.predict(W2)
    fig, ax = plt.subplots(1, 1)
    ax.plot(X[:, 0], Y, ".", label="expected")
    ax.plot(X[:, 0], yp, ".", label="predicted")
    ax.legend()
    ax.set_title("Régression linéaire par morceaux\nsur un nuage linéaire par morceaux\n" +
                 "réduction du nombre de segments\nR2=%f" % results.rsquared);



.. image:: regression_lineaire_57_0.png


Le coefficient :math:`R^2` est quasiment identique pour un nombre de
segments moindre. Je me suis amusé à rendre ce code plus générique pour
comparer la première étape, le découpage en morceaux, via deux modèles,
un arbre de décision et le nouvel objet
`KBinsDiscretizer <https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.KBinsDiscretizer.html>`__
qui segmente une variable sans tenir compte de la cible. La régression
n’est plus nécessaire linéaire : `Piecewise linear
regression <http://www.xavierdupre.fr/app/mlinsights/helpsphinx/notebooks/piecewise_linear_regression.html>`__.
Je me suis également amusé à faire de même pour une classification par
morceaux
`PiecewiseClassifier <http://www.xavierdupre.fr/app/mlinsights/helpsphinx/mlinsights/mlmodel/piecewise_estimator.html?mlinsights.mlmodel.piecewise_estimator.PiecewiseClassifier>`__.
Celle-ci pose quelques soucis pratiques car toutes les classes ne sont
pas forcément représentées dans chaque compartiment…