12_principal_components_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook deals with Principal Component Analysis (PCA), which is a method for finding a new (and hopefully more efficient) representation of continuously valued features.  PCA is used for two primary tasks.  The first is the analysis of combinations of independent variables in an effort to describe the patterns underlying features in a more concise way.  The second and more ubiquitous use is for dimensionality reduction: throwing away data that is redundant or non-explanatory.  Let's look at an example for what we mean by this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import statements\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.figsize'] = (10,10)\n",
    "\n",
    "# Let's generate a synthetic dataset, with two features, and one label\n",
    "x_1 = np.random.randn(21) + 1.0\n",
    "x_2 = x_1.copy()\n",
    "X = np.column_stack((x_1,x_2))\n",
    "y = (x_1>1) + (x_2>1)\n",
    "\n",
    "plt.scatter(X[:,0],X[:,1],c=y,cmap=plt.cm.prism)\n",
    "plt.xlabel('x_1')\n",
    "plt.ylabel('x_2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ignoring the class labels for a moment, what is the dimensionality of the above dataset?  Ostensibly, it is two dimensional because the data varies in both $x_1$ and $x_2$.  However, what would happen if we simply rotated the coordinate system by 45 degrees (either clockwise or counter-clockwise)?  This rotation matrix is given by \n",
    "$$\n",
    "R = \\begin{bmatrix} \\cos \\theta & -\\sin \\theta \\\\\n",
    "                               \\sin \\theta & \\cos\\theta \\end{bmatrix}\n",
    "                              $$\n",
    "If we multiply our data matrix $X$ by this matrix, we get the following new data matrix $X_{new}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.deg2rad(-45)\n",
    "R = np.array([[np.cos(theta),-np.sin(theta)],\n",
    "              [np.sin(theta),np.cos(theta)]])\n",
    "\n",
    "X_new = np.dot(X,R.T)\n",
    "print(X_new)\n",
    "\n",
    "plt.scatter(X_new[:,0],X_new[:,1],c=y,cmap=plt.cm.prism)\n",
    "plt.xlabel('x\\'_1')\n",
    "plt.ylabel('x\\'_2')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, for the right choice of coordinates, this dataset becomes one dimensional.  What was the right choice of coordinates?  In this case, our new axes (expressed in terms of the original ones) were \n",
    "$$\n",
    "x'_1 = a x_1 + b x_2\n",
    "$$\n",
    "$$\n",
    "x'2 =  c x_2 - d x_1,\n",
    "$$\n",
    "with $a,b,c,d=\\frac{1}{\\sqrt{2}}$.  These new axes are called *basis* vectors.  What can we say about them?  First, they are a linear combination of the original coordinate system's basis vectors.  Second, they are orthogonal to one another (as is the original coordinate system): $x'1^T x'2 = 0$.  Finally, they are normalized: $x'1^T x'1 = 1$.  \n",
    "\n",
    "In fact, these properties would have held regardless of our choice of $\\theta$.  We just happend to select $\\theta$ such that it got rid of all the information on one of the dimensions.  Why is this useful?  Imagine we wanted to perform classification on this dataset.  A good classifier in the original coordinates would be:\n",
    "$$\n",
    "y_{pred} = \\begin{cases} 1\\;\\mathrm{if}\\;x_1<1,x_2<1 \\\\\n",
    "                  0\\;\\mathrm{else} \\end{cases}\n",
    "$$\n",
    "In the new coordinate system, a good classifier is \n",
    "$$\n",
    "y_{pred} = \\begin{cases} 1\\;\\mathrm{if}\\;x'_1<\\frac{1}{\\sqrt{2}} \\\\\n",
    "                  0\\;\\mathrm{else} \\end{cases}\n",
    "$$\n",
    "Note that the classifier doesn't depend on $x'2$, and we have to do have as many operations: with some preprocessing, we have made our model more efficient by a factor of 2, and reduced the number of parameters by the same.  \n",
    "\n",
    "Can we do this for an arbitrary dataset though?  Consider a similar dataset but with noise added in both $x_1$ and $x_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's generate a synthetic dataset, with two features, and one label\n",
    "x_1 = np.random.randn(21) + 1.0\n",
    "x_2 = x_1.copy()\n",
    "X = np.column_stack((x_1,x_2))\n",
    "X += np.random.randn(*X.shape)*0.1\n",
    "y = (x_1>1) + (x_2>1)\n",
    "\n",
    "plt.scatter(X[:,0],X[:,1],c=y,cmap=plt.cm.prism)\n",
    "plt.xlabel('x_1')\n",
    "plt.ylabel('x_2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what should we do?  There's no longer a clear rotation that reduces the dimensionality of the problem.  However, we can see that not all dimensions have the same amount of information.  Can we still come up with a rotation that gives a coordinate system that is better than this one?  First we need a definition for *better*.  \n",
    "\n",
    "We can take some inspiration from the previous (much easier) dataset.  In that case, we came up with the rotation that allowed us to explain *all* the variability in the data with a single axis.  Clearly we won't be able to do that in this case, but perhaps we can come up with the rotation that allows us to explain as much variability in the data *as possible* with a single axis.  But how can we quantify variability?  Here, we can retreat to a familiar assumption, namely that the data is distributed according to a multivariate normal (MVN).  If we do that, then it's quite simple to quantify the variability in the data through the sample covariance $\\Sigma$.  Let's compute that for this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xbar = X.mean(axis=0)   # Compute the sample mean\n",
    "Sigma = 1./(len(X)-1)*np.dot((X-Xbar).T,(X-Xbar))\n",
    "print (Sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this distribution look like?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0,'./libraries/')\n",
    "from plot_ellipse import plot_ellipse\n",
    "\n",
    "plt.scatter(X[:,0],X[:,1],c=y,cmap=plt.cm.prism)\n",
    "plot_ellipse(Xbar,Sigma,alpha=0.3)\n",
    "plt.xlabel('x_1')\n",
    "plt.ylabel('x_2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this information, it's relatively obvious what the rotation that we come up with should be: the major axis of this ellipse.  Stated alternatively, what's the rotation that makes the covariance matrix diagonal?  To rotate a matrix, we have to rotate both the rows and the columns, and so we need a matrix $R$, that does this:\n",
    "$$\n",
    "R^T \\Sigma R = \\Lambda,\n",
    "$$\n",
    "where $\\Lambda$ is a diagonal matrix with the variance along the ellipse axes on the diagonals.  As it turns out, this is called the eigen-decomposition of a matrix (remember eigenvalues and vectors from linear algebra?).  Thus we can find the rotation matrix $R$ by finding the eigenvectors of the covariance matrix.  Even better, we can compute the variance associated with each one of the resulting axes by finding the eigenvalues.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamda,R = np.linalg.eig(Sigma)  # Find the eigenvalues and vectors of the covariance matrix\n",
    "idx = lamda.argsort()[::-1]     # Sort them, largest to smallest\n",
    "lamda = lamda[idx]\n",
    "R = R[:,idx]\n",
    "\n",
    "X_new = np.dot(X-Xbar,R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_new[:,0],X_new[:,1],c=y,cmap=plt.cm.prism)\n",
    "plot_ellipse(np.zeros((2)),np.diag(lamda),alpha=0.3)\n",
    "plt.xlabel('x\\'_1')\n",
    "plt.ylabel('x\\'_2')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the covariance matrix is diagonal, and most of the variability is on the $x'_1$ axis.  In fact, so much of the variability is explained by $x'_1$, that maybe we don't even need to keep $x'_2$.  We can determine precisely how much of the variability in the data is explained by each axis by normalizing the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variance_fraction = lamda/lamda.sum()\n",
    "print(variance_fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More than 99\\% of the variance is explained by the first principal component!  Thus we could ignore that second dimension and still make a good classifier.  What are we doing by ignoring that second component?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new[:,1] = 0  # Set the information associated with the second PCA to zero\n",
    "X_simplified = np.dot(X_new,R.T) + Xbar\n",
    "plt.scatter(X_simplified[:,0],X_simplified[:,1],c=y,cmap=plt.cm.prism)\n",
    "plt.scatter(X[:,0],X[:,1],marker='^')\n",
    "plt.xlabel('x_1')\n",
    "plt.ylabel('x_2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the points now fall on a line in the original coordinate system.  This is sort of like linear regression, but with an important distinction: in linear regression the thing that was being minimized was the vertical difference between the line and the dependent variable.  In this case, we're drawing the line that minimizes the distance *perpendicular* (or orthogonal) to the line.  This is the result of the fact that this algorithm is unsupervised: no variable is any more special than any other in this case.  However, this does help us to come up with a second interpretation of the first principal component: it is the line that passes closest to the data points.\n",
    "\n",
    "To review, there are two completely equivalent ways to interpret the first principal component:  1, it is the axis along which the maximum amount of variability in the data falls, and 2) it is the line which minimizes the orthogonal distance between itself and the data points.  \n",
    "\n",
    "Note that in this case, because we're considering rotations, we didn't have a choice regarding what the second principal component should be, because it was completely specified by the requirement that it be orthogonal to the first.  This can be seen by the fact that there is only one parameter ($\\theta$) in a 2D rotation matrix.  This is not the case for higher dimensions.  However, the procedure can be applied recursively: the second principal component is the direction that maximizes the explained variance, excluding the variance explained by the first principal component.  Fortunately, the math doesn't change; we still find the principal components by computing the sample covariance, finding its eigenvectors and eigenvalues, and sorting these from largest to smallest.  The eigenvalues from largest to smallest still represent the amount of variance accounted for in the direction of their associated eigenvector/principal component.  \n",
    "\n",
    "Let's see how this applies to a familiar dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "X,X_test,y,y_test = train_test_split(X,y,test_size=0.33)\n",
    "\n",
    "classes = [0,1,2]\n",
    "\n",
    "m = X.shape[0]\n",
    "n = X.shape[1]\n",
    "N = len(classes)\n",
    "\n",
    "fig,axs = plt.subplots(nrows=4,ncols=4)\n",
    "for i in range(4):\n",
    "    for j in range(4):\n",
    "        if i>j:\n",
    "            axs[i,j].scatter(X[:,i],X[:,j],c=y,cmap=plt.cm.prism)\n",
    "            axs[i,j].set_xlabel(iris['feature_names'][i])\n",
    "            axs[i,j].set_ylabel(iris['feature_names'][j])\n",
    "        else: # delete redundant plots\n",
    "            fig.delaxes(axs[i,j])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The iris dataset is 4D.  Let's use PCA to see how many dimensions it *really* has.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xbar = X.mean(axis=0)   # Compute the sample mean\n",
    "Sigma = 1./(len(X)-1)*np.dot((X-Xbar).T,(X-Xbar))\n",
    "print (Sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamda,R = np.linalg.eig(Sigma)  # Find the eigenvalues and vectors of the covariance matrix\n",
    "idx = lamda.argsort()[::-1]     # Sort them, largest to smallest\n",
    "lamda = lamda[idx]\n",
    "R = R[:,idx]\n",
    "\n",
    "X_new = np.dot(X-Xbar,R)\n",
    "X_new_test = np.dot(X_test-Xbar,R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we can look at the values of $\\lambda$ to determine how much of the information content lies in each PC by normalizing it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lamda/lamda.sum()) #Variance fraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more useful metric when deciding what dimensions to keep might be the *cumulative* variance fraction.  As this number approaches one, all of the variability of the data is explained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.cumsum(lamda/lamda.sum())) #Cumulative variance fraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is pretty interesting.  Less than half a percent of the information content of the dataset is explained by the last component.  This is strong evidence that this dataset is actually close to 3D, rather than 4D.  Realistically, the third PC isn't all that interesting either, explaining less than 2\\% of the total data variance.  Let's plot the transformed dataset in the first two dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_new[:,0],X_new[:,1],c=y,cmap=plt.cm.prism)\n",
    "plt.xlabel('PCA 1')\n",
    "plt.ylabel('PCA 2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you stuck a piece of paper through 4D space in a way that it was closest to all the points, then drew a dot on the paper where each point was closest, this is what you would get.  Remember that this has more than 97% of the information content of the original dataset, but with only half the dimensions.  You can easily run a classifier on these transformed observations.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.model_selection import train_test_split\n",
    "nb = GaussianNB()\n",
    "\n",
    "\n",
    "### Confusion matrix for PCA'd iris dataset; two components ###\n",
    "X_new_2 = X_new[:,:2]\n",
    "X_new_2_test = X_new_test[:,:2]\n",
    "nb.fit(X_new_2,y)\n",
    "y_pred = nb.predict(X_new_2_test)\n",
    "print(\"confusion matrix with 2 component PCA\")\n",
    "print(confusion_matrix(y_test,y_pred))\n",
    "\n",
    "### Confusion marix for PCA'd iris dataset; three components ###\n",
    "X_new_3 = X_new[:,:3]\n",
    "X_new_3_test = X_new_test[:,:3]\n",
    "nb.fit(X_new_3,y)\n",
    "y_pred = nb.predict(X_new_3_test)\n",
    "print(\"confusion matrix with 3 component PCA\")\n",
    "print(confusion_matrix(y_test,y_pred))\n",
    "\n",
    "### Confusion marix for PCA'd iris dataset; four components ###\n",
    "nb.fit(X_new,y)\n",
    "y_pred = nb.predict(X_new_test)\n",
    "print(\"confusion matrix with 4 component PCA\")\n",
    "print(confusion_matrix(y_test,y_pred))\n",
    "\n",
    "### Confusion marix for iris dataset; four components ###\n",
    "nb.fit(X,y)\n",
    "y_pred = nb.predict(X_test)\n",
    "print(\"confusion matrix with non-transformed data\")\n",
    "print(confusion_matrix(y_test,y_pred))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the somewhat interesting behavior here with Naive Bayes, in that we get different classifier accuracies between the four component PCA (just different axes, no information lost) and untransformed naive Bayes: this is a result of the naive assumption: since we're not modelling correlations between the features, a rotated dataset produces a different result.  While it didn't work quite as well in this case, it could have just as easily gone the other way.\n",
    "\n",
    "Finally, it can be quite interesting to look at the PCs themselves, particularly the important ones, to see which combinations of features are explaining most of the variability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (R[:,0])  # Print the first PC\n",
    "print (iris.feature_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is like a recipe for explaining variability in the data.  It's saying \"if you want to find the direction that spreads the data out the most, add 35% of sepal length, subtract 7% sepal width, add 85% petal length and 37% petal width\".  \n",
    "\n",
    "While these examples are informative, it's difficult to see the power of PCA on low-dimensional datasets like this.  At the same time, it's difficult to understand what the PCs represent in high dimensional datasets (even the iris example).  An exception to this is *image* data, where the spatial correlation of pixels lets us visualize principal components, even in high dimensions.  To see this, let's look at the MNIST dataset, which is a set of handwritten digits that are 28x28 pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "\n",
    "mnist = datasets.fetch_mldata('MNIST original')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot some random examples of MNIST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(nrows=5,ncols=5)\n",
    "for r in axs:\n",
    "    for ax in r:\n",
    "        ax.imshow(mnist.data[np.random.randint(len(mnist.data)),:].reshape((28,28)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is the first time we've really messed with *image* data, it's worth thinking for just a moment about what the features are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we'll use sklearn's fast implementation of the PCA.  To begin, we'll compute the first 50 components of the PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition\n",
    "\n",
    "pca = decomposition.PCA(n_components=50,copy=True)  \n",
    "pca.fit(mnist.data)\n",
    "X = pca.transform(mnist.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that each PC is just a linear combination of the input features.  Since the input features in this case are pixel values, we can transform the PCs back into 28x28 images, and plot the intensity of each.  This produces something called an eigendigit.  Let's plot the first 9:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(nrows=3,ncols=3)\n",
    "counter = 0\n",
    "for r in axs:\n",
    "    for ax in r:\n",
    "        ax.imshow(pca.components_[counter,:].reshape((28,28)))\n",
    "        counter+=1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These represent the best possible basis for describing digits, which is to say that if you only had image that you had to scale to describe all of the digits in MNIST, that O thing would be the one.  If you could only use linear combinations of two images to describe all of MNIST, it would be the first two.  This seems absurd, but think about what the original basis set is for images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(nrows=3,ncols=3)\n",
    "counter = 0\n",
    "for r in axs:\n",
    "    for ax in r:\n",
    "        basis = np.zeros((28,28))\n",
    "        basis[counter,0] = 1\n",
    "        ax.imshow(basis)\n",
    "        counter+=1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the PCA basis is more informative!  Just how informative? Let's look at the cumulative variance accounted for by the PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (np.cumsum(pca.explained_variance_ratio_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That first component explains a full 10% of data variability.  The first 50 components explain 82% of the variability.  Note that this data is obviously far more multi-faceted and complex than the iris dataset.  Not only are there more classes (10 digits), but also there is more intra-class variability (the variability between how different people write sixes, for example).  How many PCs do we need in order to get up to 95% variance?  sklearn lets us specify a variance ratio, and it will keep computing PCs until it gets there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = decomposition.PCA(0.95,copy=True,whiten=False)\n",
    "pca.fit(mnist.data)\n",
    "X = pca.transform(mnist.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cumulative_variance_ratio = np.cumsum(pca.explained_variance_ratio_)\n",
    "plt.plot(1-cumulative_variance_ratio)\n",
    "plt.ylabel('Unexplained variance remaining')\n",
    "plt.show()\n",
    "154/784"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is known as a scree plot (for the shape of a mountainside).  It represents how much information you're throwing away for a given model simplification.  At L=154 PCs, we're throwing away 5% of the data.  But that's a huge compression! We still have 95% of the information at $154/784 \\approx 20\\%$ of the number of features.\n",
    "\n",
    "When we transform our data into representation by the PCs, our new data representation is a vector of length 154, which isn't conducive to visualization.  However, we can reconstruct what our data looks like under the basis simplification, just like we did when we zeroed out a column in the synthetic example, then transformed back to the original coordinates.  We do this by taking a linear combination of the PCs, with the coefficients given by the entries in the data matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(nrows=5,ncols=2)\n",
    "fig.set_size_inches(12,30)\n",
    "\n",
    "data_index = 15000\n",
    "for i,n_components in enumerate([1,5,10,50,150]):\n",
    "    X_reconstructed = 0\n",
    "    counter = 0 \n",
    "    for c,l in zip(pca.components_,X[data_index]):\n",
    "        if counter<n_components:\n",
    "            X_reconstructed += c*l\n",
    "        counter += 1\n",
    "    X_reconstructed += pca.mean_\n",
    "    X_reconstructed = X_reconstructed.reshape((28,28))\n",
    "\n",
    "    axs[i,0].imshow(mnist.data[data_index,:].reshape((28,28)),vmin=0,vmax=255)\n",
    "    axs[i,1].imshow(X_reconstructed,vmin=0,vmax=255)\n",
    "    axs[i,1].set_title(\"L=\"+str(n_components))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use this method to come up with a fairly efficient classifier for these digits in a very similar way to what we did before.  We can even use Gaussian Naive Bayes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "\n",
    "mnist = datasets.fetch_mldata('MNIST original')\n",
    "X,X_test,y,y_test = train_test_split(mnist.data,mnist.target,test_size=0.33)\n",
    "\n",
    "pca = decomposition.PCA(35,copy=True,whiten=True)\n",
    "pca.fit(X)\n",
    "X_pca = pca.transform(X)\n",
    "X_pca_test = pca.transform(X_test)\n",
    "\n",
    "nb = GaussianNB()\n",
    "### Confusion marix for mnist dataset###\n",
    "nb.fit(X_pca,y)\n",
    "y_pred = nb.predict(X_pca_test)\n",
    "c = confusion_matrix(y_test,y_pred)\n",
    "accuracy = np.diagonal(c).sum()/c.sum()\n",
    "print(\"Confusion matrix: \\n\",c)\n",
    "print(\"Classification accuracy: \",accuracy)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just out of curiosity, what did this classifier get wrong?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_wrong_value = np.random.choice(np.where(y_pred!=y_test)[0])\n",
    "fig,axs = plt.subplots(nrows=1,ncols=1)\n",
    "axs.imshow(X_test[random_wrong_value].reshape((28,28)))\n",
    "print(\"Classified as :\",y_pred[random_wrong_value])\n",
    "print(\"Actual class :\",y_test[random_wrong_value])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.choice?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}