Chapter_12_sec_12.5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.5 Lab: Unsupervised Learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from numpy.linalg import svd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.cluster import KMeans\n",
    "from scipy.cluster import hierarchy\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.5.1 Principal Components Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "USArrests = pd.read_csv('./data/USArrests.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "USArrests.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pandas has a built-in function to get the mean and variance of each column\n",
    "print(USArrests.mean())\n",
    "print(USArrests.var())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc = StandardScaler()\n",
    "X = pd.DataFrame(sc.fit_transform(USArrests), index=USArrests.index, columns=USArrests.columns)\n",
    "# The loading vectors (i.e. these are the projection of the data onto the principal components)\n",
    "pca_loadings = pd.DataFrame(PCA().fit(X).components_.T, index=USArrests.columns, columns=['V1', 'V2', 'V3', 'V4'])\n",
    "pca_loadings\n",
    "\n",
    "\"\"\" \n",
    "Depends on the version of python/module, you may see a flipped loading vector in signs. \n",
    "This is normal because the orientation of the principal components is not deterministic. \n",
    "\"\"\"\n",
    "# X1=pd.DataFrame(sc.inverse_transform(X), index=USArrests.index, columns=USArrests.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the PCA model and transform X to get the principal components\n",
    "pca = PCA()\n",
    "df_plot = pd.DataFrame(pca.fit_transform(X), columns=['PC1', 'PC2', 'PC3', 'PC4'], index=X.index)\n",
    "df_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig , ax1 = plt.subplots(figsize=(9,7))\n",
    "\n",
    "ax1.set_xlim(-3.5,3.5)\n",
    "ax1.set_ylim(-3.5,3.5)\n",
    "\n",
    "# plot Principal Components 1 and 2\n",
    "for i in df_plot.index:\n",
    "    ax1.annotate(i, (df_plot.PC1.loc[i], -df_plot.PC2.loc[i]), ha='center')\n",
    "\n",
    "# plot reference lines\n",
    "ax1.hlines(0,-3.5,3.5, linestyles='dotted', colors='grey')\n",
    "ax1.vlines(0,-3.5,3.5, linestyles='dotted', colors='grey')\n",
    "\n",
    "ax1.set_xlabel('First Principal Component')\n",
    "ax1.set_ylabel('Second Principal Component')\n",
    "    \n",
    "# plot Principal Component loading vectors, using a second y-axis.\n",
    "ax2 = ax1.twinx().twiny() \n",
    "\n",
    "ax2.set_ylim(-1,1)\n",
    "ax2.set_xlim(-1,1)\n",
    "ax2.tick_params(axis='y', colors='orange')\n",
    "ax2.set_xlabel('Principal Component loading vectors', color='orange')\n",
    "\n",
    "# plot labels for vectors. Variable 'a' is a small offset parameter to separate arrow tip and text.\n",
    "a = 1.07  \n",
    "for i in pca_loadings[['V1', 'V2']].index:\n",
    "    ax2.annotate(i, (pca_loadings.V1.loc[i]*a, -pca_loadings.V2.loc[i]*a), color='orange')\n",
    "\n",
    "# plot vectors\n",
    "ax2.arrow(0,0,pca_loadings.V1[0], -pca_loadings.V2[0])\n",
    "ax2.arrow(0,0,pca_loadings.V1[1], -pca_loadings.V2[1])\n",
    "ax2.arrow(0,0,pca_loadings.V1[2], -pca_loadings.V2[2])\n",
    "ax2.arrow(0,0,pca_loadings.V1[3], -pca_loadings.V2[3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in previous chapter, we talked about PCR. In that case, we could use the downstream task's (i.e. regression RMSE) \n",
    "# performance to select the hyperparameters (i.e. # number of PCs).\n",
    "# here let us use the portion of explained variance to select the number of PCs. Those info is available in the pca object.\n",
    "print(pca.explained_variance_)\n",
    "print(pca.explained_variance_ratio_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,5))\n",
    "\n",
    "plt.plot([1,2,3,4], pca.explained_variance_ratio_, '-o', label='Individual component')\n",
    "plt.plot([1,2,3,4], np.cumsum(pca.explained_variance_ratio_), '-s', label='Cumulative')\n",
    "\n",
    "plt.ylabel('Proportion of Variance Explained')\n",
    "plt.xlabel('Principal Component')\n",
    "plt.xlim(0.75,4.25)\n",
    "plt.ylim(0,1.05)\n",
    "plt.xticks([1,2,3,4])\n",
    "plt.legend(loc=2)\n",
    "plt.show()\n",
    "\n",
    "\"\"\"\n",
    "In this case, if we want to preserve 80% of variance of the data, we need to select 2 PCs.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.5.2 Matrix Completion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "I am happy seeing this lab added. SVD seems pretty heavy in theory/math, but this has lots of application in real problems, \n",
    "such as recommendation systems, clustering, outlier smoothing, and so on.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run sigular value decomposition on the data (SVD)\n",
    "u, s, vh = svd(X, full_matrices=False)\n",
    "u.shape, s.shape, vh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this vh will be the principal components similar to pca.components_ (up to an unimportant sign flip)\n",
    "# The matrix u is equivalent to the matrix of standardized scores, and the standard deviations are in the vector s.\n",
    "print(vh)\n",
    "print ('-------')\n",
    "print(pca.components_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruction based on full SVD\n",
    "np.allclose(X, np.dot(u * s, vh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruction based on reduced SVD\n",
    "num_components = 3\n",
    "recovered = pd.DataFrame(np.dot(u[:, :num_components] * s[:num_components,], vh[:num_components,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(recovered.head(n=2))\n",
    "print(X.head(n=2))\n",
    "\n",
    "\"\"\"\n",
    "Change the num_components from 1 to 4 and see how the reconstruction error changes.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.5.3 Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "K-means clustering "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate data\n",
    "np.random.seed(21)\n",
    "X = np.random.standard_normal((50,2))\n",
    "X[:25,0] = X[:25,0]+3\n",
    "X[:25,1] = X[:25,1]-4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_clusters = 2\n",
    "km1 = KMeans(n_clusters=n_clusters, n_init=20)\n",
    "km1.fit(X)\n",
    "\n",
    "n_clusters = 3\n",
    "km2 = KMeans(n_clusters=n_clusters, n_init=20)\n",
    "km2.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(km1.labels_)\n",
    "print(dir(km1)) # we can use dir to see other saved attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(14,5))\n",
    "\n",
    "ax1.scatter(X[:,0], X[:,1], s=40, c=km1.labels_, cmap=plt.cm.prism) \n",
    "ax1.set_title('K-Means Clustering Results with K=2')\n",
    "ax1.scatter(km1.cluster_centers_[:,0], km1.cluster_centers_[:,1], marker='+', s=100, c='k', linewidth=2)\n",
    "\n",
    "ax2.scatter(X[:,0], X[:,1], s=40, c=km2.labels_, cmap=plt.cm.prism) \n",
    "ax2.set_title('K-Means Clustering Results with K=3')\n",
    "ax2.scatter(km2.cluster_centers_[:,0], km2.cluster_centers_[:,1], marker='+', s=100, c='k', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hierarchical Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1,ax2,ax3) = plt.subplots(3,1, figsize=(15,18))\n",
    "\n",
    "for linkage, cluster, ax in zip([hierarchy.complete(X), hierarchy.average(X), hierarchy.single(X)], ['c1','c2','c3'],\n",
    "                                [ax1,ax2,ax3]):\n",
    "    cluster = hierarchy.dendrogram(linkage, ax=ax, color_threshold=0)\n",
    "\n",
    "ax1.set_title('Complete Linkage')\n",
    "ax2.set_title('Average Linkage')\n",
    "ax3.set_title('Single Linkage')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.5.4 NCI60 Data Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA on the NCI60 Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I was not able to make the json work, so I went back to R and saved the data and label separately.\n",
    "X = pd.read_csv('./data/NCI60_data.csv')\n",
    "y = pd.read_csv('./data/NCI60_labs.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca2 = PCA()\n",
    "X_standardized = StandardScaler().fit_transform(X)\n",
    "df2_plot = pd.DataFrame(pca2.fit_transform(X_standardized))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(15,6))\n",
    "\n",
    "color_idx = pd.factorize(y.iloc[:, 0])[0]\n",
    "cmap = plt.cm.hsv\n",
    "\n",
    "# left plot\n",
    "ax1.scatter(df2_plot.iloc[:,0], -df2_plot.iloc[:,1], c=color_idx, cmap=cmap, alpha=0.5, s=50)\n",
    "ax1.set_ylabel('Principal Component 2')\n",
    "\n",
    "# right plot\n",
    "ax2.scatter(df2_plot.iloc[:,0], df2_plot.iloc[:,2], c=color_idx, cmap=cmap, alpha=0.5, s=50)\n",
    "ax2.set_ylabel('Principal Component 3')\n",
    "\n",
    "# custom legend for the classes (y) since we do not create scatter plots per class (which could have their own labels).\n",
    "handles = []\n",
    "labels = pd.factorize(y.iloc[:, 0].unique())\n",
    "norm = mpl.colors.Normalize(vmin=0.0, vmax=14.0)\n",
    "\n",
    "for i, v in zip(labels[0], labels[1]):\n",
    "    handles.append(mpl.patches.Patch(color=cmap(norm(i)), label=v, alpha=0.5))\n",
    "\n",
    "ax2.legend(handles=handles, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# xlabel for both plots\n",
    "for ax in fig.axes:\n",
    "    ax.set_xlabel('Principal Component 1')   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame([df2_plot.iloc[:,:5].std(axis=0, ddof=0).array,\n",
    "              pca2.explained_variance_ratio_[:5],\n",
    "              np.cumsum(pca2.explained_variance_ratio_[:5])],\n",
    "             index=['Standard Deviation', 'Proportion of Variance', 'Cumulative Proportion'],\n",
    "             columns=['PC1', 'PC2', 'PC3', 'PC4', 'PC5'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2_plot.iloc[:,:10].var(axis=0, ddof=0).plot(kind='bar', rot=0)\n",
    "plt.ylabel('Variances')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig , (ax1,ax2) = plt.subplots(1,2, figsize=(15,5))\n",
    "\n",
    "# left plot\n",
    "ax1.plot(pca2.explained_variance_ratio_, '-o')\n",
    "ax1.set_ylabel('Proportion of Variance Explained')\n",
    "ax1.set_ylim(ymin=-0.01)\n",
    "\n",
    "# right plot\n",
    "ax2.plot(np.cumsum(pca2.explained_variance_ratio_), '-ro')\n",
    "ax2.set_ylabel('Cumulative Proportion of Variance Explained')\n",
    "ax2.set_ylim(ymax=1.05)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    ax.set_xlabel('Principal Component')\n",
    "    ax.set_xlim(-1,65) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clustering the Observations of the NCI60 Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc = StandardScaler()\n",
    "X_standardized = pd.DataFrame(sc.fit_transform(X), index=y.iloc[:, 0], columns=X.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize=(20,20))\n",
    "\n",
    "for linkage, cluster, ax in zip([hierarchy.complete(X_standardized), hierarchy.average(X), hierarchy.single(X_standardized)],\n",
    "                                ['c1','c2','c3'],\n",
    "                                [ax1,ax2,ax3]):\n",
    "    cluster = hierarchy.dendrogram(linkage, labels=X_standardized.index, orientation='right', color_threshold=0, leaf_font_size=10, ax=ax)\n",
    "\n",
    "ax1.set_title('Complete Linkage')\n",
    "ax2.set_title('Average Linkage')\n",
    "ax3.set_title('Single Linkage')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,20))\n",
    "cut4 = hierarchy.dendrogram(hierarchy.complete(X_standardized),\n",
    "                            labels=X_standardized.index, orientation='right', color_threshold=140, leaf_font_size=10)\n",
    "plt.vlines(140,0,plt.gca().yaxis.get_data_interval()[1], colors='r', linestyles='dashed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kmeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(21)\n",
    "km3 = KMeans(n_clusters=4, n_init=50)\n",
    "km3.fit(X_standardized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "km3.labels_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine with PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,20))\n",
    "pca_cluster = hierarchy.dendrogram(hierarchy.complete(X_standardized), labels=X_standardized.index,\n",
    "orientation='right', color_threshold=100, leaf_font_size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hierarchy based on Principal Components 1 to 5\n",
    "plt.figure(figsize=(10,20))\n",
    "pca_cluster = hierarchy.dendrogram(hierarchy.complete(df2_plot.iloc[:,:5]), labels=X_standardized.index,\n",
    "orientation='right', color_threshold=100, leaf_font_size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 12"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}