code/PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal Components Analysis\n",
    "\n",
    "In this notebook, we shall see how we can use PCA to reduce dimensionality. Also, with the included plots, hopefully an intuition as to\n",
    "how it works can be built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the required libraries\n",
    "import numpy as np\n",
    "import sklearn as sk\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We'll create some data and see what PCA discovers in them. For this, we will draw samples from a Gaussian distribution\n",
    "# Our data is centered about origin\n",
    "mean = [0, 0, 0, 0]\n",
    "# Our data will have covariance as so. See that the last two features are correlated!\n",
    "cov = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 4, 4], [0, 0, 4, 4]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(500, 4)\n"
     ]
    }
   ],
   "source": [
    "# We will generate 500 points lying in 4D space, generated by the gaussian distribution with given mean and covariance\n",
    "X = np.random.multivariate_normal(mean, cov, 500)\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that our data is composed of 4 \"features\" now. The last two have higher variance (4), and are also correlated. That means if we do PCA, we should\n",
    "\"discover\" 1 direction along which the spread is maximum, then two directions which have reduced spread (first two features are unrelated, and don't spread), \n",
    "and the fourth direction will not have any spread (since features are perfectly correlated!). Feel free to play with different settings of the Gaussian distribution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,\n",
       "  svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4lHed9/H3N5MT4RQgoUBICGcIlQKN2INWalsLbQEf\nt6vg5aG7rn2q4mmr+9BHbYF6rodrt8uqWPVRdxW73XUNLS21Z6ttJRxKTULaQDmEY8r5HBK+zx8z\nwBCGZAiT3HP4vK4rV+ae+WXme3fKZ35zzz3fn7k7IiKSXrKCLkBERBJP4S4ikoYU7iIiaUjhLiKS\nhhTuIiJpSOEuIpKGFO4iImlI4S4ikoYU7iIiaSg7qAcuKiry8vLyoB5eRCQlrVq16i13L+5oXGDh\nXl5eTnV1dVAPLyKSksxsczzjdFhGRCQNKdxFRNKQwl1EJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQN\nxRXuZjbdzOrNrMHM5se4vczMnjWzNWa2zsxuSXypYXU7DvLtJ9aj5QFFRC6sw3A3sxCwGJgBVABz\nzayizbCvAA+7+2RgDvBviS70tJc37uGHz23gD7W7uuohRERSXjwz96lAg7tvdPdmYCkwu80YB/pE\nLvcFtieuxHN9+KphjLmsF/c/Vsvxk61d9TAiIiktnnAvAbZGbTdGrou2APiwmTUCy4HPJKS6GHJC\nWdw3cwJb9x7joT9u7KqHERFJafGEu8W4ru0B77nA/3P3ocAtwK/M7Lz7NrM7zazazKqbmpouvtqI\na0cVMX3CIBY/u4EdB451+n5ERNJVPOHeCJRGbQ/l/MMuHwceBnD3l4B8oKjtHbn7EnevdPfK4uIO\nm5q168u3jqfVnW8uX39J9yMiko7iCfeVwGgzG25muYQ/MK1qM2YLcAOAmY0nHO6dn5rHobR/AXdd\nN4KqV7ezctPernwoEZGU02G4u3sLMA9YAdQRPiumxswWmdmsyLC7gU+Y2avAb4A7vBvOVfzktFEM\n6ZvPfb+vofWUTo0UETktrn7u7r6c8Ael0dfdG3W5Frg2saV1rEduiHtuGc9nfrOG367cyofeUdbd\nJYiIJKWU/4bqbRMHM3V4fx5YsZ4DR08GXY6ISFJI+XA3MxbMnMCBYyf5wVOvB12OiEhSSPlwB6gY\n0ocPvaOMX728mfqdh4IuR0QkcGkR7gB33zSWXnnZLFxWo74zIpLx0ibc+/XM5e73juHPG/awomZn\n0OWIiAQqbcId4ENTyxg3qDf3P1qnvjMiktHSKtyzI31ntu0/xpIX1HdGRDJXWoU7wNUjB3Dr2wbz\nb881sG2/+s6ISGZKu3AHuOeWcbjDN5bXBV2KiEgg0jLch/Yr4JPTRvLYuh28vHFP0OWIiHS7tAx3\ngLvePZKSwh4sqKqhpfVU0OWIiHSrtA33/JwQX751POt3HuI3K7d2/AciImkkbcMdYMblg7h6xAC+\n92Q9+482B12OiEi3SetwNzPum1XBwWMn+f4f1HdGRDJHWoc7wLhBffjIVcP495c3U7fjYNDliIh0\ni7jC3cymm1m9mTWY2fwYt//AzNZGfl43s/2JL7XzvnDTGPr2yGFBlfrOiEhm6DDczSwELAZmABXA\nXDOriB7j7l9w90nuPgl4EPjvrii2swoLcrn7vWN55c29LH9NfWdEJP3FM3OfCjS4+0Z3bwaWArPb\nGT+X8FJ7SWXu1DLGD+7D1x+r5Viz+s6ISHqLJ9xLgOhzCRsj153HzIYBw4FnLr20xAplGQtnTWD7\ngeP86PkNQZcjItKl4gl3i3HdhQ5czwEecfeYU2Mzu9PMqs2suqmpKd4aE2bq8P7MvGIIP3p+A437\njnb744uIdJd4wr0RKI3aHgpsv8DYObRzSMbdl7h7pbtXFhcXx19lAt0zYxxm6jsjIuktnnBfCYw2\ns+Fmlks4wKvaDjKzsUA/4KXElphYQwp78Olpo1j+2k7+3PBW0OWIiHSJDsPd3VuAecAKoA542N1r\nzGyRmc2KGjoXWOopcK7hJ64bwdB+PVi4rFZ9Z0QkLVlQWVxZWenV1dWBPDbAE3/dyV3/voqFsybw\nsWvKA6tDRORimNkqd6/saFzaf0P1Qm6ecBnvHFXE956sZ+8R9Z0RkfSSseFuZtw3s4Ijza1878n6\noMsREUmojA13gNGX9eajVw/j13/ZQs32A0GXIyKSMBkd7gCfv3EM/QpyWVhVq74zIpI2Mj7c+/bI\n4Us3j+Uvm/aybN2OoMsREUmIjA93gA9UljJhSB++8VgdR5tbgi5HROSSKdw523dm58Hj/PA59Z0R\nkdSncI+oLO/P+yYN4ccvbGTLHvWdEZHUpnCPMn/GeLKzjK8vrw26FBGRS6JwjzKobz6fvn4UK2p2\n8cc3ur9rpYhIoijc2/j4O4dT1r+AhctqOam+MyKSohTubeTnhPjqbRU07D7Mr17aHHQ5IiKdonCP\n4cbxA3nX6CJ+8NTr7Dl8IuhyREQumsI9htN9Z441t/Jd9Z0RkRSkcL+AUQN7c8c15SxduZXXGtV3\nRkRSi8K9HZ+9cTQDeuayYFmN+s6ISEqJK9zNbLqZ1ZtZg5nNv8CYD5hZrZnVmNmvE1tmMPrk5/BP\nN49j1eZ9/H7thZaNFRFJPh2Gu5mFgMXADKACmGtmFW3GjAbuAa519wnA57ug1kDcfuVQJg7tyzcf\nr+PICfWdEZHUEM/MfSrQ4O4b3b0ZWArMbjPmE8Bid98H4O67E1tmcLKyjPtmTmDXwRMsfrYh6HJE\nROIST7iXAFujthsj10UbA4wxsz+Z2ctmNj3WHZnZnWZWbWbVTU2p8w3QK4f14/1TSnjoj2+y6a0j\nQZcjItKheMLdYlzX9tPFbGA0MA2YCzxkZoXn/ZH7EnevdPfK4uLii601UPOnjyMnZHztsbqgSxER\n6VA84d4IlEZtDwXafrrYCPze3U+6+5tAPeGwTxsD++TzmRtG81TdLp5/PXXedYhIZoon3FcCo81s\nuJnlAnOAqjZj/ge4HsDMiggfptmYyEKTwd9dW075gAIWLquhuUV9Z0QkeXUY7u7eAswDVgB1wMPu\nXmNmi8xsVmTYCmCPmdUCzwJfcvc9XVV0UPKyQ9w7s4KNTUf45Uubgi5HROSCLKgv51RWVnp1dXUg\nj32p7vj5X1i1aR/PfHEaxb3zgi5HRDKIma1y98qOxukbqp3w1dsqON7SygMr1gddiohITAr3ThhZ\n3Iu/v3Y4D1c38urW/UGXIyJyHoV7J817zyiKeuWxYFkNp06p74yIJBeFeyf1zs9h/oxxrNmyn9+t\n2RZ0OSIi51C4X4L3Ty7hitJCvvXEeg6r74yIJBGF+yXIyjIWzppA06ETPPjMG0GXIyJyhsL9Ek0q\nLeT2K4fysxffZGPT4aDLEREBFO4J8U/Tx5KXHVLfGRFJGgr3BBjYO5/P3TCaZ9bv5tn1adPtWERS\nmMI9QT52TTkjinuy6NFa9Z0RkcAp3BMkNzuLe2+r4M23jvDzP70ZdDkikuEU7gk0bexAbhg3kH95\n+g12HzwedDkiksEU7gn21dsqONnqfPuJ+qBLEZEMpnBPsPKinvz9O4fzX6sbWbNlX9DliEiGUrh3\ngXnvGcXA3nksqFLfGREJRlzhbmbTzazezBrMbH6M2+8wsyYzWxv5+YfEl5o6euVlc88t43i18QCP\nrG4MuhwRyUAdhruZhYDFwAygAphrZhUxhv7W3SdFfh5KcJ0p532TSphSVsh3nqjn4PGTQZcjIhkm\nnpn7VKDB3Te6ezOwFJjdtWWlPjNjwawJ7DlyggefVt8ZEele8YR7CbA1arsxcl1bf2Nm68zsETMr\nTUh1KW7i0EI+cGUpP//TJhp2q++MiHSfeMLdYlzX9lPCZUC5u08EngJ+EfOOzO40s2ozq25qarq4\nSlPUl6aPpUdOiEWP1hLUerUiknniCfdGIHomPhTYHj3A3fe4+4nI5k+AK2PdkbsvcfdKd68sLi7u\nTL0pp6hXHp+7cTQvvN7E03XqOyMi3SOecF8JjDaz4WaWC8wBqqIHmNngqM1ZgNojRvnYNeWMGtiL\n+x+r5URLa9DliEgG6DDc3b0FmAesIBzaD7t7jZktMrNZkWGfNbMaM3sV+CxwR1cVnIpyQlncN7OC\nzXuO8tMX1XdGRLqeBXUcuLKy0qurqwN57KB84pfV/KnhLZ794jQu65MfdDkikoLMbJW7V3Y0Tt9Q\n7UZfvbWCllPOtx5fH3QpIpLmFO7dqGxAAZ9413B+t2YbqzbvDbocEUljCvdu9qlpoxjUJ58FVbW0\nqu+MiHQRhXs36xnpO/PatgP8Z/XWjv9ARKQTFO4BmHXFEN5e3o8HVtRz4Jj6zohI4incA2Bm3Ddz\nAnuPNvPPT6nvjIgknsI9IJeX9GXO28v45UubeGPXoaDLEZE0o3AP0BffO4aC3BALl6nvjIgklsI9\nQAN65fGFm8bwYsNbPFm7K+hyRCSNKNwD9uGrhjHmsl587bFajp9U3xkRSQyFe8DCfWcmsHXvMR76\n48agyxGRNKFwTwLXjipi+oRBLH52AzsOHAu6HBFJAwr3JPHlW8dzyp1vLlffGRG5dAr3JFHav4D/\nfd0Iql7dzspN6jsjIpdG4Z5EPjltFEP65nPf72vUd0ZELonCPYn0yA1xzy3jqd1xkKUrtwRdjoik\nsLjC3cymm1m9mTWY2fx2xt1uZm5mHTaSl9humziYqcP7890V9Rw4qr4zItI5HYa7mYWAxcAMoAKY\na2YVMcb1JrzE3iuJLjKTmBkLZk7gwLGT/OCp14MuR0RSVDwz96lAg7tvdPdmYCkwO8a4+4HvAMcT\nWF9GqhjShw+9o4xfvbyZ+p3qOyMiFy+ecC8BohuPN0auO8PMJgOl7v5oe3dkZneaWbWZVTc1NV10\nsZnk7pvG0isvm4XLatR3RkQuWjzhbjGuO5M2ZpYF/AC4u6M7cvcl7l7p7pXFxcXxV5mB+vXM5e73\njuHPG/awomZn0OWISIqJJ9wbgdKo7aHA9qjt3sDlwHNmtgm4CqjSh6qX7kNTyxg3qDf3P1qnvjMi\nclHiCfeVwGgzG25mucAcoOr0je5+wN2L3L3c3cuBl4FZ7l7dJRVnkOxI35lt+4/x4+fVd0ZE4tdh\nuLt7CzAPWAHUAQ+7e42ZLTKzWV1dYKa7euQAbn3bYH74fAPb9qvvjIjEJ67z3N19ubuPcfeR7v71\nyHX3untVjLHTNGtPrHtuGQfAN5bXBVyJiKQKfUM1BQztV8Bd7x7JY+t28PLGPUGXIyIpQOGeIu56\n90hKCnuwoKqGltZTQZcjIklO4Z4i8nNCfPnW8azfeYjfrNza8R+ISEZTuKeQGZcP4uoRA/jek/Xs\nO9IcdDkiksQU7inEzLhvVgUHj53k+39Q3xkRuTCFe4oZN6gPH7lqGP/xymbqdhwMuhwRSVIK9xT0\nhZvG0LdHDguq1HdGRGJTuKegwoJc7n7vWF55cy/LX1PfGRE5n8I9Rc2dWsb4wX34+mO1HGtW3xkR\nOZfCPUWFsoyFsyaw/cBxfvj8hqDLEZEko3BPYVOH92fmFUP48fMb2Lr3aNDliEgSUbinuHtmjCPL\nTH1nROQcCvcUN6SwB5+aNpLH/7qTPze8FXQ5IpIkFO5p4BPXjaC0fw8WLqtV3xkRARTuaSE/J8SX\nb6mgftch/uOVLUGXIyJJIK5wN7PpZlZvZg1mNj/G7XeZ2WtmttbMXjSzisSXKu25ecJlvHNUEd97\nsp696jsjkvE6DHczCwGLgRlABTA3Rnj/2t3f5u6TgO8A3094pdIuM+O+mRUcaW7lu0/WB12OiAQs\nnpn7VKDB3Te6ezOwFJgdPcDdo5uc9AT0nfgAjL6sNx+9ehi/+csW/rrtQNDliEiA4gn3EiC6gXhj\n5LpzmNmnzWwD4Zn7ZxNTnlysz984hn4FuSxcpr4zIpksnnC3GNedlxruvtjdRwL/B/hKzDsyu9PM\nqs2suqmp6eIqlbj07ZHDl24ey8pN+1i2bkfQ5YhIQOIJ90agNGp7KLC9nfFLgffFusHdl7h7pbtX\nFhcXx1+lXJQPVJZyeUkfvvFYHUebW4IuR0QCEE+4rwRGm9lwM8sF5gBV0QPMbHTU5q3AG4krUS5W\nKMtYMHMCOw8e54fPqe+MSCbqMNzdvQWYB6wA6oCH3b3GzBaZ2azIsHlmVmNma4F/BD7WZRVLXCrL\n+/O+SUP48Qsb2bJHfWdEMo0F9aFbZWWlV1dXB/LYmWLngeO853vP8c5RRSz5aGXQ5YhIApjZKnfv\n8B+0vqGaxgb1zefT14/iydpd/PENfYAtkkkU7mnu4+8czrABBSxcVstJ9Z0RyRgK9zSXnxPiK7dW\n0LD7ML96aXPQ5YhIN1G4Z4Abxw/kujHF/OCp19lz+ETQ5YhIN1C4ZwAz497bKjimvjMiGUPhniFG\nDezFHdeUs3TlVl5rVN8ZkXSncM8gn71xNAN65rJAfWdE0p7CPYP0yc/hn24ex6rN+/j92vY6SIhI\nqlO4Z5jbrxzKxKF9+ebjdRw5ob4zIulK4Z5hsrKMBbMmsOvgCRY/2xB0OSLSRRTuGWhKWT/eP6WE\nh/74JpveOhJ0OSLSBRTuGWr+9HHkhIyvPVYXdCki0gUU7hlqYJ98PnPDaJ6q28Vz9buDLkdEEkzh\nnsH+7tpyygcUsOjRWppb1HdGJJ0o3DNYXnaIe2dWsLHpCL98aVPQ5YhIAsUV7mY23czqzazBzObH\nuP0fzazWzNaZ2dNmNizxpUpXeM+4y7h+bDH//NQbNB1S3xmRdNFhuJtZCFgMzAAqgLlmVtFm2Bqg\n0t0nAo8A30l0odJ1vnpbBcdbWnlgxfqgSxGRBIln5j4VaHD3je7eTHgB7NnRA9z9WXc/vZbby4QX\n0ZYUMaK4F39/7XAerm7k1a37gy5HRBIgnnAvAbZGbTdGrruQjwOPX0pR0v3mvWcURb3yWLCshlOn\n1HdGJNXFE+4W47qY//rN7MNAJfDABW6/08yqzay6qUnLviWT3vk5zJ8xjjVb9vO7NduCLkdELlE8\n4d4IlEZtDwXO6zplZjcCXwZmuXvMT+bcfYm7V7p7ZXFxcWfqlS70/sklTCot5FtPrOew+s6IpLR4\nwn0lMNrMhptZLjAHqIoeYGaTgR8TDnZ9IyZFne4703ToBA8+80bQ5YjIJegw3N29BZgHrADqgIfd\nvcbMFpnZrMiwB4BewH+a2Vozq7rA3UmSm1RayN9eOZSfvfgmG5sOB12OiHSSBbVoQ2VlpVdXVwfy\n2NK+3YeO857vPs/U4f352R1vD7ocEYliZqvcvbKjcfqGqpxnYO98PnfDaJ5Zv5tn1u8KuhwR6QSF\nu8T0sWvKGVHck/sfrVPfGZEUpHCXmHKzs7j3tgrefOsIP//Tm0GXIyIXSeEuFzRt7EBuHD+Qf3n6\nDXYfPB50OSJyERTu0q6v3FrByVbn20/UB12KiFwEhbu0q7yoJx9/13D+a3Uja7bsC7ocEYmTwl06\n9OnrRzGwdx4LqtR3RiRVKNylQ73ysrnnlnG82niAR1Y3Bl2OiMRB4S5xed+kEqaUFfKdJ9Zz8PjJ\noMsRkQ4o3CUuZuG+M3uONPPg0+o7I5LsFO4St4lDC/lgZSk//9MmGnar74xIMlO4y0X54s1j6ZEb\nYtGjtQTVl0hEOqZwl4tS1CuPz984hhdeb+LpOnV3FklWCne5aB+9ehijBvbi/sdqOdHSGnQ5IhKD\nwl0uWk4oi/tmVrB5z1F++qL6zogkI4W7dMq7RhdzU8Vl/OszDew8oL4zIskmrnA3s+lmVm9mDWY2\nP8bt15nZajNrMbPbE1+mJKOv3lpByynn20+sD7oUEWmjw3A3sxCwGJgBVABzzayizbAtwB3ArxNd\noCSvsgEF3PmuEfxuzTZWbd4bdDkiEiWemftUoMHdN7p7M7AUmB09wN03ufs6QKs6ZJhPXT+SQX3y\nWVBVS6v6zogkjew4xpQAW6O2G4F3dE05kmoKcsN9Zz63dC3vfuBZ8nNCZGcZoSyL+p1Fdij6uqzw\n79DZMTlZWedsnzMusp0ditqOcX/Zodh/d3Z8POPOv87Mgv7PLHLR4gn3WP9nd2qKZmZ3AncClJWV\ndeYuJAnNumIIjfuOUb/zEK2nnJZTpyK/Pfy7Nfxz7GTrme1Y4062Oq2nTp39u8jvoN8RhNqEfU4o\nq82L19kXsfDt526f+0LUzovd6e1Q5MUuarvDF7sYL2KhLKOwIIdRxb3IDunciUwTT7g3AqVR20OB\n7Z15MHdfAiwBqKys1Hv4NGFmfPr6UV12/+5nw77llNPaGuMF5NTZF4aW1qjtVo857mTr+X937guP\n09La9oXm3O3zbm+NfX8nWlpj1BDZh1bn5Jmx5+9TIhTkhpg4tC9Tyvoxuawfk8sKKeqVl5D7luQV\nT7ivBEab2XBgGzAH+FCXViUSxSxyKCUUdCXdy9055dASeZE650WhzQtR+J1P2xeHU+w+eIK1W/ez\ness+lryw8cwLRln/AqaUFTK5rB9TyvoxbnBvcjS7Tysdhru7t5jZPGAFEAJ+5u41ZrYIqHb3KjN7\nO/A7oB8w08wWuvuELq1cJM2ZGSGDUFaIvHimYRfwvsklABw/2cpr2w6wZss+Vm/ez5837OF/1obf\nhOfnZDGxpJDJZwK/kIF98hOxGxIQC6r5U2VlpVdXVwfy2CISfmew/cDxM2G/Zus+arYdpLk1fNJb\nSWEPJpcVRg7nFDJhSF9yszW7D5qZrXL3yo7GXcJ8QERSmZlRUtiDksIe3DZxCBCe3dfuOMjqzftY\ns3U/qzfv49F1OwDIzc7i8iF9zhy7nzKskMF9ewS5C9IOzdxFpF07I7P702G/btsBmlvCs/tBffKZ\nMqyQyaXhsJ8wpC/5ORn24Ug308xdRBJiUN98ZrxtMDPeNhiA5pZT1O04yOot+1izJfxh7fLXdgKQ\nEzIqhvRlcmkhU4aFj92XFPbQdwUCoJm7iFyy3YeOs2bL/jNhv65xP8dPhmf3xb3zzjkz520lfemR\nq9l9Z2nmLiLdZmDvfG6eMIibJwwC4GTrKep3Hgp/WLtlP2u27GNFzS4AsrOM8YP7nPNhbVn/As3u\nE0wzdxHpFnsOnz3nfs2W/by6dT9HmsOLvQzomXvmNMzJZYVcMbSQnpdy/mca08xdRJLKgF553DD+\nMm4YfxkAraec13cdOufY/VORpRuzDMYO6hN1OKeQ4UU9Nbu/CJq5i0jS2H+0mTVb90eO3+9j7Zb9\nHDrRAkBhQQ6TS88eu7+itC+983MCrrj7aeYuIimnsCCX68cO5PqxAwE4dcppaDocPu9+S/iLVs+9\n3oQ7mMGYgb3POXY/srgXWVma3YNm7iKSYg4cO8mrW8+embNmyz4OHg/P7nvnZzOp9GzYTy7tR9+C\n9Jrda+YuImmpb48crhtTzHVjioHw7H7jW0fOOTPnwWfe4HRTzZHFPZlS1o8pw8KBP3pgb0IZMLvX\nzF1E0s7hEy2sizozZ/WWfew7ehKAXnnZXFHa95zZfb+euQFXHD/N3EUkY/XKy+aaUUVcM6oICDdJ\n27zn6Dlh/2/PbTizEMzwop7ndMQce1nvlF/gRDN3EclIR5tbWNd44Jxj928dbgbOLnBy+sycZFrg\nRDN3EZF2FORmc9WIAVw1YgAQnt037jt2zuz+J20WOIk+M2f84D5JvcBJXOFuZtOBfya8WMdD7v6t\nNrfnAb8ErgT2AB90902JLVVEpOuYGaX9CyjtX8DsSbEXOHlpwx5+H1ngJC87K2r5wnDoJ9MCJx0e\nljGzEPA6cBPh9VRXAnPdvTZqzKeAie5+l5nNAf6Xu3+wvfvVYRkRSTXxLnBy+th9xZA+5CV4fchE\nHpaZCjS4+8bIHS8FZgO1UWNmAwsilx8B/tXMzIM6oC8i0gViLXByoqWVmu3tL3ASfex+SGH3LHAS\nT7iXAFujthuBd1xoTGTN1QPAAOCtRBQpIpKs8rJD4fPoy/qdua7tAif//vJmfvrim0B4gZN7bhl3\n5tBPV4kn3GOd7d92Rh7PGMzsTuBOgLKysjgeWkQk9VxogZPTX7Qa2Lvrj83HE+6NQGnU9lBg+wXG\nNJpZNtAX2Nv2jtx9CbAEwsfcO1OwiEiqyc3O4orSQq4oLeSOa7vnMeM5j2clMNrMhptZLjAHqGoz\npgr4WOTy7cAzOt4uIhKcDmfukWPo84AVhE+F/Jm715jZIqDa3auAnwK/MrMGwjP2OV1ZtIiItC+u\n89zdfTmwvM1190ZdPg78bWJLExGRzkrer1eJiEinKdxFRNKQwl1EJA0p3EVE0pDCXUQkDQXWz93M\nmoDNnfzzItKntYH2Jfmky36A9iVZXcq+DHP34o4GBRbul8LMquPpipYKtC/JJ132A7Qvyao79kWH\nZURE0pDCXUQkDaVquC8JuoAE0r4kn3TZD9C+JKsu35eUPOYuIiLtS9WZu4iItCOpw93MpptZvZk1\nmNn8GLfnmdlvI7e/Ymbl3V9lfOLYlzvMrMnM1kZ+/iGIOjtiZj8zs91m9tcL3G5m9i+R/VxnZlO6\nu8Z4xbEv08zsQNRzcm+scUEzs1Ize9bM6sysxsw+F2NMSjwvce5Lqjwv+Wb2FzN7NbIvC2OM6boM\nc/ek/CHcXngDMALIBV4FKtqM+RTwo8jlOcBvg677EvblDuBfg641jn25DpgC/PUCt98CPE54da6r\ngFeCrvkS9mUa8GjQdcaxH4OBKZHLvQkvaN/2/6+UeF7i3JdUeV4M6BW5nAO8AlzVZkyXZVgyz9zP\nLMzt7s3A6YW5o80GfhG5/Ahwg5nFWvIvaPHsS0pw9xeIscpWlNnALz3sZaDQzAZ3T3UXJ459SQnu\nvsPdV0cuHwLqCK9rHC0lnpc49yUlRP5bH45s5kR+2n7I2WUZlszhHmth7rZP8jkLcwOnF+ZONvHs\nC8DfRN4yP2JmpTFuTwXx7muquDrytvpxM5sQdDEdibytn0x4lhgt5Z6XdvYFUuR5MbOQma0FdgN/\ncPcLPi+JzrBkDveELcydBOKpcxlQ7u4Tgac4+2qealLlOYnHasJf9b4CeBD4n4DraZeZ9QL+C/i8\nux9se3MoiHd7AAABdElEQVSMP0na56WDfUmZ58XdW919EuG1p6ea2eVthnTZ85LM4X4xC3PT3sLc\nSaDDfXH3Pe5+IrL5E+DKbqot0eJ53lKCux88/bbaw6uR5ZhZUcBlxWRmOYTD8D/c/b9jDEmZ56Wj\nfUml5+U0d98PPAdMb3NTl2VYMod7Oi3M3eG+tDn+OYvwscZUVAV8NHJ2xlXAAXffEXRRnWFmg04f\n/zSzqYT/vewJtqrzRWr8KVDn7t+/wLCUeF7i2ZcUel6KzawwcrkHcCOwvs2wLsuwuNZQDYKn0cLc\nce7LZ81sFtBCeF/uCKzgdpjZbwifrVBkZo3AfYQ/KMLdf0R4rd1bgAbgKPB3wVTasTj25Xbgk2bW\nAhwD5iTp5OFa4CPAa5HjuwD/FyiDlHte4tmXVHleBgO/MLMQ4Regh9390e7KMH1DVUQkDSXzYRkR\nEekkhbuISBpSuIuIpCGFu4hIGlK4i4ikIYW7iEgaUriLiKQhhbuISBr6//y3IQNQJH5XAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f36fc60bc90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now, we will print how much \"variance\" each of the principal components explains. Roughly speaking, it tells us how much % of the spread lies\n",
    "# along each direction\n",
    "plt.plot(pca.explained_variance_ratio_)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see above, the first component explains 80% of the total data spread. The second and third component explains roughly 10% each, \n",
    "and the last component explains nothing! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA on real dataset\n",
    "\n",
    "The MNIST [dataset](http://yann.lecun.com/exdb/mnist/) is a famous set of photographs, each showing one handwritten digit. We will work with a processed\n",
    "version of the data set. You can find the dataset in my github [repository](https://github.com/govg/acass/tree/master/data). Please download the different\n",
    "datasets. The folder for each data set will also contain python files that show how to convert from raw data -> python array. Please take a look at those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's work with some real data!\n",
    "import os, sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11022, 784)\n"
     ]
    }
   ],
   "source": [
    "# The MNIST dataset defines images of the size 28x28. We treat the entire image as a row!\n",
    "X = np.load(\"../data/mnist/Xtrain.npy\")\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,\n",
       "  svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Since we assume images to be very correlated, we can now do PCA and see the effect of it on reducing dimensions!\n",
    "# Note that this may time some time depending on how powerful your machine is. On my machine, doing PCA on the entire\n",
    "# MNIST dataset that was loaded took a few seconds.\n",
    "pca2 = PCA()\n",
    "pca2.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGtJJREFUeJzt3X2UXPV93/H3Z2a0K+0K9LgEWQ+WxNOxOHCArAVuUuyY\nmginQU4rYsk5x7glVROHk7pum8Jpi22ak1NSH+Om4dTWCaQYah5M41bFchRsTHrq2LKWZ4QsvAiM\n1gJLQg9YD7Da3W//uHfFMHtn9kra3RnufF7n7Nk7v/ubud/dmf3Mb39zHxQRmJlZeyg1uwAzM5s6\nDn0zszbi0DczayMOfTOzNuLQNzNrIw59M7M24tA3M2sjDn0zszbi0DczayOVZhdQa/78+bF06dJm\nl2Fm9q7y+OOP74uInvH6tVzoL126lL6+vmaXYWb2riLpp3n6eXrHzKyNOPTNzNqIQ9/MrI049M3M\n2ohD38ysjTj0zczaiEPfzKyNFCb0j7w1xJf+ZgdP7TrY7FLMzFpWYUL/zePD/Nmj/Twz4NA3M6un\nMKFfkgAYHvGF3s3M6ilO6JeS0Hfmm5nVV5zQTzKfEae+mVldhQn98omRvkPfzKyewoT+6Jy+B/pm\nZvUVJvQ1Or3jkb6ZWV2FCf3y6EjfQ30zs7pyhb6kVZJ2SOqXdFPG+islPSFpSNKaqvZLJP1A0jZJ\nz0j6+EQWX+3ELpse6ZuZ1TVu6EsqA3cA1wArgHWSVtR0ewX4FPD1mvajwCcj4kJgFfBlSbNPt+gs\n3mXTzGx8eS6XuBLoj4idAJLuB1YDz492iIiX03Uj1XeMiBeqlndL2gP0AJNy2GxJnt4xM2skz/TO\nQmBX1e2BtO2kSFoJdAAvZqxbL6lPUt/evXtP9qFPKJfkD3LNzBrIE/rKaDupZJW0ALgH+CcRMVK7\nPiI2RERvRPT29Ix7MfdG2/GcvplZA3lCfwBYXHV7EbA77wYknQl8C/j3EfHDkyvv5JQlnPlmZvXl\nCf2twHmSlknqANYCG/M8eNr/m8DXIuIbp15mPiX5hGtmZo2MG/oRMQTcCGwGtgMPRsQ2SbdKuhZA\n0vslDQDXAV+VtC29+28DVwKfkvRU+nXJpPwkJHvweE7fzKy+PHvvEBGbgE01bbdULW8lmfapvd+9\nwL2nWWNuJU/vmJk1VJgjcsHTO2Zm4ylU6HuXTTOzxgoV+pJD38yskUKFflliZMxRAGZmNqpQoV+S\nT7hmZtZIsULfc/pmZg0VK/Qln3DNzKyBQoV+svdOs6swM2tdhQp9eU7fzKyhQoV+csI1h76ZWT2F\nCv2Sd9k0M2uoUKHv6R0zs8YKFfrlkqd3zMwaKVTolySfcM3MrIFihb532TQza6hYoS98RK6ZWQOF\nCv2yz7JpZtZQoULfc/pmZo0VK/RLeE7fzKyBYoW+T7hmZtZQoULfl0s0M2usUKEviWFnvplZXblC\nX9IqSTsk9Uu6KWP9lZKekDQkaU3Nuusl/ST9un6iCs9SEj4i18ysgXFDX1IZuAO4BlgBrJO0oqbb\nK8CngK/X3Hcu8DngcmAl8DlJc06/7GzeZdPMrLE8I/2VQH9E7IyIQeB+YHV1h4h4OSKeAWrPcfnr\nwCMRsT8iDgCPAKsmoO5Mkhj2WTbNzOrKE/oLgV1VtwfStjxy3VfSekl9kvr27t2b86HHKpfw3jtm\nZg3kCX1ltOVN1lz3jYgNEdEbEb09PT05H3qsckk+tbKZWQN5Qn8AWFx1exGwO+fjn859T1rJc/pm\nZg3lCf2twHmSlknqANYCG3M+/mbgaklz0g9wr07bJoUPzjIza2zc0I+IIeBGkrDeDjwYEdsk3Srp\nWgBJ75c0AFwHfFXStvS++4H/SPLGsRW4NW2bFJ7eMTNrrJKnU0RsAjbVtN1StbyVZOom6753AXed\nRo25+Rq5ZmaNFeqI3HLJ59M3M2ukUKHvUyubmTVWrND3CdfMzBoqVOgnp2FodhVmZq2rUKFfEp7e\nMTNroFihX/J++mZmjRQq9MvyfvpmZo0UK/T9Qa6ZWUOFCn354Cwzs4YKFfrlEp7eMTNroFih74Oz\nzMwaKlTol0rJ6ft9nVwzs2zFCn0loe/RvplZtkKFfjkd6Xte38wsW6FCf3Sk7z14zMyyFSr0y+lP\n4331zcyyFSr0T8zpO/TNzDIVMvR9/h0zs2yFCv0TH+Q69M3MMhUq9NPM9zn1zczqKFbop6nvD3LN\nzLIVKvTLPjjLzKyhXKEvaZWkHZL6Jd2Usb5T0gPp+i2Slqbt0yTdLelZSdsl3Tyx5b+TR/pmZo2N\nG/qSysAdwDXACmCdpBU13W4ADkTEucDtwG1p+3VAZ0RcBPwy8M9H3xAmgw/OMjNrLM9IfyXQHxE7\nI2IQuB9YXdNnNXB3uvwQcJUkAQF0S6oAM4BB4I0JqTzD6MFZ3k/fzCxbntBfCOyquj2QtmX2iYgh\n4BAwj+QN4AjwKvAK8MWI2H+aNdflE66ZmTWWJ/SV0VabqvX6rASGgfcAy4B/JWn5mA1I6yX1Serb\nu3dvjpKylX1qZTOzhvKE/gCwuOr2ImB3vT7pVM4sYD/wCeCvI+J4ROwBvg/01m4gIjZERG9E9Pb0\n9Jz8T5HyaRjMzBrLE/pbgfMkLZPUAawFNtb02Qhcny6vAR6NZLj9CvBhJbqBK4AfT0zpY3l6x8ys\nsXFDP52jvxHYDGwHHoyIbZJulXRt2u1OYJ6kfuCzwOhunXcAM4HnSN48/jIinpngn+GE0ekd771j\nZpatkqdTRGwCNtW03VK1/CbJ7pm19zuc1T5ZKr6IiplZQ4U6Ird04oRrHuqbmWUpVOifGOk7883M\nMhUq9Ec/yB3ySN/MLFOhQr9S9ge5ZmaNFCr0PdI3M2usUKFf8Vk2zcwaKlToj+6nPzTs0Dczy1LI\n0PdI38wsWyFDf8inYTAzy1TI0Pe5d8zMshUr9H3CNTOzhooV+p7eMTNrqFCh//bBWQ59M7MshQr9\nsjzSNzNrpFih7102zcwaKmTo++AsM7NshQx9j/TNzLIVMvQ9p29mlq2Qoe/99M3MshUr9H1wlplZ\nQ8UKfY/0zcwaKlToS6Ikh76ZWT25Ql/SKkk7JPVLuiljfaekB9L1WyQtrVp3saQfSNom6VlJ0yeu\n/LEqpRLD3nvHzCzTuKEvqQzcAVwDrADWSVpR0+0G4EBEnAvcDtyW3rcC3Av8XkRcCHwIOD5h1Wco\nlTzSNzOrJ89IfyXQHxE7I2IQuB9YXdNnNXB3uvwQcJUkAVcDz0TE0wAR8XpEDE9M6dkqpZJD38ys\njjyhvxDYVXV7IG3L7BMRQ8AhYB5wPhCSNkt6QtIfnX7JjXlO38ysvkqOPspoq03Ven0qwK8C7weO\nAt+V9HhEfPcdd5bWA+sBlixZkqOk+irlEkMjI6f1GGZmRZVnpD8ALK66vQjYXa9POo8/C9iftv9t\nROyLiKPAJuCy2g1ExIaI6I2I3p6enpP/KaqUJIad+WZmmfKE/lbgPEnLJHUAa4GNNX02Ateny2uA\nRyMigM3AxZK60jeDDwLPT0zp2SolMeyRvplZpnGndyJiSNKNJAFeBu6KiG2SbgX6ImIjcCdwj6R+\nkhH+2vS+ByR9ieSNI4BNEfGtSfpZgOQALY/0zcyy5ZnTJyI2kUzNVLfdUrX8JnBdnfveS7Lb5pQo\ne6RvZlZXoY7IhXR6xzvvmJllKlzolzzSNzOrq3Chn3yQ66G+mVmWwoV+ssumQ9/MLEvhQr9Sduib\nmdVTuNAvSb5coplZHYUL/UpJvjC6mVkdhQv9UkkMeZ9NM7NMhQt9j/TNzOorXOiXS57TNzOrp5Ch\nP+LQNzPLVLzQ9947ZmZ1FS/0fUSumVldDn0zszbi0DczayPFDH3vsmlmlqmQoe+Ds8zMshUu9H1w\nlplZfYULfR+cZWZWXyFD3wdnmZllK17o++AsM7O6ihf6pZJH+mZmdeQKfUmrJO2Q1C/ppoz1nZIe\nSNdvkbS0Zv0SSYcl/euJKbu+cgmP9M3M6hg39CWVgTuAa4AVwDpJK2q63QAciIhzgduB22rW3w58\n+/TLHV+5VGJoZGQqNmVm9q6TZ6S/EuiPiJ0RMQjcD6yu6bMauDtdfgi4SpIAJH0M2Alsm5iSG+vu\nKHN8OBgccvCbmdXKE/oLgV1VtwfStsw+ETEEHALmSeoG/i3whdMvNZ+Z0ysAHHlraKo2aWb2rpEn\n9JXRVjtpXq/PF4DbI+Jwww1I6yX1Serbu3dvjpLqm9mZhP5hh76Z2RiVHH0GgMVVtxcBu+v0GZBU\nAWYB+4HLgTWS/hSYDYxIejMi/rz6zhGxAdgA0Nvbe1qfwp6RjvR/8aZD38ysVp7Q3wqcJ2kZ8DNg\nLfCJmj4bgeuBHwBrgEcjIoC/P9pB0ueBw7WBP9Fmdk4DPNI3M8sybuhHxJCkG4HNQBm4KyK2SboV\n6IuIjcCdwD2S+klG+Gsns+hGRuf0D791vFklmJm1rDwjfSJiE7Cppu2WquU3gevGeYzPn0J9J210\nTt/TO2ZmYxXuiNw5Xcn0zoEjg02uxMys9RQu9Gd3dSDB/qOe3jEzq1W40C+XxOwZ09h/5K1ml2Jm\n1nIKF/oAc7s72O/pHTOzMQoZ+vO6O9n7C4/0zcxqFTL0l87v4qV9R5tdhplZyylk6J/TM5N9h9/i\nkD/MNTN7h0KG/qI5XQC8+saxJldiZtZaChn6s9N99Q96pG9m9g6FDP1ZMxz6ZmZZChn6oyP9N445\n9M3MqhU09DsAOHjM++qbmVUrZOh3d5Qpl+TpHTOzGoUMfUl0d5Q5Ojjc7FLMzFpKIUMfoLuz4uvk\nmpnVKGzod3mkb2Y2RmFDf2ZnhSODHumbmVUrbOh3dVQ4+pZH+mZm1Qob+t2dZY/0zcxqFDb0uzoq\nntM3M6tR2NDv7ixz6NhxIqLZpZiZtYzChv6lS+aw/8ggT+462OxSzMxaRq7Ql7RK0g5J/ZJuyljf\nKemBdP0WSUvT9o9IelzSs+n3D09s+fVdtmQ2AAMHfHplM7NR44a+pDJwB3ANsAJYJ2lFTbcbgAMR\ncS5wO3Bb2r4P+M2IuAi4Hrhnogofz8zO5KRrh9/0h7lmZqPyjPRXAv0RsTMiBoH7gdU1fVYDd6fL\nDwFXSVJEPBkRu9P2bcB0SZ0TUfh4Zk6vAHD4LZ9/x8xsVJ7QXwjsqro9kLZl9omIIeAQMK+mzz8G\nnoyIKbliede0MpJH+mZm1So5+iijrXaXmIZ9JF1IMuVzdeYGpPXAeoAlS5bkKGl8pZKY2VHhFz7/\njpnZCXlG+gPA4qrbi4Dd9fpIqgCzgP3p7UXAN4FPRsSLWRuIiA0R0RsRvT09PSf3EzQwc7pPumZm\nVi1P6G8FzpO0TFIHsBbYWNNnI8kHtQBrgEcjIiTNBr4F3BwR35+oovOa2VnhjWMOfTOzUeOGfjpH\nfyOwGdgOPBgR2yTdKunatNudwDxJ/cBngdHdOm8EzgX+g6Sn0q+zJvynqOP8XzqDvp8eYHjEB2iZ\nmUG+OX0iYhOwqabtlqrlN4HrMu73x8Afn2aNp+yD5/fwrWdf5ZX9R1k2v7tZZZiZtYzCHpELcM5Z\nMwHYufdwkysxM2sNxQ79nmR0v/3VN5pciZlZayh06M/u6uCSxbP5622vNbsUM7OWUOjQB7h8+Vxe\neO0wQ8MjzS7FzKzpCh/658yfyeDwiE+8ZmZGO4T+Wcm8/s59/jDXzKzwob98frIHz4t7jjS5EjOz\n5it86M/p7mBudwePvbDHV9Eys7ZX+NAH+J3Ll/D9/td5aZ9H+2bW3toi9H+7Nzlf3GM79ja5EjOz\n5mqL0F88t4vlPd089oJD38zaW1uEPsCHzj+LH+583adaNrO21jah/9GLzmZwaIQvf+eFZpdiZtY0\nbRP6vUvn8o8uW8id/+8lHnn+580ux8ysKdom9AE+95sXcsHZZ3Lrw9t8jn0za0ttFfqzZkzjxl87\nl137j/GZB55qdjlmZlOurUIf4DcuXsCnP3QO/+fp3Wz4v5mX7DUzK6y2C32A3//QOaxcNpc/2fRj\nNj5de413M7PiasvQP2P6NL72T1fS+945/OF9T/L5jds4OuhdOc2s+Noy9AGmTytz7+9ezicuX8J/\n/7uX+bUvPsY3+nZxbHC42aWZmU2atg19SIL/T37rIu77Z1dw9pnT+TcPPcMH//P3uO9Hr/ggLjMr\nJLXamSd7e3ujr69vyrc7MhJ8+7nX+MrfvsizPztER6XElefN54rl87hi+Tzet+BMyiVNeV1mZnlI\nejwiesfrV8n5YKuA/wKUgb+IiP9Us74T+Brwy8DrwMcj4uV03c3ADcAw8IcRsfkkfo4pUyqJ37h4\nAR+96Gz+7sXXeeT5n/O9HXv4zvY9AJw5vcKlS+Zw4XvO5IKzz2DJ3C4Wzeli/swOJL8ZmNm7w7gj\nfUll4AXgI8AAsBVYFxHPV/X5NHBxRPyepLXAb0XExyWtAO4DVgLvAb4DnB8RdSfOmzXSr+fVQ8fY\nsnM/P9z5Ok/tOkj/nsMMVR3Y1VkpsXDODBbOnsHc7g7mdCXn75/T3cHcrg5mzZjGjI4yXenXjI4y\nM6aV6eqo+D8HM5swEznSXwn0R8TO9IHvB1YDz1f1WQ18Pl1+CPhzJcPf1cD9EfEW8JKk/vTxfpD3\nB2m2BbNm8LFLF/KxSxcC8NbQMC/vO8qu/UcZOHCUnx08xsCBY+w+eIyfvn6UA0cG+UXOzwM6KqX0\nDaD6zaBMZ6VMpSwqpRKVkqiUxbRyiXJJTBttLytdV2JaSZTTtneuz+qXPFapJEoCkXxHUJIQnFhH\nuk56u68EGu2bfs/Td7S/qrb5jr4lkm2P19f/VZmdljyhvxDYVXV7ALi8Xp+IGJJ0CJiXtv+w5r4L\nT7naFtBZKXPB2Wdwwdln1O0zODTCwaODHDh6nINHBzl2fJhjg8McOz7M0cHa5aEx7ceODzP05gjH\nh4PhkeD4yAhDw8HQ8AhDI8HQSHB8OGkbXd9iH800Ve37Qu3bRNYbx9g+tesz3mzG3U7WXdSwT9Zb\nWm29Y/qcymNMwHaz339PbjuT9Tsau378tql4brKaqh/jfQvO5L+uuzRjSxMnT+hn/ay1EVOvT577\nImk9sB5gyZIlOUpqbR2VEmedOZ2zzpw+ZdscTt8IhkeCoeGqN4rq76PrhkcYiSACRgIiIvlO0pa0\nR9KHZH1tX0i+v/04ydNa+7iRPm5130gfd2Qk/T5e3/TxRjLe2ca01PSpXZ/15hg1vWr7ZL2fju0z\n/p3G1tJ4u9n3abzdPG/+Y7ab2efktzP+7zrHY0zCdsc8Nxl3Gve5GfsI475OsqbOx/mVsHjOjIwt\nTaw8oT8ALK66vQioPYx1tM+ApAowC9if875ExAZgAyRz+nmLt7eVS6JcKje7DDNrcXn2098KnCdp\nmaQOYC2wsabPRuD6dHkN8Ggkb3MbgbWSOiUtA84DfjQxpZuZ2ckad6SfztHfCGwm2WXzrojYJulW\noC8iNgJ3AvekH9TuJ3ljIO33IMmHvkPAHzTac8fMzCaXD84yMyuAvLtstvVpGMzM2o1D38ysjTj0\nzczaiEPfzKyNOPTNzNpIy+29I2kv8NPTeIj5wL4JKmeitWptrVoXuLZT5dpOzbu5tvdGRM94D9Jy\noX+6JPXl2W2pGVq1tlatC1zbqXJtp6YdavP0jplZG3Hom5m1kSKG/oZmF9BAq9bWqnWBaztVru3U\nFL62ws3pm5lZfUUc6ZuZWR2FCX1JqyTtkNQv6aYmbP8uSXskPVfVNlfSI5J+kn6fk7ZL0p+ltT4j\n6bJJrm2xpO9J2i5pm6R/0Sr1SZou6UeSnk5r+0LavkzSlrS2B9LTepOepvuBtLYtkpZOVm3p9sqS\nnpT0cIvV9bKkZyU9JakvbWv685lub7akhyT9OH3NfaAVapN0Qfr7Gv16Q9JnWqG2dHv/Mv0beE7S\nfenfxsS/3pIrE727v0hO+fwisBzoAJ4GVkxxDVcClwHPVbX9KXBTunwTcFu6/FHg2yRXFrsC2DLJ\ntS0ALkuXzyC50P2KVqgv3cbMdHkasCXd5oPA2rT9K8Dvp8ufBr6SLq8FHpjk391nga8DD6e3W6Wu\nl4H5NW1Nfz7T7d0N/G663AHMbpXaqmosA68B722F2kguI/sSMKPqdfapyXi9Tfovd4qewA8Am6tu\n3wzc3IQ6lvLO0N8BLEiXFwA70uWvAuuy+k1Rnf8b+Eir1Qd0AU+QXIN5H1CpfX5JruvwgXS5kvbT\nJNWzCPgu8GHg4fSPv+l1pdt4mbGh3/TnEzgzDS+1Wm019VwNfL9VauPt64zPTV8/DwO/Phmvt6JM\n72RdvL0VLsD+SxHxKkD6/ay0vWn1pv8GXkoyom6J+tIplKeAPcAjJP+1HYyIoYztn6gtXX8ImDdJ\npX0Z+CNgJL09r0XqguTqqn8j6XEl15iG1ng+lwN7gb9Mp8X+QlJ3i9RWbS1wX7rc9Noi4mfAF4FX\ngFdJXj+PMwmvt6KEfq4LsLeQptQraSbwP4HPRMQbjbpmtE1afRExHBGXkIysVwLva7D9KalN0j8E\n9kTE49XNza6ryq9ExGXANcAfSLqyQd+prK1CMs353yLiUuAIyZRJPVP+t5DOi18LfGO8rhltk1Jb\n+jnCamAZ8B6gm+S5rbf9U66tKKGf6wLsTfBzSQsA0u970vYpr1fSNJLA/x8R8VetVh9ARBwEHiOZ\nP50tafRyntXbP1Fbun4WySU6J9qvANdKehm4n2SK58stUBcAEbE7/b4H+CbJm2UrPJ8DwEBEbElv\nP0TyJtAKtY26BngiIn6e3m6F2v4B8FJE7I2I48BfAX+PSXi9FSX081y8vRmqLxh/Pclc+mj7J9O9\nA64ADo3+ezkZJInkOsbbI+JLrVSfpB5Js9PlGSQv/u3A94A1dWobrXkN8GikE5sTKSJujohFEbGU\n5PX0aET8TrPrApDULemM0WWS+ennaIHnMyJeA3ZJuiBtuorkGtlNr63KOt6e2hmtodm1vQJcIakr\n/Xsd/b1N/Ottsj8wmaovkk/aXyCZD/53Tdj+fSRzccdJ3oVvIJlj+y7wk/T73LSvgDvSWp8Feie5\ntl8l+dfvGeCp9OujrVAfcDHwZFrbc8Atafty4EdAP8m/4Z1p+/T0dn+6fvkUPLcf4u29d5peV1rD\n0+nXttHXeys8n+n2LgH60uf0fwFzWqi2LuB1YFZVW6vU9gXgx+nfwT1A52S83nxErplZGynK9I6Z\nmeXg0DczayMOfTOzNuLQNzNrIw59M7M24tA3M2sjDn0zszbi0DczayP/H0FwdbAcIKwGAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f36fc2d2490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(pca2.explained_variance_ratio_)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is amazing! It tells us that of 784 possible diretions, only the first 100 or so are really important. Let's try out a very toy classification problem then, \n",
    "to see if this is actually the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11022,)\n"
     ]
    }
   ],
   "source": [
    "Y = np.load(\"../data/mnist/Ytrain.npy\")\n",
    "print(Y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xtr = pca2.transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11022, 99)\n"
     ]
    }
   ],
   "source": [
    "# Choose only the first 100 features\n",
    "Xtr = Xtr[:,0:100]\n",
    "print(Xtr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2010, 784)\n",
      "(2010,)\n"
     ]
    }
   ],
   "source": [
    "# Now let's load the test data!\n",
    "Xts = np.load(\"../data/mnist/Xtest.npy\")\n",
    "Yts = np.load(\"../data/mnist/Ytest.npy\")\n",
    "print(Xts.shape)\n",
    "print(Yts.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now compare the effect of PCA, by training a classifier on both our original data, as well as the transformed data. If the latter gives us as good \n",
    "performance as the former, then doing PCA was actually useful! Let's use a Random Forest classifier here, of course, you can use whatever you feel like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# The parameters are : 20 trees, each tree going up to 20 depth, and allow training of 3 trees in parallel. If you have a less powerful machine, you might want to reduce\n",
    "# the n_jobs parameter to 1 or 2. \n",
    "clf1 = RandomForestClassifier(n_estimators = 30, max_depth = 10, n_jobs = 3)\n",
    "clf2 = RandomForestClassifier(n_estimators = 30, max_depth = 10, n_jobs = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',\n",
       "            max_depth=10, max_features='auto', max_leaf_nodes=None,\n",
       "            min_impurity_split=1e-07, min_samples_leaf=1,\n",
       "            min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
       "            n_estimators=30, n_jobs=3, oob_score=False, random_state=None,\n",
       "            verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train the first classifier on the original data\n",
    "clf1.fit(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',\n",
       "            max_depth=10, max_features='auto', max_leaf_nodes=None,\n",
       "            min_impurity_split=1e-07, min_samples_leaf=1,\n",
       "            min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
       "            n_estimators=30, n_jobs=3, oob_score=False, random_state=None,\n",
       "            verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train the second classifier on the transformed data, only the first 100 directions\n",
    "clf2.fit(Xtr, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.993532338308\n"
     ]
    }
   ],
   "source": [
    "# Find the accuracy of the first classifier on the test data\n",
    "print((Yts == clf1.predict(Xts)).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.946766169154\n"
     ]
    }
   ],
   "source": [
    "# Find the accuracy of the second classifier on the test data. Notice that we need to transform the testing data according to our PCA model before feeding it\n",
    "# into the classifier\n",
    "print((Yts == clf2.predict(pca2.transform(Xts))).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alas! Doing PCA did not help in this situation. It could be because of multiple reasons. One simple reason could be that our Random Forest classifier does not \n",
    "try to exploit the geometry of the problem by finding seperating lines. In this case, if the features were correlated and had some sort of linear seperation, doing\n",
    "PCA actually destroys that geometry somewhat. This should tell us that using PCA blindly can also be an issue. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python2",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}