.. _PCArst:

==================================
PCA (Principal Component Analysis)
==================================


.. only:: html

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


This notebook shows how to plot a PCA with sciki-learn and statsmodels,
with or without normalization.

.. code:: ipython3

    %matplotlib inline

.. code:: ipython3

    import matplotlib.pyplot as plt
    plt.style.use('ggplot')

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





More about
`PCA <https://en.wikipedia.org/wiki/Principal_component_analysis>`__:
`Implementing a Principal Component Analysis (PCA) in Python step by
step <http://sebastianraschka.com/Articles/2014_pca_step_by_step.html>`__.

Download data
-------------

.. code:: ipython3

    import pyensae.datasource
    pyensae.datasource.download_data("auto-mpg.data", url="https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/")




.. parsed-literal::
    'auto-mpg.data'



.. code:: ipython3

    import pandas
    df = pandas.read_fwf("auto-mpg.data", encoding="utf-8",
                        names="mpg cylinders displacement horsepower weight acceleration year origin name".split())
    df["name"] = df["name"].apply(lambda s : s.strip(' "'))
    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>mpg</th>
          <th>cylinders</th>
          <th>displacement</th>
          <th>horsepower</th>
          <th>weight</th>
          <th>acceleration</th>
          <th>year</th>
          <th>origin</th>
          <th>name</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>18.0</td>
          <td>8</td>
          <td>307.0</td>
          <td>130.0</td>
          <td>3504.0</td>
          <td>12.0</td>
          <td>70</td>
          <td>1</td>
          <td>chevrolet chevelle malibu</td>
        </tr>
        <tr>
          <th>1</th>
          <td>15.0</td>
          <td>8</td>
          <td>350.0</td>
          <td>165.0</td>
          <td>3693.0</td>
          <td>11.5</td>
          <td>70</td>
          <td>1</td>
          <td>buick skylark 320</td>
        </tr>
        <tr>
          <th>2</th>
          <td>18.0</td>
          <td>8</td>
          <td>318.0</td>
          <td>150.0</td>
          <td>3436.0</td>
          <td>11.0</td>
          <td>70</td>
          <td>1</td>
          <td>plymouth satellite</td>
        </tr>
        <tr>
          <th>3</th>
          <td>16.0</td>
          <td>8</td>
          <td>304.0</td>
          <td>150.0</td>
          <td>3433.0</td>
          <td>12.0</td>
          <td>70</td>
          <td>1</td>
          <td>amc rebel sst</td>
        </tr>
        <tr>
          <th>4</th>
          <td>17.0</td>
          <td>8</td>
          <td>302.0</td>
          <td>140.0</td>
          <td>3449.0</td>
          <td>10.5</td>
          <td>70</td>
          <td>1</td>
          <td>ford torino</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    df.dtypes




.. parsed-literal::
    mpg             float64
    cylinders         int64
    displacement    float64
    horsepower       object
    weight          float64
    acceleration    float64
    year              int64
    origin            int64
    name             object
    dtype: object



We remove missing values:

.. code:: ipython3

    df[df.horsepower == "?"]






.. 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>mpg</th>
          <th>cylinders</th>
          <th>displacement</th>
          <th>horsepower</th>
          <th>weight</th>
          <th>acceleration</th>
          <th>year</th>
          <th>origin</th>
          <th>name</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>32</th>
          <td>25.0</td>
          <td>4</td>
          <td>98.0</td>
          <td>?</td>
          <td>2046.0</td>
          <td>19.0</td>
          <td>71</td>
          <td>1</td>
          <td>ford pinto</td>
        </tr>
        <tr>
          <th>126</th>
          <td>21.0</td>
          <td>6</td>
          <td>200.0</td>
          <td>?</td>
          <td>2875.0</td>
          <td>17.0</td>
          <td>74</td>
          <td>1</td>
          <td>ford maverick</td>
        </tr>
        <tr>
          <th>330</th>
          <td>40.9</td>
          <td>4</td>
          <td>85.0</td>
          <td>?</td>
          <td>1835.0</td>
          <td>17.3</td>
          <td>80</td>
          <td>2</td>
          <td>renault lecar deluxe</td>
        </tr>
        <tr>
          <th>336</th>
          <td>23.6</td>
          <td>4</td>
          <td>140.0</td>
          <td>?</td>
          <td>2905.0</td>
          <td>14.3</td>
          <td>80</td>
          <td>1</td>
          <td>ford mustang cobra</td>
        </tr>
        <tr>
          <th>354</th>
          <td>34.5</td>
          <td>4</td>
          <td>100.0</td>
          <td>?</td>
          <td>2320.0</td>
          <td>15.8</td>
          <td>81</td>
          <td>2</td>
          <td>renault 18i</td>
        </tr>
        <tr>
          <th>374</th>
          <td>23.0</td>
          <td>4</td>
          <td>151.0</td>
          <td>?</td>
          <td>3035.0</td>
          <td>20.5</td>
          <td>82</td>
          <td>1</td>
          <td>amc concord dl</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    final = df[df.horsepower != '?'].copy()
    final["horsepower"] = final["horsepower"].astype(float)

.. code:: ipython3

    final.to_csv("auto-mpg.data.csv", sep="\t", index=False, encoding="utf-8")

.. code:: ipython3

    final.shape




