{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "80f16ff6",
   "metadata": {},
   "source": [
    "\n",
    "# Chapter 12\n",
    "\n",
    " # Lab: Unsupervised Learning\n",
    "In this lab we demonstrate PCA and clustering on several datasets.\n",
    "As in other labs, we import some of our libraries at this top\n",
    "level. This makes the code more readable, as scanning the first few\n",
    "lines of the notebook tell us what libraries are used in this\n",
    "notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "24559be0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:08.163850Z",
     "iopub.status.busy": "2023-07-26T05:18:08.163726Z",
     "iopub.status.idle": "2023-07-26T05:18:09.209261Z",
     "shell.execute_reply": "2023-07-26T05:18:09.208715Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.datasets import get_rdataset\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from ISLP import load_data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59b24a4b",
   "metadata": {},
   "source": [
    "We also collect the new imports\n",
    "needed for this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "06fff57d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:09.212356Z",
     "iopub.status.busy": "2023-07-26T05:18:09.212070Z",
     "iopub.status.idle": "2023-07-26T05:18:09.327611Z",
     "shell.execute_reply": "2023-07-26T05:18:09.327246Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import \\\n",
    "     (KMeans,\n",
    "      AgglomerativeClustering)\n",
    "from scipy.cluster.hierarchy import \\\n",
    "     (dendrogram,\n",
    "      cut_tree)\n",
    "from ISLP.cluster import compute_linkage\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f091de06",
   "metadata": {},
   "source": [
    "## Principal Components Analysis\n",
    "In this lab, we perform PCA on  `USArrests`, a data set in the\n",
    "`R` computing environment.\n",
    "We retrieve the data using `get_rdataset()`, which can fetch data from\n",
    "many standard `R` packages.\n",
    "\n",
    "The rows of the data set contain the 50 states, in alphabetical order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f425e07e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:09.329909Z",
     "iopub.status.busy": "2023-07-26T05:18:09.329763Z",
     "iopub.status.idle": "2023-07-26T05:18:10.559326Z",
     "shell.execute_reply": "2023-07-26T05:18:10.558692Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Murder</th>\n",
       "      <th>Assault</th>\n",
       "      <th>UrbanPop</th>\n",
       "      <th>Rape</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alabama</th>\n",
       "      <td>13.2</td>\n",
       "      <td>236</td>\n",
       "      <td>58</td>\n",
       "      <td>21.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Alaska</th>\n",
       "      <td>10.0</td>\n",
       "      <td>263</td>\n",
       "      <td>48</td>\n",
       "      <td>44.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arizona</th>\n",
       "      <td>8.1</td>\n",
       "      <td>294</td>\n",
       "      <td>80</td>\n",
       "      <td>31.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arkansas</th>\n",
       "      <td>8.8</td>\n",
       "      <td>190</td>\n",
       "      <td>50</td>\n",
       "      <td>19.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>California</th>\n",
       "      <td>9.0</td>\n",
       "      <td>276</td>\n",
       "      <td>91</td>\n",
       "      <td>40.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Colorado</th>\n",
       "      <td>7.9</td>\n",
       "      <td>204</td>\n",
       "      <td>78</td>\n",
       "      <td>38.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Connecticut</th>\n",
       "      <td>3.3</td>\n",
       "      <td>110</td>\n",
       "      <td>77</td>\n",
       "      <td>11.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Delaware</th>\n",
       "      <td>5.9</td>\n",
       "      <td>238</td>\n",
       "      <td>72</td>\n",
       "      <td>15.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Florida</th>\n",
       "      <td>15.4</td>\n",
       "      <td>335</td>\n",
       "      <td>80</td>\n",
       "      <td>31.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Georgia</th>\n",
       "      <td>17.4</td>\n",
       "      <td>211</td>\n",
       "      <td>60</td>\n",
       "      <td>25.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Hawaii</th>\n",
       "      <td>5.3</td>\n",
       "      <td>46</td>\n",
       "      <td>83</td>\n",
       "      <td>20.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Idaho</th>\n",
       "      <td>2.6</td>\n",
       "      <td>120</td>\n",
       "      <td>54</td>\n",
       "      <td>14.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Illinois</th>\n",
       "      <td>10.4</td>\n",
       "      <td>249</td>\n",
       "      <td>83</td>\n",
       "      <td>24.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indiana</th>\n",
       "      <td>7.2</td>\n",
       "      <td>113</td>\n",
       "      <td>65</td>\n",
       "      <td>21.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Iowa</th>\n",
       "      <td>2.2</td>\n",
       "      <td>56</td>\n",
       "      <td>57</td>\n",
       "      <td>11.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kansas</th>\n",
       "      <td>6.0</td>\n",
       "      <td>115</td>\n",
       "      <td>66</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kentucky</th>\n",
       "      <td>9.7</td>\n",
       "      <td>109</td>\n",
       "      <td>52</td>\n",
       "      <td>16.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Louisiana</th>\n",
       "      <td>15.4</td>\n",
       "      <td>249</td>\n",
       "      <td>66</td>\n",
       "      <td>22.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maine</th>\n",
       "      <td>2.1</td>\n",
       "      <td>83</td>\n",
       "      <td>51</td>\n",
       "      <td>7.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maryland</th>\n",
       "      <td>11.3</td>\n",
       "      <td>300</td>\n",
       "      <td>67</td>\n",
       "      <td>27.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Massachusetts</th>\n",
       "      <td>4.4</td>\n",
       "      <td>149</td>\n",
       "      <td>85</td>\n",
       "      <td>16.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Michigan</th>\n",
       "      <td>12.1</td>\n",
       "      <td>255</td>\n",
       "      <td>74</td>\n",
       "      <td>35.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Minnesota</th>\n",
       "      <td>2.7</td>\n",
       "      <td>72</td>\n",
       "      <td>66</td>\n",
       "      <td>14.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mississippi</th>\n",
       "      <td>16.1</td>\n",
       "      <td>259</td>\n",
       "      <td>44</td>\n",
       "      <td>17.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Missouri</th>\n",
       "      <td>9.0</td>\n",
       "      <td>178</td>\n",
       "      <td>70</td>\n",
       "      <td>28.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Montana</th>\n",
       "      <td>6.0</td>\n",
       "      <td>109</td>\n",
       "      <td>53</td>\n",
       "      <td>16.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nebraska</th>\n",
       "      <td>4.3</td>\n",
       "      <td>102</td>\n",
       "      <td>62</td>\n",
       "      <td>16.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nevada</th>\n",
       "      <td>12.2</td>\n",
       "      <td>252</td>\n",
       "      <td>81</td>\n",
       "      <td>46.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>New Hampshire</th>\n",
       "      <td>2.1</td>\n",
       "      <td>57</td>\n",
       "      <td>56</td>\n",
       "      <td>9.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>New Jersey</th>\n",
       "      <td>7.4</td>\n",
       "      <td>159</td>\n",
       "      <td>89</td>\n",
       "      <td>18.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>New Mexico</th>\n",
       "      <td>11.4</td>\n",
       "      <td>285</td>\n",
       "      <td>70</td>\n",
       "      <td>32.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>New York</th>\n",
       "      <td>11.1</td>\n",
       "      <td>254</td>\n",
       "      <td>86</td>\n",
       "      <td>26.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>North Carolina</th>\n",
       "      <td>13.0</td>\n",
       "      <td>337</td>\n",
       "      <td>45</td>\n",
       "      <td>16.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>North Dakota</th>\n",
       "      <td>0.8</td>\n",
       "      <td>45</td>\n",
       "      <td>44</td>\n",
       "      <td>7.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ohio</th>\n",
       "      <td>7.3</td>\n",
       "      <td>120</td>\n",
       "      <td>75</td>\n",
       "      <td>21.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Oklahoma</th>\n",
       "      <td>6.6</td>\n",
       "      <td>151</td>\n",
       "      <td>68</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Oregon</th>\n",
       "      <td>4.9</td>\n",
       "      <td>159</td>\n",
       "      <td>67</td>\n",
       "      <td>29.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Pennsylvania</th>\n",
       "      <td>6.3</td>\n",
       "      <td>106</td>\n",
       "      <td>72</td>\n",
       "      <td>14.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rhode Island</th>\n",
       "      <td>3.4</td>\n",
       "      <td>174</td>\n",
       "      <td>87</td>\n",
       "      <td>8.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>South Carolina</th>\n",
       "      <td>14.4</td>\n",
       "      <td>279</td>\n",
       "      <td>48</td>\n",
       "      <td>22.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>South Dakota</th>\n",
       "      <td>3.8</td>\n",
       "      <td>86</td>\n",
       "      <td>45</td>\n",
       "      <td>12.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tennessee</th>\n",
       "      <td>13.2</td>\n",
       "      <td>188</td>\n",
       "      <td>59</td>\n",
       "      <td>26.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Texas</th>\n",
       "      <td>12.7</td>\n",
       "      <td>201</td>\n",
       "      <td>80</td>\n",
       "      <td>25.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Utah</th>\n",
       "      <td>3.2</td>\n",
       "      <td>120</td>\n",
       "      <td>80</td>\n",
       "      <td>22.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Vermont</th>\n",
       "      <td>2.2</td>\n",
       "      <td>48</td>\n",
       "      <td>32</td>\n",
       "      <td>11.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Virginia</th>\n",
       "      <td>8.5</td>\n",
       "      <td>156</td>\n",
       "      <td>63</td>\n",
       "      <td>20.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Washington</th>\n",
       "      <td>4.0</td>\n",
       "      <td>145</td>\n",
       "      <td>73</td>\n",
       "      <td>26.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>West Virginia</th>\n",
       "      <td>5.7</td>\n",
       "      <td>81</td>\n",
       "      <td>39</td>\n",
       "      <td>9.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Wisconsin</th>\n",
       "      <td>2.6</td>\n",
       "      <td>53</td>\n",
       "      <td>66</td>\n",
       "      <td>10.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Wyoming</th>\n",
       "      <td>6.8</td>\n",
       "      <td>161</td>\n",
       "      <td>60</td>\n",
       "      <td>15.6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Murder  Assault  UrbanPop  Rape\n",
       "Alabama           13.2      236        58  21.2\n",
       "Alaska            10.0      263        48  44.5\n",
       "Arizona            8.1      294        80  31.0\n",
       "Arkansas           8.8      190        50  19.5\n",
       "California         9.0      276        91  40.6\n",
       "Colorado           7.9      204        78  38.7\n",
       "Connecticut        3.3      110        77  11.1\n",
       "Delaware           5.9      238        72  15.8\n",
       "Florida           15.4      335        80  31.9\n",
       "Georgia           17.4      211        60  25.8\n",
       "Hawaii             5.3       46        83  20.2\n",
       "Idaho              2.6      120        54  14.2\n",
       "Illinois          10.4      249        83  24.0\n",
       "Indiana            7.2      113        65  21.0\n",
       "Iowa               2.2       56        57  11.3\n",
       "Kansas             6.0      115        66  18.0\n",
       "Kentucky           9.7      109        52  16.3\n",
       "Louisiana         15.4      249        66  22.2\n",
       "Maine              2.1       83        51   7.8\n",
       "Maryland          11.3      300        67  27.8\n",
       "Massachusetts      4.4      149        85  16.3\n",
       "Michigan          12.1      255        74  35.1\n",
       "Minnesota          2.7       72        66  14.9\n",
       "Mississippi       16.1      259        44  17.1\n",
       "Missouri           9.0      178        70  28.2\n",
       "Montana            6.0      109        53  16.4\n",
       "Nebraska           4.3      102        62  16.5\n",
       "Nevada            12.2      252        81  46.0\n",
       "New Hampshire      2.1       57        56   9.5\n",
       "New Jersey         7.4      159        89  18.8\n",
       "New Mexico        11.4      285        70  32.1\n",
       "New York          11.1      254        86  26.1\n",
       "North Carolina    13.0      337        45  16.1\n",
       "North Dakota       0.8       45        44   7.3\n",
       "Ohio               7.3      120        75  21.4\n",
       "Oklahoma           6.6      151        68  20.0\n",
       "Oregon             4.9      159        67  29.3\n",
       "Pennsylvania       6.3      106        72  14.9\n",
       "Rhode Island       3.4      174        87   8.3\n",
       "South Carolina    14.4      279        48  22.5\n",
       "South Dakota       3.8       86        45  12.8\n",
       "Tennessee         13.2      188        59  26.9\n",
       "Texas             12.7      201        80  25.5\n",
       "Utah               3.2      120        80  22.9\n",
       "Vermont            2.2       48        32  11.2\n",
       "Virginia           8.5      156        63  20.7\n",
       "Washington         4.0      145        73  26.2\n",
       "West Virginia      5.7       81        39   9.3\n",
       "Wisconsin          2.6       53        66  10.8\n",
       "Wyoming            6.8      161        60  15.6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USArrests = get_rdataset('USArrests').data\n",
    "USArrests\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2942890",
   "metadata": {},
   "source": [
    "The columns of the data set contain the four variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b127d014",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.561457Z",
     "iopub.status.busy": "2023-07-26T05:18:10.561326Z",
     "iopub.status.idle": "2023-07-26T05:18:10.564289Z",
     "shell.execute_reply": "2023-07-26T05:18:10.563966Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Murder', 'Assault', 'UrbanPop', 'Rape'], dtype='object')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USArrests.columns\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68718331",
   "metadata": {},
   "source": [
    "We first briefly examine the data. We notice that the variables have vastly different means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c7343f72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.566649Z",
     "iopub.status.busy": "2023-07-26T05:18:10.566503Z",
     "iopub.status.idle": "2023-07-26T05:18:10.570319Z",
     "shell.execute_reply": "2023-07-26T05:18:10.569890Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Murder        7.788\n",
       "Assault     170.760\n",
       "UrbanPop     65.540\n",
       "Rape         21.232\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USArrests.mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "171a1ee0",
   "metadata": {},
   "source": [
    "Dataframes have several useful methods for computing\n",
    "column-wise summaries. We can also examine the\n",
    "variance of the four variables using the `var()`  method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "34501140",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.572610Z",
     "iopub.status.busy": "2023-07-26T05:18:10.572485Z",
     "iopub.status.idle": "2023-07-26T05:18:10.575894Z",
     "shell.execute_reply": "2023-07-26T05:18:10.575543Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Murder        18.970465\n",
       "Assault     6945.165714\n",
       "UrbanPop     209.518776\n",
       "Rape          87.729159\n",
       "dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "USArrests.var()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5634db88",
   "metadata": {},
   "source": [
    "Not surprisingly, the variables also have vastly different variances.\n",
    "The `UrbanPop` variable measures the percentage of the population\n",
    "in each state living in an urban area, which is not a comparable\n",
    "number to the number of rapes in each state per 100,000 individuals.\n",
    "PCA looks for derived variables that account for most of the variance in the data set.\n",
    "If we do not scale the variables before performing PCA, then the principal components\n",
    "would mostly be driven by the\n",
    "`Assault` variable, since it has by far the largest\n",
    "variance.  So if the variables are measured in different units or vary widely in scale, it is recommended to standardize the variables to have standard deviation one before performing PCA.\n",
    "Typically we set the means to zero as well.\n",
    "\n",
    "This scaling can be done via the `StandardScaler()` transform imported above. We first `fit` the\n",
    "scaler, which computes the necessary means and standard\n",
    "deviations and then apply it to our data using the\n",
    "`transform` method. As before, we combine these steps using the `fit_transform()` method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "daf119e8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.578199Z",
     "iopub.status.busy": "2023-07-26T05:18:10.578025Z",
     "iopub.status.idle": "2023-07-26T05:18:10.581928Z",
     "shell.execute_reply": "2023-07-26T05:18:10.581446Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler(with_std=True,\n",
    "                        with_mean=True)\n",
    "USArrests_scaled = scaler.fit_transform(USArrests)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd6e5f5b",
   "metadata": {},
   "source": [
    "Having scaled the data, we can then\n",
    "perform principal components analysis using the `PCA()` transform\n",
    "from the `sklearn.decomposition` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a0eda7c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.584181Z",
     "iopub.status.busy": "2023-07-26T05:18:10.584070Z",
     "iopub.status.idle": "2023-07-26T05:18:10.586440Z",
     "shell.execute_reply": "2023-07-26T05:18:10.586052Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "pcaUS = PCA()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5758ee5",
   "metadata": {},
   "source": [
    "(By default, the `PCA()`  transform centers the variables to have\n",
    "mean zero though it does not scale them.) The transform `pcaUS`\n",
    "can be used to find the PCA\n",
    "`scores` returned by `fit()`. Once the `fit` method has been called, the `pcaUS` object also contains a number of useful quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1430fb3c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.588957Z",
     "iopub.status.busy": "2023-07-26T05:18:10.588809Z",
     "iopub.status.idle": "2023-07-26T05:18:10.593315Z",
     "shell.execute_reply": "2023-07-26T05:18:10.592971Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>PCA()</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">PCA</label><div class=\"sk-toggleable__content\"><pre>PCA()</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "PCA()"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.fit(USArrests_scaled)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56890b11",
   "metadata": {},
   "source": [
    "After fitting, the `mean_` attribute corresponds to the means\n",
    "of the variables. In this case, since we centered and scaled the data with\n",
    "`scaler()` the means will all be 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6131d8d1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.595542Z",
     "iopub.status.busy": "2023-07-26T05:18:10.595388Z",
     "iopub.status.idle": "2023-07-26T05:18:10.598527Z",
     "shell.execute_reply": "2023-07-26T05:18:10.598057Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-7.10542736e-17,  1.38777878e-16, -4.39648318e-16,  8.59312621e-16])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.mean_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1c34e7",
   "metadata": {},
   "source": [
    "The scores can be computed using the `transform()` method\n",
    "of `pcaUS` after it has been fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "08246aad",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.600656Z",
     "iopub.status.busy": "2023-07-26T05:18:10.600536Z",
     "iopub.status.idle": "2023-07-26T05:18:10.603024Z",
     "shell.execute_reply": "2023-07-26T05:18:10.602498Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "scores = pcaUS.transform(USArrests_scaled)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "574ed47f",
   "metadata": {},
   "source": [
    "We will plot these scores a bit further down.\n",
    "The `components_` attribute provides the principal component loadings:\n",
    "each row of `pcaUS.components_` contains the corresponding\n",
    "principal component loading vector.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b682b632",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.605166Z",
     "iopub.status.busy": "2023-07-26T05:18:10.605000Z",
     "iopub.status.idle": "2023-07-26T05:18:10.607923Z",
     "shell.execute_reply": "2023-07-26T05:18:10.607498Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.53589947,  0.58318363,  0.27819087,  0.54343209],\n",
       "       [ 0.41818087,  0.1879856 , -0.87280619, -0.16731864],\n",
       "       [-0.34123273, -0.26814843, -0.37801579,  0.81777791],\n",
       "       [ 0.6492278 , -0.74340748,  0.13387773,  0.08902432]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.components_ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc8bfc99",
   "metadata": {},
   "source": [
    "The `biplot`  is a common visualization method used with\n",
    "PCA. It is not built in as a standard\n",
    "part of `sklearn`, though there are python\n",
    "packages that do produce such plots. Here we\n",
    "make a simple biplot manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c165e990",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.610624Z",
     "iopub.status.busy": "2023-07-26T05:18:10.610442Z",
     "iopub.status.idle": "2023-07-26T05:18:10.736284Z",
     "shell.execute_reply": "2023-07-26T05:18:10.735735Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34480 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "i, j = 0, 1 # which components\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "ax.scatter(scores[:,0], scores[:,1])\n",
    "ax.set_xlabel('PC%d' % (i+1))\n",
    "ax.set_ylabel('PC%d' % (j+1))\n",
    "for k in range(pcaUS.components_.shape[1]):\n",
    "    ax.arrow(0, 0, pcaUS.components_[i,k], pcaUS.components_[j,k])\n",
    "    ax.text(pcaUS.components_[i,k],\n",
    "            pcaUS.components_[j,k],\n",
    "            USArrests.columns[k])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaef3e98",
   "metadata": {},
   "source": [
    "Notice that this figure is a reflection of Figure 12.1 through the $y$-axis. Recall that the\n",
    "principal components are only unique up to a sign change, so we can\n",
    "reproduce that figure by flipping the\n",
    "signs of the second set of scores and loadings.\n",
    "We also increase the length of the arrows to emphasize the loadings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "848c9f35",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.738926Z",
     "iopub.status.busy": "2023-07-26T05:18:10.738771Z",
     "iopub.status.idle": "2023-07-26T05:18:10.848441Z",
     "shell.execute_reply": "2023-07-26T05:18:10.847972Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33808 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scale_arrow = s_ = 2\n",
    "scores[:,1] *= -1\n",
    "pcaUS.components_[1] *= -1 # flip the y-axis\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "ax.scatter(scores[:,0], scores[:,1])\n",
    "ax.set_xlabel('PC%d' % (i+1))\n",
    "ax.set_ylabel('PC%d' % (j+1))\n",
    "for k in range(pcaUS.components_.shape[1]):\n",
    "    ax.arrow(0, 0, s_*pcaUS.components_[i,k], s_*pcaUS.components_[j,k])\n",
    "    ax.text(s_*pcaUS.components_[i,k],\n",
    "            s_*pcaUS.components_[j,k],\n",
    "            USArrests.columns[k])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e380e98d",
   "metadata": {},
   "source": [
    "The standard deviations of the principal component scores are as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "34fdfe21",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.852220Z",
     "iopub.status.busy": "2023-07-26T05:18:10.851752Z",
     "iopub.status.idle": "2023-07-26T05:18:10.859190Z",
     "shell.execute_reply": "2023-07-26T05:18:10.857389Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.5908673 , 1.00496987, 0.6031915 , 0.4206774 ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scores.std(0, ddof=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c24b8ca2",
   "metadata": {},
   "source": [
    "The variance of each score can be extracted directly from the `pcaUS` object via\n",
    "the `explained_variance_` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "31b43c57",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.861366Z",
     "iopub.status.busy": "2023-07-26T05:18:10.861239Z",
     "iopub.status.idle": "2023-07-26T05:18:10.864540Z",
     "shell.execute_reply": "2023-07-26T05:18:10.863998Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.53085875, 1.00996444, 0.36383998, 0.17696948])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.explained_variance_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2dcf543",
   "metadata": {},
   "source": [
    "The proportion of variance explained by each principal \n",
    "component (PVE) is stored as `explained_variance_ratio_`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "68e47d3a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.867030Z",
     "iopub.status.busy": "2023-07-26T05:18:10.866863Z",
     "iopub.status.idle": "2023-07-26T05:18:10.869949Z",
     "shell.execute_reply": "2023-07-26T05:18:10.869455Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.62006039, 0.24744129, 0.0891408 , 0.04335752])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.explained_variance_ratio_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "831b19b7",
   "metadata": {},
   "source": [
    "We see that the first principal component explains 62.0% of the\n",
    "variance in the data, the next principal component explains 24.7%\n",
    "of the variance, and so forth.\n",
    "We can plot the PVE explained by each component, as well as the cumulative PVE. We first\n",
    "plot the proportion of variance explained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e87fe198",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:10.872333Z",
     "iopub.status.busy": "2023-07-26T05:18:10.872213Z",
     "iopub.status.idle": "2023-07-26T05:18:11.030288Z",
     "shell.execute_reply": "2023-07-26T05:18:11.029440Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "fig, axes = plt.subplots(1, 2, figsize=(15, 6))\n",
    "ticks = np.arange(pcaUS.n_components_)+1\n",
    "ax = axes[0]\n",
    "ax.plot(ticks,\n",
    "        pcaUS.explained_variance_ratio_,\n",
    "        marker='o')\n",
    "ax.set_xlabel('Principal Component');\n",
    "ax.set_ylabel('Proportion of Variance Explained')\n",
    "ax.set_ylim([0,1])\n",
    "ax.set_xticks(ticks)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ed96395",
   "metadata": {},
   "source": [
    "Notice the use of `%%capture`, which suppresses the displaying of the partially completed figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "409fb0c6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.033250Z",
     "iopub.status.busy": "2023-07-26T05:18:11.033097Z",
     "iopub.status.idle": "2023-07-26T05:18:11.153804Z",
     "shell.execute_reply": "2023-07-26T05:18:11.153245Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 50844 characters]",
      "text/plain": [
       "<Figure size 1500x600 with 2 Axes>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax = axes[1]\n",
    "ax.plot(ticks,\n",
    "        pcaUS.explained_variance_ratio_.cumsum(),\n",
    "        marker='o')\n",
    "ax.set_xlabel('Principal Component')\n",
    "ax.set_ylabel('Cumulative Proportion of Variance Explained')\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_xticks(ticks)\n",
    "fig\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "496fb6be",
   "metadata": {},
   "source": [
    "The result is similar to that shown in Figure 12.3.  Note\n",
    "that the method `cumsum()`   computes the cumulative sum of\n",
    "the elements of a numeric vector. For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e563e41b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.157597Z",
     "iopub.status.busy": "2023-07-26T05:18:11.157216Z",
     "iopub.status.idle": "2023-07-26T05:18:11.162432Z",
     "shell.execute_reply": "2023-07-26T05:18:11.162026Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  3, 11,  8])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([1,2,8,-3])\n",
    "np.cumsum(a)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6794b1a3",
   "metadata": {},
   "source": [
    "## Matrix Completion\n",
    "\n",
    "We now re-create the analysis carried out on the `USArrests` data in\n",
    "Section 12.3.\n",
    "\n",
    "We saw  in Section 12.2.2  that solving the optimization\n",
    "problem (12.6)   on a centered data matrix $\\bf X$ is\n",
    "equivalent to computing the first $M$ principal\n",
    "components of the data.  We use our scaled\n",
    "and centered `USArrests` data as $\\bf X$ below. The *singular value decomposition* \n",
    "(SVD)  is a general algorithm for solving\n",
    "(12.6). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "f83ad0bc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.164364Z",
     "iopub.status.busy": "2023-07-26T05:18:11.164225Z",
     "iopub.status.idle": "2023-07-26T05:18:11.167031Z",
     "shell.execute_reply": "2023-07-26T05:18:11.166724Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((50, 4), (4,), (4, 4))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = USArrests_scaled\n",
    "U, D, V = np.linalg.svd(X, full_matrices=False)\n",
    "U.shape, D.shape, V.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c71e57",
   "metadata": {},
   "source": [
    "The `np.linalg.svd()` function returns three components, `U`, `D` and `V`. The matrix `V` is equivalent to the\n",
    "loading matrix from principal components (up to an unimportant sign flip). Using the `full_matrices=False` option ensures that\n",
    "for a tall matrix the shape of `U` is the same as the shape of `X`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "cb9bdc46",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.169269Z",
     "iopub.status.busy": "2023-07-26T05:18:11.169096Z",
     "iopub.status.idle": "2023-07-26T05:18:11.172167Z",
     "shell.execute_reply": "2023-07-26T05:18:11.171624Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.53589947, -0.58318363, -0.27819087, -0.54343209],\n",
       "       [ 0.41818087,  0.1879856 , -0.87280619, -0.16731864],\n",
       "       [-0.34123273, -0.26814843, -0.37801579,  0.81777791],\n",
       "       [ 0.6492278 , -0.74340748,  0.13387773,  0.08902432]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f23c101e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.174404Z",
     "iopub.status.busy": "2023-07-26T05:18:11.174265Z",
     "iopub.status.idle": "2023-07-26T05:18:11.177096Z",
     "shell.execute_reply": "2023-07-26T05:18:11.176712Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.53589947,  0.58318363,  0.27819087,  0.54343209],\n",
       "       [-0.41818087, -0.1879856 ,  0.87280619,  0.16731864],\n",
       "       [-0.34123273, -0.26814843, -0.37801579,  0.81777791],\n",
       "       [ 0.6492278 , -0.74340748,  0.13387773,  0.08902432]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcaUS.components_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5d9e0ae",
   "metadata": {},
   "source": [
    "The matrix `U` corresponds to a  *standardized* version of the PCA score matrix (each column standardized to have sum-of-squares one). If we multiply each column of `U` by the corresponding element  of `D`, we recover the PCA scores exactly (up to a meaningless sign flip)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4cc49622",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.179219Z",
     "iopub.status.busy": "2023-07-26T05:18:11.179054Z",
     "iopub.status.idle": "2023-07-26T05:18:11.182070Z",
     "shell.execute_reply": "2023-07-26T05:18:11.181699Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.98556588,  1.13339238, -0.44426879,  0.15626714],\n",
       "       [-1.95013775,  1.07321326,  2.04000333, -0.43858344],\n",
       "       [-1.76316354, -0.74595678,  0.05478082, -0.83465292]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(U * D[None,:])[:3]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c96c9fe1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.184051Z",
     "iopub.status.busy": "2023-07-26T05:18:11.183934Z",
     "iopub.status.idle": "2023-07-26T05:18:11.186945Z",
     "shell.execute_reply": "2023-07-26T05:18:11.186473Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.98556588, -1.13339238, -0.44426879,  0.15626714],\n",
       "       [ 1.95013775, -1.07321326,  2.04000333, -0.43858344],\n",
       "       [ 1.76316354,  0.74595678,  0.05478082, -0.83465292]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scores[:3]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b7002cb",
   "metadata": {},
   "source": [
    "While it would be possible to carry out this lab using the `PCA()` estimator,\n",
    "here we use the `np.linalg.svd()` function in order to illustrate its use.\n",
    "\n",
    "We now omit 20 entries in the $50\\times 4$ data matrix at random. We do so\n",
    "by first selecting 20 rows (states) at random, and then selecting one\n",
    "of the four entries in each row at random. This ensures that every row has\n",
    "at least three observed values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "574409d6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.189376Z",
     "iopub.status.busy": "2023-07-26T05:18:11.189199Z",
     "iopub.status.idle": "2023-07-26T05:18:11.192254Z",
     "shell.execute_reply": "2023-07-26T05:18:11.191873Z"
    }
   },
   "outputs": [],
   "source": [
    "n_omit = 20\n",
    "np.random.seed(15)\n",
    "r_idx = np.random.choice(np.arange(X.shape[0]),\n",
    "                         n_omit,\n",
    "                         replace=False)\n",
    "c_idx = np.random.choice(np.arange(X.shape[1]),\n",
    "                         n_omit,\n",
    "                         replace=True)\n",
    "Xna = X.copy()\n",
    "Xna[r_idx, c_idx] = np.nan\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd4a5fdf",
   "metadata": {},
   "source": [
    "Here the array `r_idx`\n",
    "contains 20 integers from 0 to 49; this represents the states (rows of `X`) that are selected to contain missing values. And `c_idx` contains\n",
    "20 integers from 0 to 3, representing the features (columns in `X`) that contain the missing values for each of the selected states.\n",
    "\n",
    "We now write some code to implement Algorithm 12.1. \n",
    "We first write a  function that takes in a matrix, and returns an approximation to the matrix using the `svd()` function.\n",
    "This will be needed in Step 2 of Algorithm 12.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "89f190ae",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.194569Z",
     "iopub.status.busy": "2023-07-26T05:18:11.194458Z",
     "iopub.status.idle": "2023-07-26T05:18:11.197204Z",
     "shell.execute_reply": "2023-07-26T05:18:11.196762Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "def low_rank(X, M=1):\n",
    "    U, D, V = np.linalg.svd(X)\n",
    "    L = U[:,:M] * D[None,:M]\n",
    "    return L.dot(V[:M])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a04129d0",
   "metadata": {},
   "source": [
    "To conduct Step 1 of the algorithm, we initialize `Xhat` --- this is $\\tilde{\\bf X}$ in Algorithm 12.1 ---  by replacing\n",
    "the missing values with the column means of the non-missing entries. These are stored in\n",
    "`Xbar` below after running `np.nanmean()` over the row axis.\n",
    "We make a copy so that when we assign values to `Xhat` below we do not also overwrite the\n",
    "values in `Xna`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "322f339c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.199355Z",
     "iopub.status.busy": "2023-07-26T05:18:11.199208Z",
     "iopub.status.idle": "2023-07-26T05:18:11.201784Z",
     "shell.execute_reply": "2023-07-26T05:18:11.201348Z"
    }
   },
   "outputs": [],
   "source": [
    "Xhat = Xna.copy()\n",
    "Xbar = np.nanmean(Xhat, axis=0)\n",
    "Xhat[r_idx, c_idx] = Xbar[c_idx]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "357041e8",
   "metadata": {},
   "source": [
    "Before we begin Step 2, we set ourselves up to measure the progress of our\n",
    "iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7e106d1a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.204020Z",
     "iopub.status.busy": "2023-07-26T05:18:11.203890Z",
     "iopub.status.idle": "2023-07-26T05:18:11.206286Z",
     "shell.execute_reply": "2023-07-26T05:18:11.205915Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "thresh = 1e-7\n",
    "rel_err = 1\n",
    "count = 0\n",
    "ismiss = np.isnan(Xna)\n",
    "mssold = np.mean(Xhat[~ismiss]**2)\n",
    "mss0 = np.mean(Xna[~ismiss]**2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105ac73f",
   "metadata": {},
   "source": [
    "Here  `ismiss` is a logical matrix with the same dimensions as `Xna`;\n",
    "a given element is `True` if the corresponding matrix element is missing. The notation `~ismiss` negates this boolean vector. This is useful\n",
    "because it allows us to access both the missing and non-missing entries. We store the mean of the squared non-missing elements in `mss0`.\n",
    "We store the mean squared error  of the non-missing elements  of the old version of `Xhat` in `mssold` (which currently\n",
    "agrees with `mss0`). We plan to store the mean squared error of the non-missing elements of the current version of `Xhat` in `mss`, and will then\n",
    "iterate Step 2 of  Algorithm 12.1  until the *relative error*, defined as\n",
    "`(mssold - mss) / mss0`, falls below `thresh = 1e-7`.\n",
    " {Algorithm 12.1 tells us to iterate Step 2 until (12.14) is no longer decreasing. Determining whether (12.14)  is decreasing requires us only to keep track of `mssold - mss`. However, in practice, we keep track of `(mssold - mss) / mss0` instead: this makes it so that the number of iterations required for Algorithm 12.1 to converge does not depend on whether we multiplied the raw data $\\bf X$ by a constant factor.}\n",
    "\n",
    "In Step 2(a) of Algorithm 12.1, we  approximate `Xhat` using `low_rank()`; we call this `Xapp`. In Step 2(b), we  use `Xapp`  to update the estimates for elements in `Xhat` that are missing in `Xna`. Finally, in Step 2(c), we compute the relative error. These three steps are contained in the following `while` loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "7cb05ce4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.208518Z",
     "iopub.status.busy": "2023-07-26T05:18:11.208345Z",
     "iopub.status.idle": "2023-07-26T05:18:11.212476Z",
     "shell.execute_reply": "2023-07-26T05:18:11.211847Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 1, MSS:0.395, Rel.Err 5.99e-01\n",
      "Iteration: 2, MSS:0.382, Rel.Err 1.33e-02\n",
      "Iteration: 3, MSS:0.381, Rel.Err 1.44e-03\n",
      "Iteration: 4, MSS:0.381, Rel.Err 1.79e-04\n",
      "Iteration: 5, MSS:0.381, Rel.Err 2.58e-05\n",
      "Iteration: 6, MSS:0.381, Rel.Err 4.22e-06\n",
      "Iteration: 7, MSS:0.381, Rel.Err 7.65e-07\n",
      "Iteration: 8, MSS:0.381, Rel.Err 1.48e-07\n",
      "Iteration: 9, MSS:0.381, Rel.Err 2.95e-08\n"
     ]
    }
   ],
   "source": [
    "while rel_err > thresh:\n",
    "    count += 1\n",
    "    # Step 2(a)\n",
    "    Xapp = low_rank(Xhat, M=1)\n",
    "    # Step 2(b)\n",
    "    Xhat[ismiss] = Xapp[ismiss]\n",
    "    # Step 2(c)\n",
    "    mss = np.mean(((Xna - Xapp)[~ismiss])**2)\n",
    "    rel_err = (mssold - mss) / mss0\n",
    "    mssold = mss\n",
    "    print(\"Iteration: {0}, MSS:{1:.3f}, Rel.Err {2:.2e}\"\n",
    "          .format(count, mss, rel_err))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cd8edf2",
   "metadata": {},
   "source": [
    "We see that after eight iterations, the relative error has fallen below `thresh = 1e-7`, and so the algorithm terminates. When this happens, the mean squared error of the non-missing elements equals 0.381.\n",
    "\n",
    "Finally, we compute the correlation between the 20 imputed values\n",
    "and the actual values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "6f245188",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.214882Z",
     "iopub.status.busy": "2023-07-26T05:18:11.214746Z",
     "iopub.status.idle": "2023-07-26T05:18:11.217935Z",
     "shell.execute_reply": "2023-07-26T05:18:11.217596Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.711356743429736"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.corrcoef(Xapp[ismiss], X[ismiss])[0,1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f22444ff",
   "metadata": {},
   "source": [
    "In this lab, we implemented  Algorithm 12.1  ourselves for didactic purposes. However, a reader who wishes to apply matrix completion to their data might look to more specialized `Python`  implementations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7da441a1",
   "metadata": {},
   "source": [
    "## Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef4903d6",
   "metadata": {},
   "source": [
    "### $K$-Means Clustering\n",
    "\n",
    "The estimator `sklearn.cluster.KMeans()`  performs $K$-means clustering in\n",
    "`Python`.  We begin with a simple simulated example in which there\n",
    "truly are two clusters in the data: the first 25 observations have a\n",
    "mean shift relative to the next 25 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "345fb41e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.220117Z",
     "iopub.status.busy": "2023-07-26T05:18:11.219983Z",
     "iopub.status.idle": "2023-07-26T05:18:11.222355Z",
     "shell.execute_reply": "2023-07-26T05:18:11.221948Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "np.random.seed(0);\n",
    "X = np.random.standard_normal((50,2));\n",
    "X[:25,0] += 3;\n",
    "X[:25,1] -= 4;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ed3587f",
   "metadata": {},
   "source": [
    "We now perform $K$-means clustering with $K=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "3a8c21a2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.224388Z",
     "iopub.status.busy": "2023-07-26T05:18:11.224258Z",
     "iopub.status.idle": "2023-07-26T05:18:11.262343Z",
     "shell.execute_reply": "2023-07-26T05:18:11.261799Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=2,\n",
    "                random_state=2,\n",
    "                n_init=20).fit(X)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df4fb0eb",
   "metadata": {},
   "source": [
    "We specify `random_state` to make the results reproducible.  The cluster assignments of the 50 observations are contained in `kmeans.labels_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "e3e35b5d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.265599Z",
     "iopub.status.busy": "2023-07-26T05:18:11.265438Z",
     "iopub.status.idle": "2023-07-26T05:18:11.268983Z",
     "shell.execute_reply": "2023-07-26T05:18:11.268283Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,\n",
       "       1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
       "       0, 0, 0, 0, 0, 0], dtype=int32)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmeans.labels_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6632fa4c",
   "metadata": {},
   "source": [
    "The $K$-means clustering perfectly separated the observations into two\n",
    "clusters even though we did not supply any group information to\n",
    "`KMeans()`. We can plot the data, with each observation\n",
    "colored according to its cluster assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "d928650a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.271491Z",
     "iopub.status.busy": "2023-07-26T05:18:11.271353Z",
     "iopub.status.idle": "2023-07-26T05:18:11.409914Z",
     "shell.execute_reply": "2023-07-26T05:18:11.409365Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35636 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8,8))\n",
    "ax.scatter(X[:,0], X[:,1], c=kmeans.labels_)\n",
    "ax.set_title(\"K-Means Clustering Results with K=2\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f6a6d01",
   "metadata": {},
   "source": [
    "Here the observations can be easily plotted because they are\n",
    "two-dimensional. If there were more than two variables then we could\n",
    "instead perform PCA and plot the first two principal component score\n",
    "vectors to represent the clusters.\n",
    "\n",
    "In this example,   we knew that there really\n",
    "were two clusters because we generated the data. However, for real\n",
    "data, we do not know the true number of clusters, nor whether they  exist in any precise way. We could\n",
    "instead have performed $K$-means clustering on this example with\n",
    "$K=3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "92e5175c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.412487Z",
     "iopub.status.busy": "2023-07-26T05:18:11.412250Z",
     "iopub.status.idle": "2023-07-26T05:18:11.548124Z",
     "shell.execute_reply": "2023-07-26T05:18:11.547642Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36560 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kmeans = KMeans(n_clusters=3,\n",
    "                random_state=3,\n",
    "                n_init=20).fit(X)\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "ax.scatter(X[:,0], X[:,1], c=kmeans.labels_)\n",
    "ax.set_title(\"K-Means Clustering Results with K=3\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82317a30",
   "metadata": {},
   "source": [
    "When $K=3$, $K$-means clustering  splits up the two clusters.\n",
    "We have used the `n_init` argument to run the $K$-means with 20 \n",
    "initial cluster assignments (the default is 10). If a\n",
    "value of `n_init` greater than one is used, then $K$-means\n",
    "clustering will be performed using multiple random assignments in\n",
    "Step 1 of  Algorithm 12.2, and the `KMeans()` \n",
    "function will report only the best results. Here we compare using\n",
    "`n_init=1` to `n_init=20`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "4911ecc7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.551298Z",
     "iopub.status.busy": "2023-07-26T05:18:11.551122Z",
     "iopub.status.idle": "2023-07-26T05:18:11.570392Z",
     "shell.execute_reply": "2023-07-26T05:18:11.569843Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(78.06117325882116, 75.0350825910044)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmeans1 = KMeans(n_clusters=3,\n",
    "                random_state=3,\n",
    "                n_init=1).fit(X)\n",
    "kmeans20 = KMeans(n_clusters=3,\n",
    "                  random_state=3,\n",
    "                  n_init=20).fit(X);\n",
    "kmeans1.inertia_, kmeans20.inertia_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96e6df1f",
   "metadata": {},
   "source": [
    "Note that `kmeans.inertia_` is the total within-cluster sum\n",
    "of squares, which we seek to minimize by performing $K$-means\n",
    "clustering (12.17). \n",
    "\n",
    "We *strongly* recommend always running $K$-means clustering with\n",
    "a large value of `n_init`, such as 20 or 50, since otherwise an\n",
    "undesirable local optimum may be obtained.\n",
    "\n",
    "When performing $K$-means clustering, in addition to using multiple\n",
    "initial cluster assignments, it is also important to set a random seed\n",
    "using the `random_state` argument to `KMeans()`. This way, the initial\n",
    "cluster assignments in Step 1 can be replicated, and the $K$-means\n",
    "output will be fully reproducible."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "904ea858",
   "metadata": {},
   "source": [
    "### Hierarchical Clustering\n",
    "\n",
    "The `AgglomerativeClustering()`  class from\n",
    "the `sklearn.clustering` package implements hierarchical clustering.\n",
    "As its\n",
    "name is long, we use the short hand `HClust` for *hierarchical clustering*. Note that this will not change the return type\n",
    "when using this method, so instances will still be of class `AgglomerativeClustering`.\n",
    "In the following example we use the data from the previous lab to plot the hierarchical clustering\n",
    "dendrogram using complete, single, and average linkage clustering\n",
    "with Euclidean distance as the dissimilarity measure.  We begin by\n",
    "clustering observations using complete linkage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "1b42a700",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.573065Z",
     "iopub.status.busy": "2023-07-26T05:18:11.572885Z",
     "iopub.status.idle": "2023-07-26T05:18:11.577747Z",
     "shell.execute_reply": "2023-07-26T05:18:11.577317Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-2 {color: black;background-color: white;}#sk-container-id-2 pre{padding: 0;}#sk-container-id-2 div.sk-toggleable {background-color: white;}#sk-container-id-2 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-2 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-2 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-2 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-2 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-2 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-2 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-2 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-2 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-2 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-2 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-2 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-2 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-2 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-2 div.sk-item {position: relative;z-index: 1;}#sk-container-id-2 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-2 div.sk-item::before, #sk-container-id-2 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-2 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-2 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-2 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-2 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-2 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-2 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-2 div.sk-label-container {text-align: center;}#sk-container-id-2 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-2 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-2\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>AgglomerativeClustering(distance_threshold=0, linkage=&#x27;complete&#x27;,\n",
       "                        n_clusters=None)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" checked><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">AgglomerativeClustering</label><div class=\"sk-toggleable__content\"><pre>AgglomerativeClustering(distance_threshold=0, linkage=&#x27;complete&#x27;,\n",
       "                        n_clusters=None)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "AgglomerativeClustering(distance_threshold=0, linkage='complete',\n",
       "                        n_clusters=None)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HClust = AgglomerativeClustering\n",
    "hc_comp = HClust(distance_threshold=0,\n",
    "                 n_clusters=None,\n",
    "                 linkage='complete')\n",
    "hc_comp.fit(X)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22199f4d",
   "metadata": {},
   "source": [
    "This computes the entire dendrogram.\n",
    "We could just as easily perform hierarchical clustering with average or single linkage instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "50ef7eea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.580226Z",
     "iopub.status.busy": "2023-07-26T05:18:11.580003Z",
     "iopub.status.idle": "2023-07-26T05:18:11.583835Z",
     "shell.execute_reply": "2023-07-26T05:18:11.583381Z"
    }
   },
   "outputs": [],
   "source": [
    "hc_avg = HClust(distance_threshold=0,\n",
    "                n_clusters=None,\n",
    "                linkage='average');\n",
    "hc_avg.fit(X)\n",
    "hc_sing = HClust(distance_threshold=0,\n",
    "                 n_clusters=None,\n",
    "                 linkage='single');\n",
    "hc_sing.fit(X);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e8db256",
   "metadata": {},
   "source": [
    "To use a precomputed distance matrix, we provide an additional\n",
    "argument `metric=\"precomputed\"`. In the code below, the first four lines computes the $50\\times 50$ pairwise-distance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "bf7a2408",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.586413Z",
     "iopub.status.busy": "2023-07-26T05:18:11.586258Z",
     "iopub.status.idle": "2023-07-26T05:18:11.591990Z",
     "shell.execute_reply": "2023-07-26T05:18:11.591553Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-3 {color: black;background-color: white;}#sk-container-id-3 pre{padding: 0;}#sk-container-id-3 div.sk-toggleable {background-color: white;}#sk-container-id-3 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-3 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-3 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-3 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-3 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-3 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-3 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-3 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-3 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-3 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-3 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-3 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-3 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-3 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-3 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-3 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-3 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-3 div.sk-item {position: relative;z-index: 1;}#sk-container-id-3 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-3 div.sk-item::before, #sk-container-id-3 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-3 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-3 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-3 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-3 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-3 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-3 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-3 div.sk-label-container {text-align: center;}#sk-container-id-3 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-3 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-3\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>AgglomerativeClustering(distance_threshold=0, linkage=&#x27;single&#x27;,\n",
       "                        metric=&#x27;precomputed&#x27;, n_clusters=None)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-3\" type=\"checkbox\" checked><label for=\"sk-estimator-id-3\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">AgglomerativeClustering</label><div class=\"sk-toggleable__content\"><pre>AgglomerativeClustering(distance_threshold=0, linkage=&#x27;single&#x27;,\n",
       "                        metric=&#x27;precomputed&#x27;, n_clusters=None)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "AgglomerativeClustering(distance_threshold=0, linkage='single',\n",
       "                        metric='precomputed', n_clusters=None)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = np.zeros((X.shape[0], X.shape[0]));\n",
    "for i in range(X.shape[0]):\n",
    "    x_ = np.multiply.outer(np.ones(X.shape[0]), X[i])\n",
    "    D[i] = np.sqrt(np.sum((X - x_)**2, 1));\n",
    "hc_sing_pre = HClust(distance_threshold=0,\n",
    "                     n_clusters=None,\n",
    "                     metric='precomputed',\n",
    "                     linkage='single')\n",
    "hc_sing_pre.fit(D)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cd7fb7c",
   "metadata": {},
   "source": [
    "We use\n",
    "`dendrogram()` from `scipy.cluster.hierarchy` to plot the dendrogram. However,\n",
    "`dendrogram()` expects a so-called *linkage-matrix representation*\n",
    "of the clustering, which is not provided by `AgglomerativeClustering()`,\n",
    "but can be computed. The function `compute_linkage()` in the\n",
    "`ISLP.cluster` package is provided for this purpose.\n",
    "\n",
    "We can now plot the dendrograms. The numbers at the bottom of the plot\n",
    "identify each observation. The `dendrogram()` function has a default method to\n",
    "color different branches of the tree that suggests a pre-defined cut of the tree at a particular depth.\n",
    "We prefer to overwrite this default by setting this threshold to be infinite. Since we want this behavior for many dendrograms, we store these values in a dictionary `cargs` and pass this as keyword arguments using the notation `**cargs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "a118c0ab",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.595666Z",
     "iopub.status.busy": "2023-07-26T05:18:11.595400Z",
     "iopub.status.idle": "2023-07-26T05:18:11.876074Z",
     "shell.execute_reply": "2023-07-26T05:18:11.875555Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15716 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cargs = {'color_threshold':-np.inf,\n",
    "         'above_threshold_color':'black'}\n",
    "linkage_comp = compute_linkage(hc_comp)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "dendrogram(linkage_comp,\n",
    "           ax=ax,\n",
    "           **cargs);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a325a001",
   "metadata": {},
   "source": [
    "We may want to color branches of the tree above\n",
    "and below a cut-threshold differently. This can be achieved\n",
    "by changing the `color_threshold`. Let’s cut the tree at a height of 4,\n",
    "coloring links that merge above 4 in black."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "b1ff41c0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:11.878711Z",
     "iopub.status.busy": "2023-07-26T05:18:11.878497Z",
     "iopub.status.idle": "2023-07-26T05:18:12.156685Z",
     "shell.execute_reply": "2023-07-26T05:18:12.156217Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17252 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "dendrogram(linkage_comp,\n",
    "           ax=ax,\n",
    "           color_threshold=4,\n",
    "           above_threshold_color='black');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e774708c",
   "metadata": {},
   "source": [
    "To determine the cluster labels for each observation associated with a\n",
    "given cut of the dendrogram, we can use the `cut_tree()` \n",
    "function from `scipy.cluster.hierarchy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "c2752a96",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.159239Z",
     "iopub.status.busy": "2023-07-26T05:18:12.159059Z",
     "iopub.status.idle": "2023-07-26T05:18:12.163542Z",
     "shell.execute_reply": "2023-07-26T05:18:12.163127Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 2, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 2,\n",
       "        0, 2, 2, 3, 2, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2, 3, 3,\n",
       "        3, 3, 3, 3, 3, 3]])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cut_tree(linkage_comp, n_clusters=4).T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d501d220",
   "metadata": {},
   "source": [
    "This can also be achieved by providing an argument `n_clusters`\n",
    "to `HClust()`; however each cut would require recomputing\n",
    "the clustering. Similarly, trees may be cut by distance threshold\n",
    "with an argument of `distance_threshold` to `HClust()`\n",
    "or `height` to `cut_tree()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "1407f7a4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.165971Z",
     "iopub.status.busy": "2023-07-26T05:18:12.165813Z",
     "iopub.status.idle": "2023-07-26T05:18:12.170809Z",
     "shell.execute_reply": "2023-07-26T05:18:12.170332Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [1],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [0],\n",
       "       [1],\n",
       "       [0],\n",
       "       [1],\n",
       "       [1],\n",
       "       [2],\n",
       "       [1],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [1],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [1],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [1],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2],\n",
       "       [2]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cut_tree(linkage_comp, height=5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469bbad6",
   "metadata": {},
   "source": [
    "To scale the variables before performing hierarchical clustering of\n",
    "the observations, we use `StandardScaler()`  as in our PCA example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "2d74f224",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.173359Z",
     "iopub.status.busy": "2023-07-26T05:18:12.173198Z",
     "iopub.status.idle": "2023-07-26T05:18:12.457460Z",
     "shell.execute_reply": "2023-07-26T05:18:12.457026Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23412 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scaler = StandardScaler()\n",
    "X_scale = scaler.fit_transform(X)\n",
    "hc_comp_scale = HClust(distance_threshold=0,\n",
    "                       n_clusters=None,\n",
    "                       linkage='complete').fit(X_scale)\n",
    "linkage_comp_scale = compute_linkage(hc_comp_scale)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "dendrogram(linkage_comp_scale, ax=ax, **cargs)\n",
    "ax.set_title(\"Hierarchical Clustering with Scaled Features\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3caaa1a5",
   "metadata": {},
   "source": [
    "Correlation-based distances between observations can be used for\n",
    "clustering. The correlation between two observations measures the\n",
    "similarity of their feature values. {Suppose each observation has\n",
    "  $p$ features, each a single numerical value. We measure the\n",
    "  similarity of two such observations by computing the\n",
    "  correlation of these $p$ pairs of numbers.}\n",
    "With $n$ observations, the $n\\times n$ correlation matrix can then be used as a similarity (or affinity) matrix, i.e. so that one minus the correlation matrix is the dissimilarity matrix used for clustering.\n",
    "\n",
    "Note that using correlation only makes sense for\n",
    "data with at least three features since the absolute correlation\n",
    "between any two observations with measurements on two features is\n",
    "always one. Hence, we will cluster a three-dimensional data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "b7f7da12",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.461207Z",
     "iopub.status.busy": "2023-07-26T05:18:12.460949Z",
     "iopub.status.idle": "2023-07-26T05:18:12.680696Z",
     "shell.execute_reply": "2023-07-26T05:18:12.680192Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 37780 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.random.standard_normal((30, 3))\n",
    "corD = 1 - np.corrcoef(X)\n",
    "hc_cor = HClust(linkage='complete',\n",
    "                distance_threshold=0,\n",
    "                n_clusters=None,\n",
    "                metric='precomputed')\n",
    "hc_cor.fit(corD)\n",
    "linkage_cor = compute_linkage(hc_cor)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "dendrogram(linkage_cor, ax=ax, **cargs)\n",
    "ax.set_title(\"Complete Linkage with Correlation-Based Dissimilarity\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b00804e",
   "metadata": {},
   "source": [
    "## NCI60 Data Example\n",
    "Unsupervised techniques are often used in the analysis of genomic\n",
    "data. In particular, PCA and hierarchical clustering are popular\n",
    "tools.  We illustrate these techniques on the `NCI60`  cancer cell line\n",
    "microarray data, which consists of 6830 gene expression\n",
    "measurements on 64 cancer cell lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "b94424fc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.683099Z",
     "iopub.status.busy": "2023-07-26T05:18:12.682918Z",
     "iopub.status.idle": "2023-07-26T05:18:12.688821Z",
     "shell.execute_reply": "2023-07-26T05:18:12.688474Z"
    }
   },
   "outputs": [],
   "source": [
    "NCI60 = load_data('NCI60')\n",
    "nci_labs = NCI60['labels']\n",
    "nci_data = NCI60['data']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7313fea6",
   "metadata": {},
   "source": [
    "Each cell line is labeled with a cancer type. We do not make use of\n",
    "the cancer types in performing PCA and clustering, as these are\n",
    "unsupervised techniques. But after performing PCA and clustering, we\n",
    "will check to see the extent to which these cancer types agree with\n",
    "the results of these unsupervised techniques.\n",
    "\n",
    "The data has 64 rows and 6830 columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "cea54566",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.691075Z",
     "iopub.status.busy": "2023-07-26T05:18:12.690934Z",
     "iopub.status.idle": "2023-07-26T05:18:12.693705Z",
     "shell.execute_reply": "2023-07-26T05:18:12.693319Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64, 6830)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nci_data.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ebf27f1",
   "metadata": {},
   "source": [
    "We begin by examining the cancer types for the cell lines.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "4dac41bb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.695954Z",
     "iopub.status.busy": "2023-07-26T05:18:12.695815Z",
     "iopub.status.idle": "2023-07-26T05:18:12.700521Z",
     "shell.execute_reply": "2023-07-26T05:18:12.700074Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "label      \n",
       "NSCLC          9\n",
       "RENAL          9\n",
       "MELANOMA       8\n",
       "BREAST         7\n",
       "COLON          7\n",
       "LEUKEMIA       6\n",
       "OVARIAN        6\n",
       "CNS            5\n",
       "PROSTATE       2\n",
       "K562A-repro    1\n",
       "K562B-repro    1\n",
       "MCF7A-repro    1\n",
       "MCF7D-repro    1\n",
       "UNKNOWN        1\n",
       "dtype: int64"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nci_labs.value_counts()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f16f96a",
   "metadata": {},
   "source": [
    "### PCA on the NCI60 Data\n",
    "\n",
    "We first perform PCA on the data after scaling the variables (genes)\n",
    "to have standard deviation one, although here one could reasonably argue\n",
    "that it is better not to scale the genes as they are measured in the same units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "d8ebadd6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.702785Z",
     "iopub.status.busy": "2023-07-26T05:18:12.702646Z",
     "iopub.status.idle": "2023-07-26T05:18:12.866311Z",
     "shell.execute_reply": "2023-07-26T05:18:12.856469Z"
    }
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "nci_scaled = scaler.fit_transform(nci_data)\n",
    "nci_pca = PCA()\n",
    "nci_scores = nci_pca.fit_transform(nci_scaled)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a8c9932",
   "metadata": {},
   "source": [
    "We now plot the first few principal component score vectors, in order\n",
    "to visualize the data. The observations (cell lines) corresponding to\n",
    "a given cancer type will be plotted in the same color, so that we can\n",
    "see to what extent the observations within a cancer type are similar\n",
    "to each other. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "63b5efe3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:12.882679Z",
     "iopub.status.busy": "2023-07-26T05:18:12.878318Z",
     "iopub.status.idle": "2023-07-26T05:18:13.099177Z",
     "shell.execute_reply": "2023-07-26T05:18:13.098712Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 68628 characters]",
      "text/plain": [
       "<Figure size 1500x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cancer_types = list(np.unique(nci_labs))\n",
    "nci_groups = np.array([cancer_types.index(lab)\n",
    "                       for lab in nci_labs.values])\n",
    "fig, axes = plt.subplots(1, 2, figsize=(15,6))\n",
    "ax = axes[0]\n",
    "ax.scatter(nci_scores[:,0],\n",
    "           nci_scores[:,1],\n",
    "           c=nci_groups,\n",
    "           marker='o',\n",
    "           s=50)\n",
    "ax.set_xlabel('PC1'); ax.set_ylabel('PC2')\n",
    "ax = axes[1]\n",
    "ax.scatter(nci_scores[:,0],\n",
    "           nci_scores[:,2],\n",
    "           c=nci_groups,\n",
    "           marker='o',\n",
    "           s=50)\n",
    "ax.set_xlabel('PC1'); ax.set_ylabel('PC3');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87b33a9c",
   "metadata": {},
   "source": [
    "On the whole, cell lines corresponding to a single cancer type do tend to\n",
    "have similar values on the first few principal component score\n",
    "vectors. This indicates that cell lines from the same cancer type tend\n",
    "to have pretty similar gene expression levels.\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "890fdb54",
   "metadata": {},
   "source": [
    "We can also plot the percent variance\n",
    "explained by the principal components as well as the cumulative percent variance explained.\n",
    "This is similar to the plots we made earlier for the `USArrests` data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "e20c3cc1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:13.101683Z",
     "iopub.status.busy": "2023-07-26T05:18:13.101534Z",
     "iopub.status.idle": "2023-07-26T05:18:13.327522Z",
     "shell.execute_reply": "2023-07-26T05:18:13.327041Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47200 characters]",
      "text/plain": [
       "<Figure size 1500x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(15,6))\n",
    "ax = axes[0]\n",
    "ticks = np.arange(nci_pca.n_components_)+1\n",
    "ax.plot(ticks,\n",
    "        nci_pca.explained_variance_ratio_,\n",
    "        marker='o')\n",
    "ax.set_xlabel('Principal Component');\n",
    "ax.set_ylabel('PVE')\n",
    "ax = axes[1]\n",
    "ax.plot(ticks,\n",
    "        nci_pca.explained_variance_ratio_.cumsum(),\n",
    "        marker='o');\n",
    "ax.set_xlabel('Principal Component')\n",
    "ax.set_ylabel('Cumulative PVE');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb6eb4f8",
   "metadata": {},
   "source": [
    "We see that together, the first seven principal components explain\n",
    "around 40% of the variance in the data. This is not a huge amount\n",
    "of the variance. However, looking at the scree plot, we see that while\n",
    "each of the first seven principal components explain a substantial\n",
    "amount of variance, there is a marked decrease in the variance\n",
    "explained by further principal components. That is, there is an\n",
    "*elbow*  in the plot after approximately the seventh\n",
    "principal component.  This suggests that there may be little benefit\n",
    "to examining more than seven or so principal components (though even\n",
    "examining seven principal components may be difficult)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8962791f",
   "metadata": {},
   "source": [
    "### Clustering the Observations of the NCI60 Data\n",
    "\n",
    "We now perform hierarchical clustering of the cell lines in the `NCI60` data using\n",
    "complete, single, and   average linkage. Once again, the goal is to find out whether or not the observations cluster into distinct types of cancer. Euclidean\n",
    "distance is used as the dissimilarity measure. We first write a short\n",
    "function to  produce\n",
    "the three dendrograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "622de805",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:13.330190Z",
     "iopub.status.busy": "2023-07-26T05:18:13.330013Z",
     "iopub.status.idle": "2023-07-26T05:18:13.332964Z",
     "shell.execute_reply": "2023-07-26T05:18:13.332533Z"
    }
   },
   "outputs": [],
   "source": [
    "def plot_nci(linkage, ax, cut=-np.inf):\n",
    "    cargs = {'above_threshold_color':'black',\n",
    "             'color_threshold':cut}\n",
    "    hc = HClust(n_clusters=None,\n",
    "                distance_threshold=0,\n",
    "                linkage=linkage.lower()).fit(nci_scaled)\n",
    "    linkage_ = compute_linkage(hc)\n",
    "    dendrogram(linkage_,\n",
    "               ax=ax,\n",
    "               labels=np.asarray(nci_labs),\n",
    "               leaf_font_size=10,\n",
    "               **cargs)\n",
    "    ax.set_title('%s Linkage' % linkage)\n",
    "    return hc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cf5f836",
   "metadata": {},
   "source": [
    "Let’s  plot our results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "54d40449",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:13.336193Z",
     "iopub.status.busy": "2023-07-26T05:18:13.335972Z",
     "iopub.status.idle": "2023-07-26T05:18:14.886877Z",
     "shell.execute_reply": "2023-07-26T05:18:14.886190Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 171932 characters]",
      "text/plain": [
       "<Figure size 1500x3000 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=(15,30))      \n",
    "ax = axes[0]; hc_comp = plot_nci('Complete', ax)\n",
    "ax = axes[1]; hc_avg = plot_nci('Average', ax)\n",
    "ax = axes[2]; hc_sing = plot_nci('Single', ax)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1703afbb",
   "metadata": {},
   "source": [
    "We see that the\n",
    "choice of linkage certainly does affect the results\n",
    "obtained. Typically, single linkage will tend to yield *trailing*\n",
    "clusters: very large clusters onto which individual observations\n",
    "attach one-by-one. On the other hand, complete and average linkage\n",
    "tend to yield more balanced, attractive clusters. For this reason,\n",
    "complete and average linkage are generally preferred to single\n",
    "linkage.  Clearly cell lines within a single cancer type do tend to\n",
    "cluster together, although the clustering is not perfect. We will use\n",
    "complete linkage hierarchical clustering for the analysis that\n",
    "follows.\n",
    "\n",
    "We can cut the dendrogram at the height that will yield a particular\n",
    "number of clusters, say four:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "dc80afc8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:14.891531Z",
     "iopub.status.busy": "2023-07-26T05:18:14.890910Z",
     "iopub.status.idle": "2023-07-26T05:18:14.903469Z",
     "shell.execute_reply": "2023-07-26T05:18:14.903077Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Complete</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>label</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>BREAST</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CNS</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>COLON</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>K562A-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>K562B-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LEUKEMIA</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MCF7A-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MCF7D-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MELANOMA</th>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSCLC</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OVARIAN</th>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PROSTATE</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RENAL</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UNKNOWN</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Complete     0  1  2  3\n",
       "label                  \n",
       "BREAST       2  3  0  2\n",
       "CNS          3  2  0  0\n",
       "COLON        2  0  0  5\n",
       "K562A-repro  0  0  1  0\n",
       "K562B-repro  0  0  1  0\n",
       "LEUKEMIA     0  0  6  0\n",
       "MCF7A-repro  0  0  0  1\n",
       "MCF7D-repro  0  0  0  1\n",
       "MELANOMA     8  0  0  0\n",
       "NSCLC        8  1  0  0\n",
       "OVARIAN      6  0  0  0\n",
       "PROSTATE     2  0  0  0\n",
       "RENAL        8  1  0  0\n",
       "UNKNOWN      1  0  0  0"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linkage_comp = compute_linkage(hc_comp)\n",
    "comp_cut = cut_tree(linkage_comp, n_clusters=4).reshape(-1)\n",
    "pd.crosstab(nci_labs['label'],\n",
    "            pd.Series(comp_cut.reshape(-1), name='Complete'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8ecae69",
   "metadata": {},
   "source": [
    "There are some clear patterns. All the leukemia cell lines fall in\n",
    "one cluster, while the breast cancer cell lines are spread out over\n",
    "three different clusters.\n",
    "\n",
    "We can plot a cut on the dendrogram that produces these four clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "40ff59f9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:14.905992Z",
     "iopub.status.busy": "2023-07-26T05:18:14.905760Z",
     "iopub.status.idle": "2023-07-26T05:18:15.399751Z",
     "shell.execute_reply": "2023-07-26T05:18:15.399237Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 58392 characters]",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plot_nci('Complete', ax, cut=140)\n",
    "ax.axhline(140, c='r', linewidth=4);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c03368c",
   "metadata": {},
   "source": [
    "The `axhline()`  function draws a horizontal line  line on top of any\n",
    "existing set of axes. The argument `140` plots a horizontal\n",
    "line at height 140 on the dendrogram; this is a height that\n",
    "results in four distinct clusters. It is easy to verify that the\n",
    "resulting clusters are the same as the ones we obtained in\n",
    "`comp_cut`.\n",
    "\n",
    "We claimed earlier in Section 12.4.2 that\n",
    "$K$-means clustering and hierarchical clustering with the dendrogram\n",
    "cut to obtain the same number of clusters can yield very different\n",
    "results.  How do these `NCI60` hierarchical clustering results compare\n",
    "to what we get if we perform $K$-means clustering with $K=4$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "1587e83b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:15.401922Z",
     "iopub.status.busy": "2023-07-26T05:18:15.401768Z",
     "iopub.status.idle": "2023-07-26T05:18:17.931480Z",
     "shell.execute_reply": "2023-07-26T05:18:17.924183Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>K-means</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>HClust</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "K-means   0  1  2  3\n",
       "HClust              \n",
       "0        28  3  9  0\n",
       "1         7  0  0  0\n",
       "2         0  0  0  8\n",
       "3         0  9  0  0"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nci_kmeans = KMeans(n_clusters=4, \n",
    "                    random_state=0,\n",
    "                    n_init=20).fit(nci_scaled)\n",
    "pd.crosstab(pd.Series(comp_cut, name='HClust'),\n",
    "            pd.Series(nci_kmeans.labels_, name='K-means'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b49c4a73",
   "metadata": {},
   "source": [
    "We see that the four clusters obtained using hierarchical clustering\n",
    "and $K$-means clustering are somewhat different. First we note\n",
    "that the labels in the two clusterings are arbitrary. That is, swapping\n",
    "the identifier of the cluster does not\n",
    "change the clustering. We see here Cluster 3 in\n",
    "$K$-means clustering is identical to cluster 2 in hierarchical\n",
    "clustering. However, the other clusters differ: for instance,\n",
    "cluster 0 in $K$-means clustering contains a portion of the\n",
    "observations assigned to cluster 0 by hierarchical clustering, as well\n",
    "as all of the observations assigned to cluster 1 by hierarchical\n",
    "clustering.\n",
    "\n",
    "Rather than performing hierarchical clustering on the entire data\n",
    "matrix, we can also perform hierarchical clustering on the first few\n",
    "principal component score vectors, regarding these first few components\n",
    "as a less noisy version of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "b09ceeab",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:17.955945Z",
     "iopub.status.busy": "2023-07-26T05:18:17.954718Z",
     "iopub.status.idle": "2023-07-26T05:18:18.488264Z",
     "shell.execute_reply": "2023-07-26T05:18:18.487809Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Complete-PCA</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>label</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>BREAST</th>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CNS</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>COLON</th>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>K562A-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>K562B-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LEUKEMIA</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MCF7A-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MCF7D-repro</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MELANOMA</th>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSCLC</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OVARIAN</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PROSTATE</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RENAL</th>\n",
       "      <td>7</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>UNKNOWN</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Complete-PCA  0  1  2  3\n",
       "label                   \n",
       "BREAST        0  5  0  2\n",
       "CNS           2  3  0  0\n",
       "COLON         7  0  0  0\n",
       "K562A-repro   0  0  1  0\n",
       "K562B-repro   0  0  1  0\n",
       "LEUKEMIA      2  0  4  0\n",
       "MCF7A-repro   0  0  0  1\n",
       "MCF7D-repro   0  0  0  1\n",
       "MELANOMA      1  7  0  0\n",
       "NSCLC         8  1  0  0\n",
       "OVARIAN       5  1  0  0\n",
       "PROSTATE      2  0  0  0\n",
       "RENAL         7  2  0  0\n",
       "UNKNOWN       0  1  0  0"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 58816 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hc_pca = HClust(n_clusters=None,\n",
    "                distance_threshold=0,\n",
    "                linkage='complete'\n",
    "                ).fit(nci_scores[:,:5])\n",
    "linkage_pca = compute_linkage(hc_pca)\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "dendrogram(linkage_pca,\n",
    "           labels=np.asarray(nci_labs),\n",
    "           leaf_font_size=10,\n",
    "           ax=ax,\n",
    "           **cargs)\n",
    "ax.set_title(\"Hier. Clust. on First Five Score Vectors\")\n",
    "pca_labels = pd.Series(cut_tree(linkage_pca,\n",
    "                                n_clusters=4).reshape(-1),\n",
    "                       name='Complete-PCA')\n",
    "pd.crosstab(nci_labs['label'], pca_labels)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e1733b5",
   "metadata": {},
   "source": [
    "\n",
    " \n",
    " \n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
