Clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intermediate Machine Learning with scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering\n",
    "\n",
    "For an introduction to clustering algorithms, let us look at another visualization copied from the [scikit-learn documentation](http://scikit-learn.org/stable/modules/clustering.html). This is a visualization of several synthetic distributions in two dimensions to allow visualization.  We can  see that some distributions that do not have convex differentiation throw off some clutering algorithms.  Of course, given the low dimensionality, these examples are somewhat toys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%run src/plot_cluster_comparison.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher Dimensions\n",
    "\n",
    "For examples with greater dimensionality, we will look first at a couple similar datasets in the scikit-learn samples.  A small-ish set of 506 houses sold in Boston in 1978 are available.  A different, larger sample dataset of 20k houses in California, from about 20 years later, is also available.  The features contained in the two datasets have both overlaps (albeit with different names for corresponding features) and differences.\n",
    "\n",
    "Although normally clustering is used identify patterns where no *a priori* classes or target values are available, by using these datasets, we can test a certain intuition I have.  That is, my hunch would be that the way houses cluster along many feature dimensions will correspond to the clustering or distribution in the single price target.\n",
    "\n",
    "My hypothesis can well be wrong without meaning that clustering is meaningless.  The semantic sense of the clusters may simply be something other than price—perhaps some other implicit feature we have a good name for, but also perhaps something with no obvious English word corresponding to the essence of the classes.  Let us explore this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "boston = datasets.load_boston()\n",
    "print(boston.DESCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)\n",
    "boston_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distances in Parametric Space\n",
    "\n",
    "This is mostly for a later lesson, but many clustering algorithms are sensitive to the distances between points in a parameteric space in a unit-dependent way.  A later lesson on feature engineering will look at scaling data in more detail.  For some of the examples in this lesson, scaling will make little difference; however, for a few of them the absolute units of features becomes important.\n",
    "\n",
    "Therefore, in the below examples, you will see line like this in the cells that fit models:\n",
    "\n",
    "```python\n",
    "X = StandardScaler().fit_transform(the_data)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "X = StandardScaler().fit_transform(boston.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster into N Classes\n",
    "\n",
    "For the next few cells, we will use `SpectralClustering` to group our dataset into four synthetic classes.  Informally, we will call these classes \"low\", \"mid_low\", \"mid_high\", and \"high\" corresponding to our hypothesis that price is a good proxy for these classes.  From the documentation, the algorithm is described as follows:\n",
    "\n",
    "> `SpectralClustering` does a low-dimension embedding of the affinity matrix between samples, followed by a KMeans in the low dimensional space. It is especially efficient if the affinity matrix is sparse and the [`pyamg`](https://github.com/pyamg/pyamg) module is installed. `SpectralClustering` requires the number of clusters to be specified. It works well for a small number of clusters but is not advised when using many clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import cluster\n",
    "spectral = cluster.SpectralClustering(\n",
    "        n_clusters=4, eigen_solver='arpack',\n",
    "        affinity=\"nearest_neighbors\")\n",
    "\n",
    "spectral.fit(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After fitting the data, what we get is labels for each sample.  These are simply integers in `range(n_clusters)`, with no inherent semantic meaning.  I.e. 0, 1, 2, or 3 in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(spectral.labels_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we do the clustering, let us group the dataset by the cluster labels.  We will enhance the Pandas DataFrame that holds the data by adding columns for the derived category and the known price.  We also give hopefully more meaningful names to the clusters.  We only display some of the features to narrow discussion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us do minor Pandas manipulation to view clusters in an intuitive way\n",
    "boston_df['category'] = spectral.labels_\n",
    "boston_df['PRICE'] = boston.target\n",
    "boston_clusters = boston_df.groupby('category').mean().sort_values('PRICE')\n",
    "boston_clusters.index = ['low', 'mid_low', 'mid_high', 'high']\n",
    "cols = ['PRICE'] + list(boston_clusters.columns[:-1])\n",
    "boston_clusters = boston_clusters[cols]\n",
    "boston_clusters[['PRICE', 'CRIM', 'RM', 'AGE', 'DIS']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first brush, it looks like price is pretty well differentiated by cluster.  Remember that price did not participate in the clustering model; so inasmuch as they correspond, it confirms our working hypothesis.\n",
    "\n",
    "Other columns have interstingly different patterns.  Crime rate is high for the \"low\" category houses, then drops to almost none for \"mid_low\", \"mid_high\", and \"high.\"  In contrast, the number of rooms varies little by category.  Distance to employment centers is low and similar for the first three categories of houses, but much larger for the \"high\" category houses.  Age of house shows a similar, but less dramatic, pattern as distance.\n",
    "\n",
    "Let us try to get a handle on who well price is genuinely differentiated by category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = boston_df.groupby('category').PRICE.agg(['mean', 'median', 'std']).sort_values('mean')\n",
    "clusters.index = ['low', 'mid_low', 'mid_high', 'high']\n",
    "clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a visual representation of the pattern we see *pretty good* but not overwhelming differentiation. One category of houses show a wide range of prices, while other categories tend to be more narrowly centered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "sns.distplot(boston_df[boston_df.category == 0].PRICE, hist=False, rug=True)\n",
    "sns.distplot(boston_df[boston_df.category == 1].PRICE, hist=False, rug=True)\n",
    "sns.distplot(boston_df[boston_df.category == 2].PRICE, hist=False, rug=True)\n",
    "sns.distplot(boston_df[boston_df.category == 3].PRICE, hist=False, rug=True)\n",
    "ax.set_title(\"Distribution of clusters on PRICE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given our fairly successful confirmation of our hypothesis in regard to 1970s Boston housing prices  relative to independent features, let us try the same exercise with the much larger 1990s California housing prices dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "california = datasets.california_housing.fetch_california_housing()\n",
    "ca_df = pd.DataFrame(california.data, columns=california.feature_names)\n",
    "ca_df['PRICE'] = california.target\n",
    "ca_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Takes a little while to cluster 20k observations\n",
    "from sklearn import cluster\n",
    "spectral = cluster.SpectralClustering(\n",
    "        n_clusters=4, eigen_solver='arpack',\n",
    "        affinity=\"nearest_neighbors\")\n",
    "\n",
    "X = StandardScaler().fit_transform(california.data)\n",
    "spectral.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us do minor Pandas manipulation to view clusters in an intuitive way\n",
    "ca_df['category'] = spectral.labels_\n",
    "ca_clusters = ca_df.groupby('category').mean().sort_values('PRICE')\n",
    "ca_clusters.index = ['low', 'mid_low', 'mid_high', 'high']\n",
    "ca_clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The clustering here seems to be mostly independent of price.  The only things there really seem differentiated in the cluster categories are age of house (young houses at extremes of price scale), number of rooms, and especially size of the cities where the houses are located.  Even the north/south or east/west distinctions I would expect in California houses are not strongly exposed by this technique.\n",
    "\n",
    "We can look more specifically at the price distributions in groups, but it shows nothing relevant (nor does the more detailed graph; other than a somewhat bimodel distribution for all categories)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ca_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std'])\n",
    "     .sort_values('median')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "sns.distplot(ca_df[ca_df.category == 0].PRICE, \n",
    "             label=\"0\", hist=False, rug=True)\n",
    "sns.distplot(ca_df[ca_df.category == 1].PRICE, \n",
    "             label=\"1\", hist=False, rug=True)\n",
    "sns.distplot(ca_df[ca_df.category == 2].PRICE, \n",
    "             label=\"2\", hist=False, rug=True)\n",
    "sns.distplot(ca_df[ca_df.category == 3].PRICE, \n",
    "             label=\"3\", hist=False, rug=True)\n",
    "ax.set_title(\"Distribution of clusters on PRICE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster into an Unknown Number of Categories\n",
    "\n",
    "In the first pass, we chose the number of clusters we expect in the data manually.  Many models require an `n_clusters` argument to select this.  However, some of the models determine the number of clusters from the data itself.  That is not uniformly better, but may let the \"data speak\" better in many cases.\n",
    "\n",
    "Let us try our same analysis using the very good `DBSCAN` algorithm (Density-based spatial clustering of applications with noise).  This model is often the go-to choice for data scientists:\n",
    "\n",
    "> The [DBSCAN](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html) algorithm views clusters as areas of high density separated by areas of low density. Due to this rather generic view, clusters found by DBSCAN can be any shape, as opposed to k-means which assumes that clusters are convex shaped. The central component to the DBSCAN is the concept of core samples, which are samples that are in areas of high density. A cluster is therefore a set of core samples, each close to each other (measured by some distance measure) and a set of non-core samples that are close to a core sample (but are not themselves core samples).\n",
    "\n",
    "A result of this algorithm is that the number of clusters detected depends strongly on an epsilon hyperparameter that measures the acceptable distance to extend the region of a \"core\" and to reach non-core samples.  This makes scaling important to get consistent units among features, but also makes the choice of epsilon have a strong effect on number of identified clusters.  Notwithstanding the simple illustration below, the \"core\" region can be extended in any shape by identifying additional core points along any path.\n",
    "\n",
    "\"Noise points\" in DBSCAN are given the class label -1.\n",
    "\n",
    "<img src=\"img/DBSCAN.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "X = StandardScaler().fit_transform(boston.data)\n",
    "\n",
    "for epsilon in np.arange(0.1, 1.8, 0.1):\n",
    "    db = cluster.DBSCAN(eps=epsilon)\n",
    "    db.fit(X)\n",
    "    labels = np.unique(db.labels_)\n",
    "    print(\"Epsilon=%.1f; labels=[%d .. %d]\" % (epsilon, labels[0], labels[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "X = StandardScaler().fit_transform(california.data)\n",
    "\n",
    "for epsilon in np.arange(0.1, 1.8, 0.1):\n",
    "    db = cluster.DBSCAN(eps=epsilon)\n",
    "    db.fit(X)\n",
    "    labels = np.unique(db.labels_)\n",
    "    print(\"Epsilon=%.1f; labels=[%d .. %d]\" % (epsilon, labels[0], labels[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluating Price Hypothesis with DBSCAN\n",
    "\n",
    "As you can see, by choosing different epsilon values, we can get wildly different numbers of clusters using DBSCAN.  This is not a bad thing, and being able to choose the granularity of clustering can be desirable.  However, on closer examination, we notice that  with various epsilon values, we still wind up either with most data samples in category -1 for noise, or in the overly generic category 0.  It may be that further hyperparameter tuning could improve this, but we will not do that for the California dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = StandardScaler().fit_transform(california.data)\n",
    "ca_df['category'] = cluster.DBSCAN(eps=0.5).fit(X).labels_\n",
    "\n",
    "(ca_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std', 'count'])\n",
    "     .sort_values('count', ascending=False)\n",
    "     .head()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Increasing `min_samples` from the default 5 to 20 *does* improve matters somewhat. Other hyperparameter values might continue to improve this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = StandardScaler().fit_transform(california.data)\n",
    "ca_df['category'] = cluster.DBSCAN(eps=0.5, min_samples=20).fit(X).labels_\n",
    "\n",
    "(ca_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std', 'count'])\n",
    "     .sort_values('count', ascending=False)\n",
    "     .head()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the Boston dataset we can get much less problem with noise and dominance of cluster 0.  Of course, we already had reasonable confirmation of our initial price-based clustering with spectral clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = StandardScaler().fit_transform(boston.data)\n",
    "boston_df['category'] = cluster.DBSCAN(eps=1.6).fit(X).labels_\n",
    "\n",
    "(boston_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std', 'count'])\n",
    "     .sort_values('count', ascending=False)\n",
    "     .head()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exclude the noise =-1 in the graph\n",
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "for cat in range(boston_df.category.max()+1):\n",
    "    sns.distplot(boston_df[boston_df.category == cat].PRICE, \n",
    "                 label=str(cat), hist=False, rug=True)\n",
    "ax.set_title(\"Distribution of clusters on PRICE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I believe our conclusion at this point should be that there are several price-centric categories of houses in the Boston dataset, but for other classes that are clustered, different factors take precedent.  Five or six of the 10 categories have steep Gaussian slopes, that rest are distributed in ways suggesting price is not crucial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HDBSCAN: A Better DBSCAN?\n",
    "\n",
    "A colleague of mine advised me to try out an improved variant of DBSCAN called [HDBSCAN](https://hdbscan.readthedocs.io/en/latest/index.html).  This algorithm is not currently included in scikit-learn, but it follows the same API.  The basic idea in HDBSCAN is that the algorithm will dynamically adjust the epsilon value in regions of varying density.  From the [Campello, Moulavi, and Sander](http://link.springer.com/chapter/10.1007%2F978-3-642-37456-2_14) research paper that introduced it:\n",
    "\n",
    "> We propose a theoretically and practically improved density-based, hierarchical clustering method, providing a clustering hierarchy from which a simplified tree of significant clusters can be constructed. For obtaining a “flat” partition consisting of only the most significant clusters (possibly corresponding to different density thresholds), we propose a novel cluster stability measure, formalize the problem of maximizing the overall stability of selected clusters, and formulate an algorithm that computes an optimal solution to this problem. We demonstrate that our approach outperforms the current, state-of-the-art, density-based clustering methods on a wide variety of real world data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hdbscan import HDBSCAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = StandardScaler().fit_transform(boston.data)\n",
    "\n",
    "hdb = HDBSCAN()\n",
    "hdb.fit(X)\n",
    "labels = np.unique(hdb.labels_)\n",
    "print(\"Labels=[%d .. %d]\" % (labels[0], labels[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_df['category'] = hdb.labels_\n",
    "(boston_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std', 'count'])\n",
    "     .sort_values('count', ascending=False)\n",
    "     .head()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exclude the noise =-1 in the graph\n",
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "for cat in range(boston_df.category.max()+1):\n",
    "    sns.distplot(boston_df[boston_df.category == cat].PRICE, \n",
    "                 label=str(cat), hist=False, rug=True)\n",
    "ax.set_title(\"Distribution of clusters on PRICE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the relatively confirmatory Boston housing dataset, HDBSCAN produces sharper price peaks than other clustering approaches we have looked at. However, as with other techniques, it appears to identify several price centers *and also* other clusters whose identity is not chiefly price-focused.\n",
    "\n",
    "HDBSCAN both gave us automated determination of number of clusters and avoided the need for manually tweaking epsilon. But the Boston dataset has been easy, let us try against the California housing dataset.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = StandardScaler().fit_transform(california.data)\n",
    "\n",
    "hdb = HDBSCAN(min_cluster_size=20)\n",
    "hdb.fit(X)\n",
    "labels = np.unique(hdb.labels_)\n",
    "print(\"Labels=[%d .. %d]\" % (labels[0], labels[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca_df['category'] = hdb.labels_\n",
    "(ca_df\n",
    "     .groupby('category')\n",
    "     .PRICE\n",
    "     .agg(['mean', 'median', 'std', 'count'])\n",
    "     .sort_values('count', ascending=False)\n",
    "     .head()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This approach is still relatively unsucessessful.  At least it does not identify the assumed price-basis of the clusters.  It may be time to reject our hypothesis relative to this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exclude the noise =-1 in the graph\n",
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "for cat in range(ca_df.category.max()+1):\n",
    "    sns.distplot(ca_df[ca_df.category == cat].PRICE, \n",
    "                 label=str(cat), hist=False, rug=True)\n",
    "ax.set_title(\"Distribution of clusters on PRICE\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating Clustering\n",
    "\n",
    "From the [scikit-learn documentation](http://scikit-learn.org/stable/modules/clustering.htm):\n",
    "\n",
    "> Evaluating the performance of a clustering algorithm is not as trivial as counting the number of errors or the precision and recall of a supervised classification algorithm. In particular any evaluation metric should not take the absolute values of the cluster labels into account but rather if this clustering define separations of the data similar to some ground truth set of classes or satisfying some assumption such that members belong to the same class are more similar that members of different classes according to some similarity metric.\n",
    "\n",
    "The problem, of course, is that in **unsupervised** learning, you are more-or-less assured of not having a \"ground truth\" since you are trying to tease that out of the data.  In the examples above, we partially faked it by having a hypothesis about ground truth; but that hypothesis seemed supportable in one dataset and unsupportable in a different related dataset.$^1$\n",
    "\n",
    "<hr/>\n",
    "<small>$^1$<i>It is accurate to say that we have failed to validate our hypothesis, not that we have proven it false.  It might still be that some additional clustering model we have not tried—or some hyperparameters to those we tried—would produce price-linked clusters in the California dataset of similar strength to those the Boston dataset readily showed.</i></small>\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One kind of validation that is possible to perform is to see whether different clustering algorithms produce similar classes, and assign observations to the same class as each other.  If there is agreement, that lends *some* extra weight to the idea that the clusters represent an underlying phenomenon of the data.\n",
    "\n",
    "The Rand Index computes a similarity measure between two clusterings by considering all pairs of samples and counting pairs that are assigned in the same or different clusters in the predicted and true clusterings.  This method is symmetric and independent of the particular labels an algorithm uses.  It depends only on collections of samples in the same class under one regime correspondng to the same collection under another regime.  *Adjusted mutual information score* is similar to *adjusted Rand score*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit several clustering models that take n_clusters\n",
    "from sklearn import cluster\n",
    "\n",
    "X = StandardScaler().fit_transform(boston.data)\n",
    "\n",
    "models = [cluster.SpectralClustering(n_clusters=4), \n",
    "          cluster.KMeans(n_clusters=4), \n",
    "          cluster.MiniBatchKMeans(n_clusters=4), \n",
    "          cluster.AgglomerativeClustering(n_clusters=4)]\n",
    "\n",
    "for model in models:\n",
    "    model.fit(X)\n",
    "    print(model.__class__.__name__)\n",
    "    print(\"\\tFirst 5 labels\", model.labels_[:5])\n",
    "    print(\"\\t\", len(model.labels_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example above, all four clustering models selected assigned the same class to the first five samples.  However, the particular number given to that class might be any of the four possible values (depending even on random seeds).  The label numbers are not interesting, only the samples assigned to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ... then measure their pairwise similarity\n",
    "from sklearn import metrics\n",
    "from itertools import combinations\n",
    "\n",
    "for clust1, clust2 in combinations(models, 2):\n",
    "    print(clust1.__class__.__name__, \"versus\", clust2.__class__.__name__)\n",
    "    print(\"\\tRand score:\", metrics.adjusted_rand_score(clust1.labels_, clust2.labels_))\n",
    "    print(\"\\tMutual info:\", metrics.adjusted_mutual_info_score(clust1.labels_, clust2.labels_))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both measures made range between zero and one.  As one would expect, the closely related `KMeans` and `MiniBatchKMeans` are closely similar.  Other pairs are moderately similar; certainly they are much more similar than e.g. random choices of labels, which would produce scores of zero (the *adjusted* part discounts for expected random intersections)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next Lesson\n",
    "\n",
    "**Feature Engineering**: In the current lessson we looked at several clustering models.  Besides their underlying algorithmic differences, a key difference among them is whether they choose a number of clusters based on the data itself or produce a specified number of classes. \n",
    "\n",
    "In this lesson and some previous one, we have done some feature engineering in passing. In the next lesson we will address dimensionality reduction, feature engineering, and feature selection, in more specific detail.\n",
    "\n",
    "<a href=\"Features.ipynb\"><img src=\"img/open-notebook.png\" align=\"left\"/></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}