.. parsed-literal::
    (392, 9)



PCA with scikit-learn
---------------------

.. code:: ipython3

    from sklearn.decomposition import PCA
    X = final[df.columns[1:-1]]
    Y = final["mpg"]
    pca = PCA(n_components=2)
    pca.fit(X)




.. parsed-literal::
    PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
      svd_solver='auto', tol=0.0, whiten=False)



.. code:: ipython3

    out = pca.transform(X)
    out[:5]




.. parsed-literal::
    array([[536.44492922,  50.83312832],
           [730.34140206,  79.13543921],
           [470.9815846 ,  75.4476722 ],
           [466.40143367,  62.53420646],
           [481.66788465,  55.78036021]])



.. code:: ipython3

    pca.explained_variance_ratio_, pca.noise_variance_




.. parsed-literal::
    (array([0.99756151, 0.0020628 ]), 55.14787750463889)



.. code:: ipython3

    import matplotlib.pyplot as plt
    plt.plot(out[:,0], out[:,1], ".");



.. image:: PCA_18_0.png


.. code:: ipython3

    pca.components_ 




.. parsed-literal::
    array([[ 1.79262233e-03,  1.14341275e-01,  3.89670355e-02,
             9.92673415e-01, -1.35283460e-03, -1.33684138e-03,
            -5.51538021e-04],
           [ 1.33244815e-02,  9.45778439e-01,  2.98248416e-01,
            -1.20752748e-01, -3.48258394e-02, -2.38516836e-02,
            -3.24298106e-03]])



PCA with scikit-learn and normalization
---------------------------------------

.. code:: ipython3

    from sklearn.decomposition import PCA
    from sklearn.preprocessing import Normalizer
    from sklearn.pipeline import Pipeline
    
    normpca = Pipeline([('norm', Normalizer()), ('pca', PCA(n_components=2))])
    normpca.fit(X)




.. parsed-literal::
    Pipeline(memory=None,
         steps=[('norm', Normalizer(copy=True, norm='l2')), ('pca', PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
      svd_solver='auto', tol=0.0, whiten=False))])



.. code:: ipython3

    out = normpca.transform(X)
    out[:5]




.. parsed-literal::
    array([[0.02731781, 0.00012872],
           [0.03511968, 0.00666259],
           [0.03247168, 0.00632048],
           [0.0287677 , 0.0060517 ],
           [0.02758449, 0.00325874]])



.. code:: ipython3

    normpca.named_steps['pca'].explained_variance_ratio_, normpca.named_steps['pca'].noise_variance_




.. parsed-literal::
    (array([0.86819249, 0.08034075]), 4.332607718595102e-06)



.. code:: ipython3

    import matplotlib.pyplot as plt
    plt.plot(out[:,0], out[:,1], ".");



.. image:: PCA_24_0.png


.. code:: ipython3

    normpca.named_steps['pca'].components_




.. parsed-literal::
    array([[ 0.00415209,  0.92648229,  0.11272098, -0.05732771, -0.09162071,
            -0.34198745, -0.01646403],
           [ 0.01671457,  0.0789351 ,  0.85881718, -0.06957932,  0.02998247,
             0.49941847,  0.02763848]])



PCA with statsmodels
--------------------

.. code:: ipython3

    from statsmodels.sandbox.tools import pca
    xred, fact, eva, eve = pca(X, keepdim=2, normalize=False)

.. code:: ipython3

    fact[:5]




.. parsed-literal::
    array([[536.44492922, -50.83312832],
           [730.34140206, -79.13543921],
           [470.9815846 , -75.4476722 ],
           [466.40143367, -62.53420646],
           [481.66788465, -55.78036021]])



.. code:: ipython3

    eva




.. parsed-literal::
    array([732151.6743476 ,   1513.97202164])



.. code:: ipython3

    eve




.. parsed-literal::
    array([[ 1.79262233e-03, -1.33244815e-02],
           [ 1.14341275e-01, -9.45778439e-01],
           [ 3.89670355e-02, -2.98248416e-01],
           [ 9.92673415e-01,  1.20752748e-01],
           [-1.35283460e-03,  3.48258394e-02],
           [-1.33684138e-03,  2.38516836e-02],
           [-5.51538021e-04,  3.24298106e-03]])



.. code:: ipython3

    plt.plot(fact[:,0], fact[:,1], ".");



.. image:: PCA_31_0.png


PCA with statsmodels and normalization
--------------------------------------

.. code:: ipython3

    from statsmodels.sandbox.tools import pca
    from sklearn.preprocessing import normalize
    X_norm = normalize(X)
    xred, fact, eva, eve = pca(X_norm, keepdim=2, normalize=True)

.. code:: ipython3

    eva




.. parsed-literal::
    array([3.65433661e-04, 3.38164814e-05])



.. code:: ipython3

    eve




.. parsed-literal::
    array([[ -0.21720145,   2.87429329],
           [-48.46551687,  13.57394009],
           [ -5.89658384, 147.68504393],
           [  2.99888854, -11.96508998],
           [  4.79280102,   5.15588534],
           [ 17.88981698,  85.8816515 ],
           [  0.86125514,   4.75280519]])



.. code:: ipython3

    plt.plot(fact[:,0], fact[:,1], ".");



.. image:: PCA_36_0.png