13 - Anomaly Detection/anomaly_detection_class-empty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/80x15.png\" /></a>&nbsp;| Fabrice Jimenez | <a href=\"https://supaerodatascience.github.io/machine-learning/\">https://supaerodatascience.github.io/machine-learning/</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Detection - Practical Follow-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is following the progression of the Anomaly Detection class. It provides practical illustrations in Python and short exercises to understand the notions we have seen in this course."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Author: Fabrice JIMENEZ\n",
    "    \n",
    "Link to course materials: https://github.com/SupaeroDataScience/machine-learning/tree/main/13%20-%20Anomaly%20Detection\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminary loading with Google Colab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using this notebook with Google Colab, please execute first the following cells, to retrieve the GitHub repository content. Otherwise, ignore these cells and move to the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.mkdir('data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -c https://github.com/SupaeroDataScience/machine-learning/blob/main/13%20-%20Anomaly%20Detection/data/X1.npy?raw=true -O data/X1.npy\n",
    "!wget -c https://github.com/SupaeroDataScience/machine-learning/blob/main/13%20-%20Anomaly%20Detection/data/X1_1.npy?raw=true -O data/X1_1.npy\n",
    "!wget -c https://github.com/SupaeroDataScience/machine-learning/blob/main/13%20-%20Anomaly%20Detection/data/X2.npy?raw=true -O data/X2.npy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Loading Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is sampled from random gaussian distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [10, 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = np.load('data/X1.npy')\n",
    "X1_1 = np.load('data/X1_1.npy') # Variable X1 with far outlier\n",
    "X2 = np.load('data/X2.npy')\n",
    "df = pd.DataFrame()\n",
    "df['X1'] = X1\n",
    "df['X1_1'] = X1_1\n",
    "df['X2'] = X2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Don't worry about these functions at the moment, they will be used later to plot some results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotAnomalies1D(s, anomalies, threshold1, threshold2):\n",
    "    \"\"\"\n",
    "        s: Pandas Series containing all the points to plot\n",
    "        anomalies: Pandas Series containing all the points which are anomalies\n",
    "        threshold1: Float value - minimum threshold to be normal\n",
    "        threshold2: Float value - maximum threshold to be normal\n",
    "    \"\"\"\n",
    "    plt.figure()\n",
    "    plt.plot(s, [0]*len(s), 'bo')\n",
    "    plt.plot(anomalies, [0]*len(anomalies), 'ro')\n",
    "    plt.plot([threshold1]*2, [-1,1], 'g--')\n",
    "    plt.plot([threshold2]*2, [-1,1], 'g--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotAnomalies2D(df, clf_name, clf):\n",
    "    \"\"\"\n",
    "        df: Pandas DataFrame containing all the points to plot (for features X1 and X2)\n",
    "        clf_name: String value - name of the outlier detection model\n",
    "        clf: Scikit Learn model instance - the trained outlier detection model\n",
    "    \"\"\"\n",
    "    plt.figure()\n",
    "    plt.plot(df['X1'],df['X2'], 'o')\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('X2')\n",
    "    plt.xlim([df['X1'].min()-3,df['X1'].max()+3])\n",
    "    plt.ylim([df['X2'].min()-3,df['X2'].max()+3])\n",
    "    plt.title(clf_name)\n",
    "    \n",
    "    if clf_name == 'LOF':\n",
    "        ypred = clf.fit_predict(df[['X1','X2']])\n",
    "        plt.plot(df['X1'][ypred==-1],df['X2'][ypred==-1],'ro')\n",
    "    else:\n",
    "        xx, yy = np.meshgrid(np.linspace(df['X1'].min()-3,df['X1'].max()+3, 500), np.linspace(df['X2'].min()-3,df['X2'].max()+3, 500))\n",
    "        Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
    "        Z = Z.reshape(xx.shape)\n",
    "        plt.contour(xx, yy, Z, levels=[0], linewidths=2, colors='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotAnomalyScore2D(df, clf_name, clf):\n",
    "    \"\"\"\n",
    "        df: Pandas DataFrame containing all the points to plot (for features X1 and X2)\n",
    "        clf_name: String value - name of the outlier detection model\n",
    "        clf: Scikit Learn model instance - the trained outlier detection model\n",
    "    \"\"\"\n",
    "    if clf_name == 'LOF':\n",
    "        score = clf.negative_outlier_factor_\n",
    "    else:\n",
    "        score = clf.decision_function(df[['X1','X2']])\n",
    "    \n",
    "    plt.figure()\n",
    "    sc = plt.scatter(x=df['X1'],y=df['X2'], c=-score, cmap='Reds')\n",
    "    plt.colorbar(sc, label='anomaly score')\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('X2')\n",
    "    plt.title(clf_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Univariate Outlier Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, let's consider a single variable x = X1. Let's visualize the dataset in different ways..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 'X1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scatter plot\n",
    "plt.figure()\n",
    "plt.plot(df[x], [0]*len(df), 'bo')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histogram\n",
    "df[x].hist(bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Key statistical indicators\n",
    "df[x].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 1: if you had to detect outliers in this 1D dataset, how would you do?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement your own 1D outlier detection method. To do this, you need to compute:\n",
    "- <em>threshold1</em> : Float value - minimum threshold to be normal\n",
    "- <em>threshold2</em> : Float value - maximum threshold to be normal\n",
    "- <em>anomalies</em> : Pandas Series containing all the points which are anomalies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###### CODE HERE ######\n",
    "\n",
    "#######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotAnomalies1D(df[x], anomalies, threshold1, threshold2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 2: What is the impact of far outliers on the method(s) you used?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your method(s) on variable <em>X1_1</em> instead of <em>X1</em> : 1 single point was changed, to become a very large value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Outlier Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, let's consider 2 variables X1 and X2. Let's visualize the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(df['X1'],df['X2'], 'bo')\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will experiment here with the 4 different outlier detection methods we have just seen in class. These methods are directly available in Scikit-Learn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 3: What is the impact of hyperparameters on the detection boundaries / anomalies detected, for each method?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Question 3.a:</b> For example, compare the behavior of the different models, with contamination rate at same value. What models appear to be easy to tune (without taking into account contamination rate)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Question 3.b:</b> For example, compare the behavior of different contamination rates, with the same model. What difficulties are you facing?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.covariance import EllipticEnvelope\n",
    "from sklearn.ensemble import IsolationForest\n",
    "from sklearn.neighbors import LocalOutlierFactor\n",
    "from sklearn.svm import OneClassSVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_name = 'Elliptic Envelope'\n",
    "clf = EllipticEnvelope(contamination=0.15)\n",
    "clf.fit(df[['X1','X2']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_name = 'Isolation Forest'\n",
    "clf = IsolationForest(behaviour='new', n_estimators=100, contamination=0.15)\n",
    "clf.fit(df[['X1','X2']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_name = 'Local Outlier Factor'\n",
    "clf = LocalOutlierFactor(n_neighbors=5, contamination=0.15, novelty=True)\n",
    "clf.fit(df[['X1','X2']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_name = 'One Class SVM'\n",
    "clf = OneClassSVM(nu=0.15, kernel=\"rbf\", gamma=0.3) # nu corresponds to contamination\n",
    "clf.fit(df[['X1','X2']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotAnomalies2D(df, clf_name, clf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 4: What alternative do you have instead of tuning the contamination rate?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute and visualize, for each model, the continuous anomaly scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotAnomalyScore2D(df, clf_name, clf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Novelty Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we split the dataset to have on one side the normal points, and on the other side the new points which are anomalies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normal_df = df.iloc[:300].copy()\n",
    "new_anomalies = df.iloc[300:].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(normal_df['X1'], normal_df['X2'], 'bo')\n",
    "plt.plot(new_anomalies['X1'], new_anomalies['X2'], 'ro')\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of novelty detection is to learn the behavior of normal points. We can use unsupervised methods such as those we have seen earlier, but we can also use supervised learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 5: Implement your novelty detection approach using supervised learning, to predict each variable in function of the others.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do this, you have to select a supervised model $\\tilde{f}$ to learn the relationship between X1 and X2. \n",
    "\n",
    "$$\n",
    "\\tilde{X}_2 = \\tilde{f}\\left(X_1\\right)\n",
    "$$You need to compute:\n",
    "- <em>ypred_normal</em> : an array containing the predicted values for the normal dataset\n",
    "- <em>ypred_anomalies</em> : an array containing the predicted values for the new anomalies dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(normal_df['X1']).reshape(-1,1)\n",
    "y = np.array(normal_df['X2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###### CODE HERE ######\n",
    "\n",
    "#######################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the superposition of actual values and the predicted values, as scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Superposition\n",
    "plt.figure()\n",
    "plt.scatter(x=normal_df['X1'],y=normal_df['X2'],c='b')\n",
    "plt.scatter(x=normal_df['X1'],y=ypred_normal,c='cyan')\n",
    "\n",
    "plt.scatter(x=new_anomalies['X1'],y=new_anomalies['X2'],c='r')\n",
    "plt.scatter(x=new_anomalies['X1'],y=ypred_anomalies,c='magenta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at the predicted values versus the actual values, for normal and anomaly points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Value / Value plot\n",
    "plt.figure()\n",
    "plt.scatter(x=normal_df['X2'], y=ypred_normal, c='b')\n",
    "plt.scatter(x=new_anomalies['X2'], y=ypred_anomalies,c='r')\n",
    "plt.plot([-50,50], [-50,50], 'g--')\n",
    "plt.xlim([-40,40])\n",
    "plt.ylim([-21,21])\n",
    "plt.xlabel('Actual values')\n",
    "plt.ylabel('Predicted values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Question 6: To compute an anomaly score, we will use the prediction error of the model. Compute for example the mean squared error for normal points and anomaly points.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###### CODE HERE ######\n",
    "\n",
    "#######################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's visualize the prediction error as color levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "sc = plt.scatter(x=list(normal_df['X1']) + list(new_anomalies['X1']), \n",
    "            y=list(normal_df['X2']) + list(new_anomalies['X2']), \n",
    "            c=list(error_normal) + list(error_anomalies),\n",
    "            cmap='Reds')\n",
    "plt.colorbar(sc, label='Prediction error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like a pretty good anomaly score!"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}