Anomaly Detection - Part 1 - Multivariate Gaussian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "from scipy.stats import multivariate_normal\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "from sklearn.metrics import precision_recall_fscore_support\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import f1_score\n",
    "from sklearn.metrics.pairwise import euclidean_distances\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1 -  Multivariate Gaussian distribution\n",
    "You watched the video tutorial by Andrew Ng on anomaly detection using the multivariate Gaussian distribution. There is a data set on ILIAS containing credit card transactions made available by a European bank in September 2013 on [Kaggle](https://www.kaggle.com/mlg-ulb/creditcardfraud). Not surprisingly, the data set is anonymized, i.e. a PCA transformation was executed on all features (each feature therefore is a linear combination of the original but unknown features). The only exceptions are amount and time (milliseconds since first transaction). Finally, there is a class feature to indicate whether the transaction is fraudulent (value: 1) or genuine (value: 0). Implement fraud detection using the statistical approach introduced by Andrew Ng."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('g2YBWQnqOpw')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Quality Assessment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('creditcard.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check for null values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isnull().any().any()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Split fraudulent data from genuine data\n",
    "We start by separating the genuine data from the fraudulent data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# genuine = ...\n",
    "# anomalies = ...\n",
    "# print(\"Genuine Transactions:\", genuine.shape[0])\n",
    "# print(\"Anomalous Transaction:\", anomalies.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "genuine = df[df.Class == 0]\n",
    "anomalies = df[df.Class == 1]\n",
    "print(\"Genuine Transactions:\", genuine.shape[0])\n",
    "print(\"Anomalous Transaction:\", anomalies.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Amount feature\n",
    "Fortunately, that *Amount* feature has not been anonymized with the PCA transformation. So let's take a look at the *Amount* feature and see if we can find a pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genuine.Amount.loc[genuine.Amount < 500].hist(bins=100)\n",
    "print(\"Mean\", genuine.Amount.mean())\n",
    "print(\"Median\", genuine.Amount.median())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anomalies.Amount.loc[anomalies.Amount < 500].hist(bins=100)\n",
    "print(\"Mean\", anomalies.Amount.mean())\n",
    "print(\"Median\", anomalies.Amount.median())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Interestingly the median is lower but the mean is higher for the fraudulent transactions. This suggests there are some high value oriented criminals and some that focus on withdrawals \"below the radar\" to avoid detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distribution of the features\n",
    "We want to model our anomaly detection by fitting a multivariate gaussian distribution. Therefore it makes sense to plot the distribution of our features. We plot the genuine and the fraudulent transactions separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "feature_cols = genuine.drop(columns=[\"Class\"]).columns\n",
    "plt.figure(figsize=(12, len(feature_cols)*4))\n",
    "gs = gridspec.GridSpec(len(feature_cols), 2)\n",
    "\n",
    "for i, col in enumerate(tqdm(feature_cols)):\n",
    "    ax = plt.subplot(gs[i])\n",
    "    sns.distplot(genuine[col], color='g',label='Genuine')\n",
    "    sns.distplot(anomalies[col], color='r',label='Anomaly')\n",
    "    ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Based on the plots we made, drop the features that might not be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# features_to_drop = [...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "features_to_drop = ['Time', 'V8', 'V13', 'V15', 'V20', 'V22', 'V23', 'V24', 'V25', 'V26', 'V27', 'V28']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We can see that the distribution of the fraudulent transaction (class=1) is matching with the genuine transactions (class=0) for the following features:\n",
    " * *V8*, *V13*, *V15*, *V20*, *V22*, *V23*, *V24*, *V25*, *V26*, *V27*, *V28*\n",
    " * Therefore they might not be useful in finding fradulent transactions. \n",
    "* The variable *Time* is also not a useful variable since it contains the seconds elapsed between the transaction for that record and the first transaction in the dataset. The data is in increasing order, which can lead to a singular covariance matrix. Also it is not normally distributed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's drop these features in our splitted dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genuine = genuine.drop(features_to_drop, axis=1)\n",
    "anomalies = anomalies.drop(features_to_drop, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Split the data\n",
    "We split the targets from the features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genuine_X = genuine.drop(columns=[\"Class\"]).values\n",
    "genuine_y = genuine.Class.values\n",
    "\n",
    "anomalies_X = anomalies.drop(columns=[\"Class\"]).values\n",
    "anomalies_y = anomalies.Class.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use 60% of the *genuine* data as the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(genuine_X, genuine_y, test_size=0.4, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And add the anomalies to the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add anomalies to the test set\n",
    "X_test = np.concatenate([X_test, anomalies_X])\n",
    "y_test = np.concatenate([y_test, anomalies_y])\n",
    "\n",
    "print(X_test.shape)\n",
    "print(y_test.shape)\n",
    "print(\"Number of anomalies:\", y_test.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*(Note that because we further split the data we don't have to shuffle it)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can further split the test data into a test set and a validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test, X_val, y_test, y_val = train_test_split(X_test, y_test, test_size=0.5, random_state=42)\n",
    "print(X_test.shape)\n",
    "print(y_test.shape)\n",
    "print(\"Number of anomalies y_test:\", y_test.sum())\n",
    "print(X_val.shape)\n",
    "print(y_val.shape)\n",
    "print(\"Number of anomalies y_val:\", y_val.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data Normalization\n",
    "Now that we have splitted our data, we can normalize it by applying z-normalization.\n",
    "\n",
    "> Normalize the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# scaler = ...\n",
    "# X_train = ...\n",
    "# X_test = ...\n",
    "# X_val = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_test = scaler.transform(X_test)\n",
    "X_val = scaler.transform(X_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline model\n",
    "To get an idea on how good our model performs, let's implement a dummy predictor. Our dummy model should mark each transaction as genuine\n",
    "\n",
    "> Implement the dummy predictor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2,
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "def dummy_predict(X):\n",
    "    # START YOUR CODE\n",
    "    pass\n",
    "    # END YOUR CODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def dummy_predict(X):\n",
    "    return np.zeros(X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Now predict on the test set and calculate the accuracy by using the Scikit-Learn function [accuracy_score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Calculate the accuracy on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "y_pred = dummy_predict(X_test)\n",
    "accuracy_score(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow, such a high accuracy! It seems like our work is done and we can spend the rest of our day drinking coffee!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, let's be honest there must be a mistake. \n",
    "\n",
    "> Let's check how many fraudulent and genuine data we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Calculate the percentage of fraudulent data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "print(df.Class.value_counts())\n",
    "print(\"Percentage of fraudulent data:\", 100*df.Class.value_counts()[1] /len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our dataset is highly imbalanced! We could create a balanced test set and then calculate the accuracy. But there must be another way... Let's plot the confusion matrix first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm):\n",
    "    fig, (ax1) = plt.subplots(ncols=1, figsize=(5,5))\n",
    "    sns.heatmap(cm, \n",
    "                xticklabels=['Genuine', 'Fraud'],\n",
    "                yticklabels=['Genuine', 'Fraud'],\n",
    "                annot=True,ax=ax1,\n",
    "                linewidths=.2,linecolor=\"Darkblue\", cmap=\"Blues\")\n",
    "    plt.title('Confusion Matrix', fontsize=14)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_test, y_pred)\n",
    "plot_confusion_matrix(cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we have a lot of *True Positives* and only few *True Negatives*. This is due to the class imbalance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision, Recall and F1 Score\n",
    "For a skewed data set, we often use other metrics like [Precision](https://en.wikipedia.org/wiki/Precision_and_recall), [Recall](https://en.wikipedia.org/wiki/Precision_and_recall) or the harmonic mean of both, called [F1 Score](https://en.wikipedia.org/wiki/F1_score). Scikit-Learn provides the function [precision_recall_fscore_support](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.precision_recall_fscore_support.html) which calculates all these metrics at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average=\"binary\")\n",
    "print(\"Precision:\", precision)\n",
    "print(\"Recall\", recall)\n",
    "print(\"F1\", f1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's clear that our dummy anomaly detection model performs really bad! So let's stick to those metrics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the multivariate gaussian distribution\n",
    "Now that we have splitted our data set into a training, test and cross validation set, we can train our anomaly detection model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Complete the `fit` function which calculate the *mean* and the *covariance* matrix from the training set and then fits a multivariate gaussian distribution. Use the [multivariate_normal](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.multivariate_normal.html) function from SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2,
    "solution2": "shown",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "def fit(X):\n",
    "    pass\n",
    "    # START YOUR CODE HERE\n",
    "    #return gauss\n",
    "    # END YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "shown"
   },
   "outputs": [],
   "source": [
    "def fit(X):\n",
    "    mu = np.mean(X, axis=0)\n",
    "    cov = np.cov(X.T)\n",
    "    gauss = multivariate_normal(mean=mu, cov=cov)\n",
    "    return gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss = fit(X_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predict anomalies\n",
    "Let's use our fitted gaussian distribution to predict anomalies. \n",
    "\n",
    "> Use the `pdf` method of the gauss model to calculate the densities. Then mark all elements which are below the threshold `epsilon` as an anomaly. We use the value 10<sup>-20</sup> as the threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "def predict(X, gauss, eps):\n",
    "    # START YOUR CODE\n",
    "    pass\n",
    "    # END YOUR CODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def predict(X, gauss, eps):\n",
    "    densities = gauss.pdf(X)\n",
    "    y_pred = densities < eps\n",
    "    return y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 1e-20\n",
    "y_pred = predict(X_test, gauss, eps)\n",
    "precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average=\"binary\")\n",
    "print(\"Precision:\", precision)\n",
    "print(\"Recall\", recall)\n",
    "print(\"F1\", f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_test, y_pred)\n",
    "plot_confusion_matrix(cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We have a relatively high recall value, which means that our model does detect most of the anomalies. \n",
    "* The low precision value  however indicates that we have a lot of *False Positives*. *False Positives* are the genuine transaction which were classified as frauds. \n",
    "* This might be a good model if for example the model is only used to determine if the user has to enter an extra verification code. However if those cards would be blocked, this would be really bad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hyperparameter tuning\n",
    "Our threshold `eps` is the hyperparameter we would like to tune. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But first we take a look at our densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "densities = gauss.pdf(X_val)\n",
    "print(\"Min density:\", densities.min())\n",
    "print(\"Max density:\", densities.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values of our densities seem to be rather small! Luckily when we are dealing with probabilities we can always apply a logarithmic transformation as it is a strictly monotonic function. Scipy already privdes an implementation `logpdf`. \n",
    "Let's check the values of our log transformed densities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "densities = gauss.logpdf(X_val)\n",
    "print(\"Min density:\", densities.min())\n",
    "print(\"Max density:\", densities.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These values look much better. However in order to reuse our existing code we need to redefine our `predict` function. \n",
    "\n",
    "> Reimplement the `predict` function but this time use `logpdf` instead of the `pdf` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "def predict(X, gauss, eps):\n",
    "    # YOUR CODE\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def predict(X, gauss, eps):\n",
    "    densities = gauss.logpdf(X)\n",
    "    y_pred = densities < eps\n",
    "    return y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit our threshold by sorting\n",
    "We can try to find the optimal threshold by sorting our densities and then picking the density with the index of the number of anomalies on our validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected_anomalies = y_val.sum()\n",
    "sorted_densities = np.sort(gauss.logpdf(X_val))\n",
    "best_eps = (sorted_densities[expected_anomalies] + \n",
    "            0.5 * (sorted_densities[expected_anomalies + 1] - \n",
    "                   sorted_densities[expected_anomalies]))\n",
    "\n",
    "\n",
    "y_pred = predict(X_val, gauss, best_eps)\n",
    "f1 = f1_score(y_val, y_pred, average=\"binary\")\n",
    "print(\"Best epsilon {} with f1 score {}\".format(best_eps, f1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow this worked out very well on our validation set! Let's check if we can still improve our result by applying a grid search approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Grid Search Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the densities for each anomaly and assume that the maximum epsilon is equal to the maximum density. We take $-500$ as our minimum `epsilon` for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "densities.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_anomalies =  anomalies.drop(columns=[\"Class\"]).values\n",
    "densities = gauss.logpdf(X_val)\n",
    "\n",
    "max_eps = densities.max()\n",
    "print(\"Maximum epsilon:\", max_eps)\n",
    "\n",
    "min_eps = -500\n",
    "print(\"Minimum epsilon:\", min_eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means the optimal value of `epsilon` lies between -500 and -16. Let's try to find it by using a grid search approach.\n",
    "> Loop through the list of epsilons and calculate for each epsilon the *precision*, *recall* and *f1* score. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "f1_scores = []\n",
    "recall_scores = []\n",
    "precision_scores = []\n",
    "\n",
    "epsilons = np.linspace(min_eps, max_eps, num=1000)\n",
    "\n",
    "# START YOUR CODE\n",
    "#for eps in tqdm(epsilons):\n",
    "\n",
    "#best_eps = epsilons[np.argmax(f1_scores)]\n",
    "#best_f1 = np.max(f1_scores)\n",
    "\n",
    "#print(\"Best epsilon {} with f1 score {}\".format(best_eps, best_f1))\n",
    "# END YOUR CODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "f1_scores = []\n",
    "recall_scores = []\n",
    "precision_scores = []\n",
    "\n",
    "epsilons = np.linspace(min_eps, max_eps, num=1000)\n",
    "densities = gauss.logpdf(X_val)\n",
    "for eps in tqdm(epsilons):\n",
    "    y_pred = densities < eps\n",
    "    precision, recall, f1, _ = precision_recall_fscore_support(y_val, y_pred, average=\"binary\")\n",
    "    precision_scores.append(precision)\n",
    "    recall_scores.append(recall)\n",
    "    f1_scores.append(f1)\n",
    "    \n",
    "best_eps_gridsearch = epsilons[np.argmax(f1_scores)]\n",
    "best_f1 = np.max(f1_scores)\n",
    "\n",
    "print(\"Best epsilon {} with f1 score {}\".format(best_eps_gridsearch, best_f1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Visualization\n",
    "Let's visualize the *F1*, *Precision* and *Recall* curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "ax.plot(epsilons, f1_scores, label=\"f1\")\n",
    "ax.plot(epsilons, precision_scores, label=\"precision\")\n",
    "ax.plot(epsilons, recall_scores, label=\"recall\")\n",
    "ax.plot(best_eps, best_f1, marker=\"o\")\n",
    "ax.set_xlabel(\"Epsilon\")\n",
    "ax.set_ylabel(\"Score\")\n",
    "ax.legend()\n",
    "ax.set_xlim(max_eps, min_eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate the model on the test set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we have selected an epsilon, we can check how good our model performs on the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = predict(X_test, gauss, best_eps)\n",
    "precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average=\"binary\")\n",
    "print(\"Precision:\", precision)\n",
    "print(\"Recall\", recall)\n",
    "print(\"F1\", f1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the confusion matrix again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_test, y_pred)\n",
    "plot_confusion_matrix(cm)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "text_representation": {
    "extension": ".py",
    "format_name": "percent",
    "format_version": "1.2",
    "jupytext_version": "0.8.6"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"