10-dimensionality-reduction-with-pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective\n",
    "\n",
    "This notebook aims to implement **principal component analysis (PCA)** on the decathlon's data, which refers to men's  performance during two athletic meetings: [2004 Summer Olympics]\n",
    "(https://en.wikipedia.org/wiki/Athletics_at_the_2004_Summer_Olympics_%E2%80%93_Men%27s_decathlon) and 2004 [Decastar](https://en.wikipedia.org/wiki/Décastar)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Decathlon's data set\n",
    "\n",
    "The data are available at ``data/decathlon.txt``\n",
    "\n",
    "### Data set description\n",
    "\n",
    "* The dataset comprises 41 rows and 13 columns\n",
    "* The first row is a header describing the content of the columns and the remaining rows refer to the 40 observations (i.e., athletes) considered in this dataset\n",
    "* Columns 1 to 12 are continuous variables: the first ten columns correspond to the performance of the athletes for each event of the decathlon\n",
    "* Columns 11 and 12 correspond respectively to the rank and the obtained points\n",
    "* The last column is a categorical variable corresponding to the athletic meeting (2004 Olympic Games or 2004 Decastar)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the shape of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the name of the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the first five records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "decathlon.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Selecting a column by name. Observer that the returned object is also a pandas object (a *series*--a single-columned DataFrame). \n",
    "This means that we can use the `head()` function to view the first few rows only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "decathlon['100m'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also select a list of columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon[['100m', '400m']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the unique values of a feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon['Competition'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select only the rows that satisfy a given condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon[decathlon['Competition']=='OlympicG'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To *index* a row, we can use the data frame's `loc` object. This behaves like a dictionary whose keys are the data frame's index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon.loc['CLAY']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get an overview of the data through the method ```describe```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "decathlon.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Manipulating the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X = decathlon.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the types of a dataframe through its attribute ``dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decathlon.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__TODO__ Visualize a scatter plot of the athletes' performances regarding two disciplines \n",
    "(e.g., 400m vs Shot.put).\n",
    "\n",
    "_Hint_: check the method `plot` of pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting a **scatterplot** matrix of the variables.\n",
    "\n",
    "A scatterplot matrix enables us to visualize:\n",
    "\n",
    "* the density estimation for each of the features on main the diagonal\n",
    "* a scatterplot between two of the features on each of the off-diagonal\n",
    "\n",
    "In this case, each dot represents a samples "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas.plotting import scatter_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data cleaning\n",
    "\n",
    "We don't need the features _Points_, _Rank_,  and _Competition_ to represent the athletes. Thus, we can remove them through method ``drop`` (i.e., `dataframe.drop`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Implementing principal components analysis (PCA) to find the two principal components of the data\n",
    "\n",
    "Assuming that we are performing PCA on some dataset $\\boldsymbol X$ for $M$ principal components.\n",
    "\n",
    "The PCA of $\\mathbf{X} \\in \\mathbb{R}^{N \\times D}$ comprises in finding the $M$ principal components as the column of $\\mathbf{W} \\in \\mathbb{R}^{D \\times M}$ and a matrix of projected data $\\mathbf{Z} \\in \\mathbb{R}^{N \\times M}$, such that $\\mathbf{Z}\\mathbf{W}^T$ minimises the error as a reconstruction of $\\mathbf{X}$ for the choice of $M$; i.e., it is equivalently to the factorisation maximally explaining the variance in $\\mathbf{X}$. \n",
    "\n",
    "The principal components $\\mathbf{W}$, correspond to the $M$ largest eigenvectors of the empirical covariance matrix $$\\boldsymbol\\Sigma = \\frac{1}{N}\\mathbf{X}^T\\mathbf{X}$$\n",
    "\n",
    "Before computing the $M$ principal components, we need to:\n",
    "\n",
    "1. Normalize the data (i.e., mean $\\mu = 0$ and standard deviation $\\sigma = 1$). This step is known as [data normalization or feature scaling](https://en.wikipedia.org/wiki/Feature_scaling)\n",
    "2. Find the eigenvalues and the corresponding eigenvectors for the covariance matrix $S$\n",
    "3. Sort by the largest eigenvalues and the corresponding eigenvectors (`eig`)\n",
    "4. Compute the projection and the reconstruction of the data onto the spaced spanned by the top $n$ eigenvectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linalg\n",
    "\n",
    "def normalize(X):\n",
    "    \n",
    "    \"\"\"Normalize the given dataset X\n",
    "    \n",
    "    Args:\n",
    "        X: ndarray, data set\n",
    "    \n",
    "    Returns:\n",
    "        (Xbar, mean, std): tuple of ndarray, Xbar is the normalized dataset with mean 0 and standard deviation 1; \n",
    "        mean and std are respectively the mean and standard deviation.\n",
    "    \n",
    "    Note:\n",
    "        Dimensions where the standard deviation is zero their normalized data will be NaN. \n",
    "        This function handles these cases by setting `std = 1`.\n",
    "    \"\"\"\n",
    "    \n",
    "    mu = None\n",
    "    std = None\n",
    "    Xbar = None\n",
    "    \n",
    "    return Xbar, mu, std\n",
    "\n",
    "\n",
    "def eig(S):\n",
    "    \n",
    "    \"\"\"Compute the eigenvalues and corresponding eigenvectors for the covariance matrix S.\n",
    "    \n",
    "    Args:\n",
    "        S: ndarray, covariance matrix\n",
    "    \n",
    "    Returns:\n",
    "        (eigvals, eigvecs): ndarray, the eigenvalues and eigenvectors\n",
    "\n",
    "    Note:\n",
    "        the eigenvals and eigenvectors are sorted in descending order of the eigenvalues\n",
    "    \"\"\"\n",
    "    eigvals, eigvals = None, None\n",
    "    \n",
    "    return (eigvals, eigvecs)\n",
    "\n",
    "\n",
    "def projection_matrix(B):\n",
    "    \n",
    "    \"\"\"Compute the projection matrix onto the space spanned by `B`\n",
    "    \n",
    "    Args:\n",
    "        B: ndarray of dimension (D, M), the basis for the subspace\n",
    "    \n",
    "    Returns:\n",
    "        P: the projection matrix\n",
    "    \"\"\"\n",
    "    return (B @ np.linalg.inv(B.T @ B) @ B.T)\n",
    "\n",
    "\n",
    "def PCA(X, n_components):\n",
    "    \n",
    "    \"\"\" Compute the principal components of given data set.\n",
    "    \n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the data, and N is the number of data points\n",
    "        n_components: the number of principal components to use.\n",
    "    Returns:\n",
    "        X_reconstructed: ndarray of the reconstruction of X from the first `n_components` principal components.\n",
    "    \"\"\"\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__TODO:__ compute the two principal components of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = new_decathlon.values\n",
    "X_projected = PCA(X, n_components = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us display the projected data. We use a jupyter magic command to display plots inline automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create figure and axis objects\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "\n",
    "# create scatterplot on axis N.B. we record the return value to feed to the colorbar\n",
    "cax = ax.scatter(X_projected[:, 0], \n",
    "                 X_projected[:, 1], \n",
    "                 c = decathlon['Rank'], \n",
    "                 cmap=plt.get_cmap('viridis'))\n",
    "\n",
    "ax.set_xlim([-5.5, 5.5])\n",
    "ax.set_ylim([-4, 4])\n",
    "\n",
    "plt.colorbar(cax, label='Rank');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Can you identity a pattern?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Answer__:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Using scikit-learn to find the principal components\n",
    "\n",
    "Use scikit-learn to compute the principal components, and compare the results to what we have gotten so far. Check the online documentation available on: http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition, preprocessing\n",
    "\n",
    "X_scaled = preprocessing.StandardScaler().fit_transform(X)\n",
    "\n",
    "pca = decomposition.PCA(n_components = 10)\n",
    "pca.fit(X_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `pca.transform` to project the data onto its principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_projected = pca.transform(X_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pca.explained_variance_ratio_` gives the percentage of variance explained by each of the components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pca.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** How is `pca.explained_variance_ratio_` computed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Project the data onto its principal components after standardizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "\n",
    "pcs = pca.components_\n",
    "\n",
    "plt.scatter(pcs[0, :], pcs[1, :])\n",
    "for i, (x, y) in enumerate(zip(pcs[0, :], pcs[1, :])):\n",
    "    plt.text((x+0.01), (y+0.01), decathlon.columns[i], fontsize='9')\n",
    "\n",
    "# Plot a horizontal line y=0\n",
    "plt.plot([-0.7, 0.7], [0, 0], color='grey', ls='--')\n",
    "\n",
    "# Plot a vertical line x=0\n",
    "plt.plot([0, 0], [-0.7, 0.7], color='grey', ls='--')\n",
    "\n",
    "plt.xlim([-0.7, 0.7])\n",
    "plt.ylim([-0.7, 0.7])\n",
    "\n",
    "plt.xlabel('Weight of the contribution to PC 1', fontsize=12)\n",
    "plt.ylabel('Weight of the contribution to PC 2', fontsize=12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.bar(np.arange(10), pca.explained_variance_ratio_, color='blue')\n",
    "plt.xlim([-1, 9])\n",
    "plt.xlabel(\"Number of PCs\", fontsize=10)\n",
    "plt.ylabel(\"Fraction of variance explained\", fontsize=10)\n",
    "plt.xticks(np.arange(0, 10), np.arange(1, 11));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To better understand the information captured by the principal components, we can consider  `pca.components_`. These are the columns of $\\mathbf{W}$ (for $M = 2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcs = pca.components_\n",
    "print(pcs[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display each row of $\\mathbf{W}$ in a 2D plot whose x-axis gives its contribution to the first component and y-axis to the second component. Note, whereas before we were visualising the projected data, $\\mathbf{Z}$, now we are visualising the projections, $\\mathbf{W}$. This indicates how the features cluster i.e. if a pair of feature projections are close, observations will tend to be similarly-valued over those features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 5))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_xlim([-0.7, 0.7])\n",
    "ax.set_ylim([-0.7, 0.7])\n",
    "\n",
    "for i, (x, y) in enumerate(zip(pcs[0, :], pcs[1, :])):\n",
    "    # plot line between origin and point (x, y)\n",
    "    ax.plot([0, x], [0, y], color='k')\n",
    "    # display the label of the point\n",
    "    ax.text(x, y, decathlon.columns[i], fontsize='9')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Based on the two previous graphs, can you find a meaning for the two principal components?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}