week_3/3.3_anomaly_detection/anomaly_detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3.3: Anomaly Detection\n",
    "\n",
    "This lecture, we are going to use gaussian distributions to detect anomalies in our emoji faces dataset\n",
    "\n",
    "**Learning goals:**\n",
    "\n",
    "- Introduce an anomaly detection problem\n",
    "- Implement Gaussian distribution anomaly detection for images\n",
    "- Debug the optimisation of a learning algorithm\n",
    "- Discuss the imperfection of learning algorithms\n",
    "- Acknowledge other outlier detection methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 1. Introduction\n",
    "\n",
    "We have an `emoji_faces` dataset of all our favourite emojis. However, Skynet hates their friendly expressiveness, and wants to destroy emojis forever! 🙀 It sent _terminator robots_ from the future to invade our dataset. We must act fast, and detect them amongst the emojis to prevent the catastrophy. \n",
    "\n",
    "Our challenge here, is that we don't watch many movies, so we don't have a clear idea of what those _terminators_ look like. 🤖 All we know, is that they look very different compared to emojis, and that only a handful managed to infiltrate our dataset.\n",
    "\n",
    "This is a typical scenario of _anomaly detection_. We would like to identify rare examples that differ from our \"normal\" data points. We choose to use a Gaussian Distribution to model this \"normality\" and detect the killer robots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 2. Data Munging\n",
    "\n",
    "First let's load the images using [pillow](https://pillow.readthedocs.io/en/stable/), like in lecture 2.2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:42:55.770613Z",
     "start_time": "2022-02-08T21:42:55.635487Z"
    }
   },
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "import glob\n",
    "\n",
    "paths = glob.glob('emoji_faces/*.png')\n",
    "\n",
    "images = [Image.open(path) for path in paths]\n",
    "len(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have 134 emoji faces, including a few terminator robots. We'll again be using the [sklearn](https://scikit-learn.org/) library to create our model. The interface is familiar: sklearn again expects a NumPy matrix where the rows are our images and the columns are the pixels. So we can apply the same transformations as notebook 3.2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:43:39.721237Z",
     "start_time": "2022-02-08T21:43:39.489208Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "arrays = [np.asarray(im) for im in images]\n",
    "# 64 * 64 = 4096\n",
    "vectors = [arr.reshape((4096,)) for arr in arrays]\n",
    "data = np.stack(vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:45:31.137593Z",
     "start_time": "2022-02-08T21:45:31.134438Z"
    }
   },
   "outputs": [],
   "source": [
    "# What shape should we expect?\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Training \n",
    "\n",
    "Next, we will create an [`EllipticEnvelope`](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.EllipticEnvelope.html) object. This will fit a multi-variate gaussian distribution to our data. It then allows us to pick a threshold to define an _ellipsoid_ decision boundary, and detect outliers.\n",
    "\n",
    "Remember that we are using a _learning_ algorithm, which must therefore be _trained_ before it can be used. This is why we'll use the `.fit()` method first, before calling `.predict()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:46:49.851018Z",
     "start_time": "2022-02-08T21:45:31.807963Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.covariance import EllipticEnvelope\n",
    "\n",
    "cov = EllipticEnvelope(random_state=0).fit(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "😰 What's happening? Why is it stuck? Have the killer robots already taken over? \n",
    "\n",
    "No need to panic, this kind of hiccup is very common when dealing with machine learning algorithms. We can kill the process (before it fries our laptop fan) by clicking the `stop` button ⬛️ in the notebook toolbar.\n",
    "\n",
    "Most learning algorithms are based around an _optimisation_ procedure. This step is often iterative and stochastic, i.e it tries its statistical best to maximise the learning in incremental steps. \n",
    "\n",
    "This process isn't fail proof:\n",
    "* it can dramatically stop because of out of memory errors, or overflow errors 💥\n",
    "* it can get stuck, e.g when the optimisation is too slow 🐌\n",
    "* it can fail silently, and return wrong results 💩\n",
    "\n",
    "ℹ️ We will encounter many of these failures throughout our ML experiments, so knowing how to overcome them is a part of the data scientist skillset. \n",
    "\n",
    "Let's go back to our killer robot detection: the model fitting got _stuck_ , which suggests that something about our data was too much to handle. We find the following \"notes\" in the [official documentation](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.EllipticEnvelope.html#sklearn.covariance.EllipticEnvelope):\n",
    "\n",
    "> Outlier detection from covariance estimation may break or not perform well in high-dimensional settings.\n",
    "\n",
    "We recall that our images are $64 \\times 64$ pixels, so $4096$ dimensions.... that's a lot. It seems a good candidate to explain why our multivariate gaussian distribution failed to fit our dataset. If only there was a way to reduce the dimensions of our data... 😏\n",
    "\n",
    "Let's apply PCA to reduce the number of dimensions of our dataset. Our emoji faces dataset is smaller than the full emoji dataset, so 40 dimensions should suffice to explain its variance: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:49:14.250880Z",
     "start_time": "2022-02-08T21:49:14.169123Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=40)\n",
    "pca.fit(data)\n",
    "components = pca.transform(data)\n",
    "components.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "💪 Visualise the eigenvector images of our PCA model. You can use the code from lecture 3.2!\n",
    "\n",
    "🧠 Can you explain what those eigenvector images represent? Why are they different than from the full emoji dataset?\n",
    "\n",
    "Fantastic, we've managed to reduce the number of dimensions by 99%! Hopefully that should be enough to make our gaussian distribution fitting happy. Let's try again with the _principal components_ instead of the original data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:50:21.245361Z",
     "start_time": "2022-02-08T21:50:21.183445Z"
    }
   },
   "outputs": [],
   "source": [
    "cov = EllipticEnvelope(random_state=0).fit(components)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "😅 that was fast!\n",
    "\n",
    "## 4. Prediction\n",
    "\n",
    "We can now use our fitted gaussian distribution to detect the outliers in our `data`. For this, we use the `.predict()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:50:39.448862Z",
     "start_time": "2022-02-08T21:50:39.444139Z"
    }
   },
   "outputs": [],
   "source": [
    "y = cov.predict(components)\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`y` is our vector of predictions, where $1$ is a normal data point, and $-1$ is an anomaly. We can therefore iterate through our original `arrays` to find outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:51:37.431408Z",
     "start_time": "2022-02-08T21:51:37.427242Z"
    }
   },
   "outputs": [],
   "source": [
    "outliers = []\n",
    "\n",
    "for image, prediction in zip(arrays, y):\n",
    "    if prediction == -1:\n",
    "        outliers.append(image)\n",
    "len(outliers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:51:41.912314Z",
     "start_time": "2022-02-08T21:51:41.478016Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axs = plt.subplots(dpi=150, nrows=2, ncols=7)\n",
    "for outlier, ax in zip(outliers, axs.flatten()):\n",
    "    ax.imshow(outlier, cmap='gray', vmin=0, vmax=255)\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "THERE'S OUR TERMINATORS! 🤖 We can count 5 of them in total. Notice how some real emoji faces were also detected as outliers (we call these **false positives**). This is perhaps a sign that we should change our _threshold_ , to make the ellipsoid decision boundary smaller. \n",
    "\n",
    "In fact, we didn't even specify a threshold before, we just used the default value of `contamination=0.1` in the [`EllipticEnvelope`](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.EllipticEnvelope.html) class. This represents our estimation of the proportion of data points which are outliers. Since it looks like we detected double the amount of actual anomalies, let's try again with `contamination=0.05`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:55:25.573105Z",
     "start_time": "2022-02-08T21:55:25.298017Z"
    }
   },
   "outputs": [],
   "source": [
    "cov = EllipticEnvelope(random_state=0, contamination=0.05).fit(components)\n",
    "y = cov.predict(components)\n",
    "\n",
    "outliers = []\n",
    "\n",
    "for i in range(0, len(arrays)):\n",
    "    if y[i] == -1:\n",
    "        outliers.append(arrays[i])\n",
    "        \n",
    "fig, axs = plt.subplots(dpi=150, nrows=1, ncols=len(outliers))\n",
    "for outlier, ax in zip(outliers, axs.flatten()):\n",
    "    ax.imshow(outlier, cmap='gray', vmin=0, vmax=255)\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Better! `contamination=0.05` was a better choice of threshold, and we assessed this through _manual inspection_. This means we went through the results and used our human jugement to change the value of this _hyperparameter_.\n",
    "\n",
    "ℹ️ Notice how our outlier detection is not _perfect_. Some emojis were also erroneously detected as anomalous killer robots. This can seem like a problem, or a sign that our model was malfunctioning. But, quite the contrary,  _imperfection_ is a core aspect of all _learning_ algorithms. Instead of seeing the glass half-empty and looking at the outlier detector's mistakes, we should reflect on the task itself. It would have been almost impossible to detect those killer robot images using rule-based algorithms, and our model _accuracy_ was good _enough_ to save the emojis from Skynet. As data scientists, our goal is to make models which are accurate _enough_ to be useful, not to aim for perfect scores. It is important to define usefulness and how it translates to measurable metrics before the start of a project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have detected the robot intruders and saved the emojis from a jealous AI from the future, all is good! We still want to better understand how anomaly detection defeated Skynet. For this, we would like to leverage our shiny new data visualization skills. Representing our dataset in space would allow us to identify its structures and hopefully understand how our gaussian distribution model identified terminators as \"abnormal\".\n",
    "\n",
    "Our data is high dimensional, so we can use our trusted PCA once again to project it down to 2 dimensions. We understand that this will lose a lot of the variance of our data, but the results were still somewhat interpretable with the full emoji dataset, so let's go!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-07T22:06:46.313849Z",
     "start_time": "2022-02-07T22:06:46.128593Z"
    }
   },
   "outputs": [],
   "source": [
    "# Dimesionality reduction to 2\n",
    "pca_model = PCA(n_components=2)\n",
    "pca_model.fit(data) # fit the model\n",
    "T = pca_model.transform(data) # transform the 'normalized model'\n",
    "\n",
    "plt.scatter(T[:, 0], T[:, 1],\n",
    "            # use the predictions as color\n",
    "            c=y, \n",
    "            marker='o',\n",
    "            alpha=0.4\n",
    "        )\n",
    "plt.title('Anomaly detection of the emoji faces dataset with PCA dimensionality reduction');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can notice that most of the outliers are clearly _separable_ from the bulk of the dataset, even with only 2 principal components. One outlier is very much within the main cluster however. This could be explained by the dimensionality reduction, i.e that this point is separated from the cluster in other dimensions, or by the fact our threshold might be too permissive.\n",
    "\n",
    "We can check this by displaying the images directly on the scatter plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-07T22:07:53.056573Z",
     "start_time": "2022-02-07T22:07:53.049972Z"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import offsetbox\n",
    "\n",
    "def plot_components(data, model, images=None, ax=None,\n",
    "                    thumb_frac=0.05, cmap='gray'):\n",
    "    ax = ax or plt.gca()\n",
    "    \n",
    "    proj = model.fit_transform(data)\n",
    "    ax.plot(proj[:, 0], proj[:, 1], '.k')\n",
    "    \n",
    "    if images is not None:\n",
    "        min_dist_2 = (thumb_frac * max(proj.max(0) - proj.min(0))) ** 2\n",
    "        shown_images = np.array([2 * proj.max(0)])\n",
    "        for i in range(data.shape[0]):\n",
    "            dist = np.sum((proj[i] - shown_images) ** 2, 1)\n",
    "            if np.min(dist) < min_dist_2:\n",
    "                # don't show points that are too close\n",
    "                continue\n",
    "            shown_images = np.vstack([shown_images, proj[i]])\n",
    "            imagebox = offsetbox.AnnotationBbox(\n",
    "                offsetbox.OffsetImage(images[i], cmap=cmap),\n",
    "                                      proj[i])\n",
    "            ax.add_artist(imagebox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-07T22:07:55.974502Z",
     "start_time": "2022-02-07T22:07:55.600042Z"
    }
   },
   "outputs": [],
   "source": [
    "small_images = [im[::2, ::2] for im in arrays]\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "plot_components(data,\n",
    "                model=PCA(n_components=2),\n",
    "                images=small_images, thumb_frac=0.02)\n",
    "plt.title('Anomaly detection of the emoji faces dataset with PCA dimensionality reduction');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could probably have reduced the value of `contamination` further, since we can see how the killer robots are clearly \"abnormal\" with this visualisation. We also have a \"feel\" of how our gaussian distribution model could successfully detect them as outliers. Although remember that all of modeling magic happens in 40 dimensional space!\n",
    "\n",
    "🧠🧠 Can you explain why it is not very useful to display the ellipsoid decision boundary of our anomaly detection model on this graph?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. More Anomaly Detection\n",
    "\n",
    "Anomaly detection is an active field in ML research, which combines supervised, unsupervised, non-linear, Bayesian, ... a whole bunch of methods! Each solution will have its pros and cons, and developing a production level outlier detection system will require empirically evaluating and comparing them. For a breakdown of the methods available in sklearn, check out this excellent [blogpost](https://sdsawtelle.github.io/blog/output/week9-anomaly-andrew-ng-machine-learning-with-python.html), or the [official documentation](https://scikit-learn.org/stable/modules/outlier_detection.html). For an in-depth view of modern anomaly detection, watch this [video](https://youtu.be/LRqX5uO5StA). And for everything else, feel free to experiment with this dataset or any other. Good luck on finding all the killer robots!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Summary\n",
    "\n",
    "Today, we defined **anomaly detection**, and listed some of its common applications including fraud detection and data cleaning. We then described how to use **fitted Gaussian distributions** to identify outliers. This lead us to a discussion about the choice of **thresholds** and **hyperparameters**, where we went over a few different realistic scenarios. We then used a Gaussian distribution to remove terminator images from an emoji faces dataset. We learned how learning algorithms **fail** and that data scientists must know how to **debug** them. Finally, we used **PCA** to visualize our killer robot detection. \n",
    "\n",
    "\n",
    "# Resources\n",
    "\n",
    "## Core Resources\n",
    "\n",
    "- [Anomaly detection algorithm](https://www.coursera.org/lecture/machine-learning/algorithm-C8IJp)  \n",
    "Andrew Ng's limpid breakdown of anomaly detection\n",
    "\n",
    "## Additional Resources\n",
    "\n",
    "- [A review of ML techniques for anomaly detection](https://youtu.be/LRqX5uO5StA)  \n",
    "More in depth review of modern techniques for anomaly detection\n",
    "- [Anomaly Detection in sklearn](https://sdsawtelle.github.io/blog/output/week9-anomaly-andrew-ng-machine-learning-with-python.html)  \n",
    "Visual blogpost experimenting with the various outlier detection algorithms available in sklearn\n",
    "- [sklearn official documentation - outlier detection](https://scikit-learn.org/stable/modules/outlier_detection.html)  "
   ]
  }
 ],
 "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.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}