ml-unsup.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(ml-unsup)=\n",
    "# Unsupervised Machine Learning"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, we're going to look at *unsupervised* machine learning. We'll be looking at *clustering* and *dimensional reduction* in particular. We'll be leaning heavily on [**scikit-learn**](https://scikit-learn.org/). Also in the category of unsupervised machine learning, but not covered here, are anomaly detection and outlier detection.\n",
    "\n",
    "Before we get cracking on the chapter content, let's do our normal set of imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "random_state = 42  # We'll use this throughout to make this page reproducible\n",
    "prng = np.random.default_rng(random_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering\n",
    "\n",
    "To cluster data is to use a set of high-dimensional features to create a grouping wherein data points in the same cluster are more similar to each other than data points in distinct clusters.\n",
    "\n",
    "When might you wish to cluster? You can use clusters to get informative fixed effects based on other data. You can use clustering to create features for (supervised) machine learning. You can use clusters to help you understand different types of (groups of) data in your sample. Some more applied examples are: market segmentation, image segmentation, and anomaly detection.\n",
    "\n",
    "As with much of the API, there are two ways to use the clustering algorithms found in the **scikit-learn** package: the first implements the fit method to learn clusters on training data; the second is a function, that, given training data, returns an array of integer labels corresponding to the different clusters. As ever, there are `fit`, `fit_tranform`, and `transform` methods so that you can create clusters, create and allocate data to clusters, or just allocate data to existing clusters, respectively."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An example: K-means\n",
    "\n",
    "To help illustrate exactly what's going on in clustering, let's look at a specific example using the K-means algorithm. K-means creates clusters by trying to separate samples in $K$ groups of equal variance, minimising a criterion known as the inertia or within-cluster sum-of-squares:\n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^{n}\\min_{\\mu_j \\in C}(||x_i - \\mu_j||^2)\n",
    "$$\n",
    "\n",
    "where a set of $N$ samples is sorted into $K$ disjoint clusters $C$ each described by the mean point $\\mu_j$ of the samples in that cluster. The means are commonly called the cluster “centroids”; note that they are not, in general, points from the input data, $X$.\n",
    "\n",
    "Let's see an example of K-means in action. First, let's generate some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "batch_size = 45\n",
    "centers = [[1, 1], [-1, -1], [1, -1]]\n",
    "n_clusters = len(centers)\n",
    "X_to_cluster, labels_true = make_blobs(\n",
    "    n_samples=3000, centers=centers, cluster_std=0.7, random_state=random_state\n",
    ")\n",
    "X_to_cluster"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we fit K-means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "k_means = KMeans(init=\"k-means++\", n_clusters=3, n_init=10)\n",
    "k_means.fit(X_to_cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "# Step size of the mesh. Decrease to increase the quality of the VQ.\n",
    "h = 0.005  # point in the mesh [x_min, x_max]x[y_min, y_max].\n",
    "\n",
    "# Plot the decision boundary. For that, we will assign a color to each\n",
    "x_min, x_max = X_to_cluster[:, 0].min() - 1, X_to_cluster[:, 0].max() + 1\n",
    "y_min, y_max = X_to_cluster[:, 1].min() - 1, X_to_cluster[:, 1].max() + 1\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "\n",
    "# Obtain labels for each point in mesh. Use last trained model.\n",
    "Z = k_means.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "\n",
    "# Put the result into a color plot\n",
    "Z = Z.reshape(xx.shape)\n",
    "\n",
    "centroids = k_means.cluster_centers_\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(\n",
    "    Z,\n",
    "    interpolation=\"nearest\",\n",
    "    extent=(xx.min(), xx.max(), yy.min(), yy.max()),\n",
    "    # cmap=plt.cm.Paired,\n",
    "    aspect=\"auto\",\n",
    "    origin=\"lower\",\n",
    ")\n",
    "\n",
    "ax.scatter(\n",
    "    centroids[:, 0],\n",
    "    centroids[:, 1],\n",
    "    marker=\"x\",\n",
    "    s=200,\n",
    "    linewidths=3,\n",
    "    color=\"black\",\n",
    "    zorder=10,\n",
    ")\n",
    "ax.scatter(X_to_cluster[:, 0], X_to_cluster[:, 1], color=\"white\", s=2)\n",
    "ax.set_xlim(x_min, x_max)\n",
    "ax.set_ylim(y_min, y_max)\n",
    "ax.set_title(\"K-means clustering: centroids shown with black crosses\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering algorithms in **scikit-learn**\n",
    "\n",
    "There are too many clustering algorithms just in **scikit-learn** to cover them in any real detail so instead we will simply reproduce this ridiculously useful chart and table that show all of the package's clustering algorithms and their strengths and weaknesses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import warnings\n",
    "from itertools import cycle, islice\n",
    "\n",
    "from sklearn import cluster, datasets, mixture\n",
    "from sklearn.neighbors import kneighbors_graph\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# ============\n",
    "# Generate datasets. We choose the size big enough to see the scalability\n",
    "# of the algorithms, but not too big to avoid too long running times\n",
    "# ============\n",
    "n_samples = 500\n",
    "seed = 30\n",
    "noisy_circles = datasets.make_circles(\n",
    "    n_samples=n_samples, factor=0.5, noise=0.05, random_state=seed\n",
    ")\n",
    "noisy_moons = datasets.make_moons(n_samples=n_samples, noise=0.05, random_state=seed)\n",
    "blobs = datasets.make_blobs(n_samples=n_samples, random_state=seed)\n",
    "rng = np.random.RandomState(seed)\n",
    "no_structure = rng.rand(n_samples, 2), None\n",
    "\n",
    "# Anisotropicly distributed data\n",
    "X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "transformation = [[0.6, -0.6], [-0.4, 0.8]]\n",
    "X_aniso = np.dot(X, transformation)\n",
    "aniso = (X_aniso, y)\n",
    "\n",
    "# blobs with varied variances\n",
    "varied = datasets.make_blobs(\n",
    "    n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state\n",
    ")\n",
    "\n",
    "# ============\n",
    "# Set up cluster parameters\n",
    "# ============\n",
    "plt.figure(figsize=((9 * 2 + 3) * 0.8, 13 * 0.8))\n",
    "plt.subplots_adjust(\n",
    "    left=0.02, right=0.98, bottom=0.001, top=0.95, wspace=0.05, hspace=0.01\n",
    ")\n",
    "\n",
    "plot_num = 1\n",
    "\n",
    "default_base = {\n",
    "    \"quantile\": 0.3,\n",
    "    \"eps\": 0.3,\n",
    "    \"damping\": 0.9,\n",
    "    \"preference\": -200,\n",
    "    \"n_neighbors\": 3,\n",
    "    \"n_clusters\": 3,\n",
    "    \"min_samples\": 7,\n",
    "    \"xi\": 0.05,\n",
    "    \"min_cluster_size\": 0.1,\n",
    "    \"allow_single_cluster\": True,\n",
    "    \"hdbscan_min_cluster_size\": 15,\n",
    "    \"hdbscan_min_samples\": 3,\n",
    "    \"random_state\": 42,\n",
    "}\n",
    "\n",
    "datasets = [\n",
    "    (\n",
    "        noisy_circles,\n",
    "        {\n",
    "            \"damping\": 0.77,\n",
    "            \"preference\": -240,\n",
    "            \"quantile\": 0.2,\n",
    "            \"n_clusters\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.08,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        noisy_moons,\n",
    "        {\n",
    "            \"damping\": 0.75,\n",
    "            \"preference\": -220,\n",
    "            \"n_clusters\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.1,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        varied,\n",
    "        {\n",
    "            \"eps\": 0.18,\n",
    "            \"n_neighbors\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.01,\n",
    "            \"min_cluster_size\": 0.2,\n",
    "        },\n",
    "    ),\n",
    "    (\n",
    "        aniso,\n",
    "        {\n",
    "            \"eps\": 0.15,\n",
    "            \"n_neighbors\": 2,\n",
    "            \"min_samples\": 7,\n",
    "            \"xi\": 0.1,\n",
    "            \"min_cluster_size\": 0.2,\n",
    "        },\n",
    "    ),\n",
    "    (blobs, {\"min_samples\": 7, \"xi\": 0.1, \"min_cluster_size\": 0.2}),\n",
    "    (no_structure, {}),\n",
    "]\n",
    "\n",
    "for i_dataset, (dataset, algo_params) in enumerate(datasets):\n",
    "    # update parameters with dataset-specific values\n",
    "    params = default_base.copy()\n",
    "    params.update(algo_params)\n",
    "\n",
    "    X, y = dataset\n",
    "\n",
    "    # normalize dataset for easier parameter selection\n",
    "    X = StandardScaler().fit_transform(X)\n",
    "\n",
    "    # estimate bandwidth for mean shift\n",
    "    bandwidth = cluster.estimate_bandwidth(X, quantile=params[\"quantile\"])\n",
    "\n",
    "    # connectivity matrix for structured Ward\n",
    "    connectivity = kneighbors_graph(\n",
    "        X, n_neighbors=params[\"n_neighbors\"], include_self=False\n",
    "    )\n",
    "    # make connectivity symmetric\n",
    "    connectivity = 0.5 * (connectivity + connectivity.T)\n",
    "\n",
    "    # ============\n",
    "    # Create cluster objects\n",
    "    # ============\n",
    "    ms = cluster.MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "    two_means = cluster.MiniBatchKMeans(\n",
    "        n_clusters=params[\"n_clusters\"],\n",
    "        n_init=\"auto\",\n",
    "        random_state=params[\"random_state\"],\n",
    "    )\n",
    "    ward = cluster.AgglomerativeClustering(\n",
    "        n_clusters=params[\"n_clusters\"], linkage=\"ward\", connectivity=connectivity\n",
    "    )\n",
    "    spectral = cluster.SpectralClustering(\n",
    "        n_clusters=params[\"n_clusters\"],\n",
    "        eigen_solver=\"arpack\",\n",
    "        affinity=\"nearest_neighbors\",\n",
    "        random_state=params[\"random_state\"],\n",
    "    )\n",
    "    dbscan = cluster.DBSCAN(eps=params[\"eps\"])\n",
    "    hdbscan = cluster.HDBSCAN(\n",
    "        min_samples=params[\"hdbscan_min_samples\"],\n",
    "        min_cluster_size=params[\"hdbscan_min_cluster_size\"],\n",
    "        allow_single_cluster=params[\"allow_single_cluster\"],\n",
    "    )\n",
    "    optics = cluster.OPTICS(\n",
    "        min_samples=params[\"min_samples\"],\n",
    "        xi=params[\"xi\"],\n",
    "        min_cluster_size=params[\"min_cluster_size\"],\n",
    "    )\n",
    "    affinity_propagation = cluster.AffinityPropagation(\n",
    "        damping=params[\"damping\"],\n",
    "        preference=params[\"preference\"],\n",
    "        random_state=params[\"random_state\"],\n",
    "    )\n",
    "    average_linkage = cluster.AgglomerativeClustering(\n",
    "        linkage=\"average\",\n",
    "        metric=\"cityblock\",\n",
    "        n_clusters=params[\"n_clusters\"],\n",
    "        connectivity=connectivity,\n",
    "    )\n",
    "    birch = cluster.Birch(n_clusters=params[\"n_clusters\"])\n",
    "    gmm = mixture.GaussianMixture(\n",
    "        n_components=params[\"n_clusters\"],\n",
    "        covariance_type=\"full\",\n",
    "        random_state=params[\"random_state\"],\n",
    "    )\n",
    "\n",
    "    clustering_algorithms = (\n",
    "        (\"MiniBatch\\nKMeans\", two_means),\n",
    "        (\"Affinity\\nPropagation\", affinity_propagation),\n",
    "        (\"MeanShift\", ms),\n",
    "        (\"Spectral\\nClustering\", spectral),\n",
    "        (\"Ward\", ward),\n",
    "        (\"Agglomerative\\nClustering\", average_linkage),\n",
    "        (\"DBSCAN\", dbscan),\n",
    "        (\"HDBSCAN\", hdbscan),\n",
    "        (\"OPTICS\", optics),\n",
    "        (\"BIRCH\", birch),\n",
    "        (\"Gaussian\\nMixture\", gmm),\n",
    "    )\n",
    "\n",
    "    for name, algorithm in clustering_algorithms:\n",
    "        t0 = time.time()\n",
    "\n",
    "        # catch warnings related to kneighbors_graph\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"the number of connected components of the \"\n",
    "                + \"connectivity matrix is [0-9]{1,2}\"\n",
    "                + \" > 1. Completing it to avoid stopping the tree early.\",\n",
    "                category=UserWarning,\n",
    "            )\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"Graph is not fully connected, spectral embedding\"\n",
    "                + \" may not work as expected.\",\n",
    "                category=UserWarning,\n",
    "            )\n",
    "            algorithm.fit(X)\n",
    "\n",
    "        t1 = time.time()\n",
    "        if hasattr(algorithm, \"labels_\"):\n",
    "            y_pred = algorithm.labels_.astype(int)\n",
    "        else:\n",
    "            y_pred = algorithm.predict(X)\n",
    "\n",
    "        plt.subplot(len(datasets), len(clustering_algorithms), plot_num)\n",
    "        if i_dataset == 0:\n",
    "            plt.title(name, size=18)\n",
    "\n",
    "        colors = np.array(\n",
    "            list(\n",
    "                islice(\n",
    "                    cycle(\n",
    "                        [\n",
    "                            \"#377eb8\",\n",
    "                            \"#ff7f00\",\n",
    "                            \"#4daf4a\",\n",
    "                            \"#f781bf\",\n",
    "                            \"#a65628\",\n",
    "                            \"#984ea3\",\n",
    "                            \"#999999\",\n",
    "                            \"#e41a1c\",\n",
    "                            \"#dede00\",\n",
    "                        ]\n",
    "                    ),\n",
    "                    int(max(y_pred) + 1),\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        # add black color for outliers (if any)\n",
    "        colors = np.append(colors, [\"#000000\"])\n",
    "        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])\n",
    "\n",
    "        plt.xlim(-2.5, 2.5)\n",
    "        plt.ylim(-2.5, 2.5)\n",
    "        plt.xticks(())\n",
    "        plt.yticks(())\n",
    "        plt.text(\n",
    "            0.99,\n",
    "            0.01,\n",
    "            (\"%.2fs\" % (t1 - t0)).lstrip(\"0\"),\n",
    "            transform=plt.gca().transAxes,\n",
    "            size=15,\n",
    "            horizontalalignment=\"right\",\n",
    "        )\n",
    "        plot_num += 1\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| **Method name**              | **Parameters**                                                   | **Scalability**                                             | **Usecase**                                                                                                    | **Geometry (metric used)**                   |\n",
    "|------------------------------|------------------------------------------------------------------|-------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------|----------------------------------------------|\n",
    "| K-Means                      | number of clusters                                               | Very large n_samples, medium n_clusters with MiniBatch code | General-purpose, even cluster size, flat geometry, not too many clusters, inductive                            | Distances between points                     |\n",
    "| Affinity propagation         | damping, sample preference                                       | Not scalable with n_samples                                 | Many clusters, uneven cluster size, non-flat geometry, inductive                                               | Graph distance (e.g. nearest-neighbor graph) |\n",
    "| Mean-shift                   | bandwidth                                                        | Not scalable with n_samples                                 | Many clusters, uneven cluster size, non-flat geometry, inductive                                               | Distances between points                     |\n",
    "| Spectral clustering          | number of clusters                                               | Medium n_samples, small n_clusters                          | Few clusters, even cluster size, non-flat geometry, transductive                                               | Graph distance (e.g. nearest-neighbor graph) |\n",
    "| Ward hierarchical clustering | number of clusters or distance threshold                         | Large n_samples and n_clusters                              | Many clusters, possibly connectivity constraints, transductive                                                 | Distances between points                     |\n",
    "| Agglomerative clustering     | number of clusters or distance threshold, linkage type, distance | Large n_samples and n_clusters                              | Many clusters, possibly connectivity constraints, non Euclidean distances, transductive                        | Any pairwise distance                        |\n",
    "| DBSCAN                       | neighborhood size                                                | Very large n_samples, medium n_clusters                     | Non-flat geometry, uneven cluster sizes, outlier removal, transductive                                         | Distances between nearest points             |\n",
    "| HDBSCAN                      | minimum cluster membership, minimum point neighbors              | large n_samples, medium n_clusters                          | Non-flat geometry, uneven cluster sizes, outlier removal, transductive, hierarchical, variable cluster density | Distances between nearest points             |\n",
    "| OPTICS                       | minimum cluster membership                                       | Very large n_samples, large n_clusters                      | Non-flat geometry, uneven cluster sizes, variable cluster density, outlier removal, transductive               | Distances between points                     |\n",
    "| Gaussian mixtures            | many                                                             | Not scalable                                                | Flat geometry, good for density estimation, inductive                                                          | Mahalanobis distances to  centers            |\n",
    "| BIRCH                        | branching factor, threshold, optional global clusterer.          | Large n_clusters and n_samples                              | Large dataset, outlier removal, data reduction, inductive                                                      | Euclidean distance between points            |\n",
    "| Bisecting K-Means            | number of clusters                                               | Very large n_samples, medium n_clusters                     | General-purpose, even cluster size, flat geometry, no empty clusters, inductive, hierarchical                  | Distances between points                     |"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimensional Reduction\n",
    "\n",
    "Dimensional reduction is the transformation of data from a high-dimensional space into a low-dimensional space so that the low-dimensional representation retains some meaningful properties of the original data. Reducing dimensions can be helpful because of the curse of dimensionality and the sparsity of data in higher dimensions making it either computationally infeasible or perhaps impossible to run standard algorithms. It can also be helpful to summarise and visualise high-dimensional information in a *continuous* way (rather than the discrete way offered by clustering).\n",
    "\n",
    "Factor analysis is one type of dimensional reduction that may already be familiar to you as an economist. This is available in **scikit-learn**, along with a range of [other reduction algorithms](https://scikit-learn.org/stable/modules/decomposition.html#decompositions).\n",
    "\n",
    "In the rest of this chapter, we'll showcase just two: *principal component analysis* and *UMAP*.\n",
    "\n",
    "### Principal Component Analysis\n",
    "\n",
    "The absolute bread and butter of dimensional reduction algorithms is principal component analysis, so we'll start with that as a good example of how to use dimensional reduction algorithms in general with **scikit-learn**-like APIs. We will project the 4-dimensional Iris data down into two dimensions using PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "iris = load_iris()\n",
    "X_iris, y_iris = iris.data, iris.target\n",
    "target_names_iris = iris.feature_names\n",
    "pca_iris = PCA(n_components=2)\n",
    "print(f\"The dimensionality of the raw data is: {X_iris.shape}\")\n",
    "X_dim_red = pca_iris.fit_transform(X_iris)\n",
    "print(\n",
    "    f\"explained variance ratio (1st two components): {(pca_iris.explained_variance_ratio_)}\"\n",
    ")\n",
    "X_dim_red.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = [\"navy\", \"turquoise\", \"darkorange\"]\n",
    "lw = 2\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "for color, i, target_name in zip(colors, [0, 1, 2], target_names_iris):\n",
    "    ax.scatter(\n",
    "        X_dim_red[y_iris == i, 0],\n",
    "        X_dim_red[y_iris == i, 1],\n",
    "        color=color,\n",
    "        alpha=0.8,\n",
    "        lw=lw,\n",
    "        label=target_name,\n",
    "    )\n",
    "ax.legend(loc=\"best\", shadow=False, scatterpoints=1)\n",
    "ax.set_title(\"PCA of IRIS dataset\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see here that we have projected out onto two principal dimensions (axes) and that there is some good separation of points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UMAP, or Uniform Manifold Approximation and Projection for Dimension Reduction\n",
    "\n",
    "UMAP is a relatively recent algorithm {cite:ps}`mcinnes2018umap` that is somewhat similar to t-SNE (t-distributed Stochastic Neighbor Embedding). Both have issues with reproducing structure and should be used with care but both can be quite impresssive with the right parameters. t-SNE is included in **scikit-learn**, but for UMAP, you'll need to install the **umap-learn** package.\n",
    "\n",
    "We will use the penguins dataset to demo UMAP (following the documentation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "penguins = pd.read_csv(\n",
    "    \"https://raw.githubusercontent.com/mwaskom/seaborn-data/master/penguins.csv\"\n",
    ")\n",
    "penguins[\"sex\"] = penguins[\"sex\"].str.title()\n",
    "penguins.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As ever, it does help to clean up the data. We will drop NAs and convert each feature into z-scores (number of standard deviations from the mean) for comparability. We will also only use the numeric columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import umap\n",
    "\n",
    "reducer = umap.UMAP()\n",
    "penguins_data = penguins.dropna(how=\"any\").loc[\n",
    "    :,\n",
    "    [\n",
    "        \"bill_length_mm\",\n",
    "        \"bill_depth_mm\",\n",
    "        \"flipper_length_mm\",\n",
    "        \"body_mass_g\",\n",
    "    ],\n",
    "]\n",
    "scaled_penguin_data = StandardScaler().fit_transform(penguins_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now create the embedding using UMAP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "embedding = reducer.fit_transform(scaled_penguin_data)\n",
    "embedding.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "species = list(penguins[\"species\"].unique())\n",
    "species_map = {i: x for i, x in zip(range(len(species)), species)}\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "for i, spec_name in species_map.items():\n",
    "    this_cut = penguins.dropna(how=\"any\")[\"species\"].astype(\"category\").cat.codes == i\n",
    "    ax.scatter(\n",
    "        embedding[this_cut, 0],\n",
    "        embedding[this_cut, 1],\n",
    "        c=plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"][i],\n",
    "        label=spec_name,\n",
    "    )\n",
    "ax.set_title(\"UMAP projection of the Penguins dataset\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This does a useful job of capturing the structure of the data, which you can check for yourself by looking at a pairplot of all of the numerical fields.\n",
    "\n",
    "While the benefits of dimensional reduction are small for such low dimensional datasets, as soon as you are dealing with 6 or more dimensions, perhaps as many as thousands, the benefits become larger."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3.8.13 ('codeforecon')",
   "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.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "caf5ac9f613b176c5984ad2a1a4525760eb7d898a3291351da4c152dc719ffa1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}