Merge branch 'renku/autosave/solange.emmenegger@hslu.ch/master/0b179d7c/5e45701c' into 'master'

Remove Ng Video

See merge request solange.emmenegger/ml-adml-hslu!2

notebooks/11B Anomaly Detection/Anomaly Detection - Part 1 - Multivariate Gaussian.ipynb

 %% Cell type:markdown id: tags:
 
 # Anomaly Detection
 
 %% Cell type:code id: tags:
 
 ``` python
 %matplotlib inline
 
 import numpy as np
 import pandas as pd
 import seaborn as sns
 import matplotlib.pyplot as plt
 import matplotlib as mpl
 
 from scipy.stats import multivariate_normal
 from sklearn.preprocessing import StandardScaler
 from sklearn.model_selection import train_test_split
 
 from sklearn.metrics import precision_recall_fscore_support
 from sklearn.metrics import accuracy_score
 from sklearn.metrics import confusion_matrix
 from sklearn.metrics import f1_score
 from sklearn.metrics.pairwise import euclidean_distances
 
 from sklearn.cluster import KMeans
 
 import matplotlib.gridspec as gridspec
 
 from tqdm.notebook import tqdm
 import ipywidgets as widgets
 
 import warnings
 warnings.filterwarnings("ignore", category=FutureWarning)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Exercise 1 -  Multivariate Gaussian distribution
-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 id: tags:
-
-``` python
-from IPython.display import YouTubeVideo
-YouTubeVideo('g2YBWQnqOpw')
-```
+In the lecture we discussed anomaly detection using the univariate 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). In this exercise we will guide you through the implementation of fraud detetion using anomaly detection with the multivariate Gaussian distribution.
 
 %% Cell type:markdown id: tags:
 
 ### Data Quality Assessment
 
 %% Cell type:code id: tags:
 
 ``` python
 df = pd.read_csv('creditcard.csv')
 df.head()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 df.shape
 ```
 
 %% Cell type:markdown id: tags:
 
 > Check for null values
 
 %% Cell type:code id: tags:
 
 ``` python
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 df.isnull().any().any()
 ```
 
 %% Cell type:markdown id: tags:
 
 #### Split fraudulent data from genuine data
 We start by separating the genuine data from the fraudulent data
 
 %% Cell type:code id: tags:
 
 ``` python
 # genuine = ...
 # anomalies = ...
 # print("Genuine Transactions:", genuine.shape[0])
 # print("Anomalous Transaction:", anomalies.shape[0])
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 genuine = df[df.Class == 0]
 anomalies = df[df.Class == 1]
 print("Genuine Transactions:", genuine.shape[0])
 print("Anomalous Transaction:", anomalies.shape[0])
 ```
 
 %% Cell type:markdown id: tags:
 
 #### Amount feature
 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 id: tags:
 
 ``` python
 genuine.Amount.loc[genuine.Amount < 500].hist(bins=100)
 print("Mean", genuine.Amount.mean())
 print("Median", genuine.Amount.median())
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 anomalies.Amount.loc[anomalies.Amount < 500].hist(bins=100)
 print("Mean", anomalies.Amount.mean())
 print("Median", anomalies.Amount.median())
 ```
 
 %% Cell type:markdown id: tags:
 
 * 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 id: tags:
 
 ### Feature Engineering
 
 %% Cell type:markdown id: tags:
 
 #### Distribution of the features
 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 id: tags:
 
 ``` python
 feature_cols = genuine.drop(columns=["Class"]).columns
 plt.figure(figsize=(12, len(feature_cols)*4))
 gs = gridspec.GridSpec(len(feature_cols), 2)
 
 for i, col in enumerate(tqdm(feature_cols)):
     ax = plt.subplot(gs[i])
     sns.distplot(genuine[col], color='g',label='Genuine')
     sns.distplot(anomalies[col], color='r',label='Anomaly')
     ax.legend()
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 > Based on the plots we made, drop the features that might not be useful.
 
 %% Cell type:code id: tags:
 
 ``` python
 # features_to_drop = [...]
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 features_to_drop = ['Time', 'V8', 'V13', 'V15', 'V20', 'V22', 'V23', 'V24', 'V25', 'V26', 'V27', 'V28']
 ```
 
 %% Cell type:markdown id: tags:
 
 * We can see that the distribution of the fraudulent transaction (class=1) is matching with the genuine transactions (class=0) for the following features:
  * *V8*, *V13*, *V15*, *V20*, *V22*, *V23*, *V24*, *V25*, *V26*, *V27*, *V28*
  * Therefore they might not be useful in finding fradulent transactions.
 * 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 id: tags:
 
 Let's drop these features in our splitted dataset
 
 %% Cell type:code id: tags:
 
 ``` python
 genuine = genuine.drop(features_to_drop, axis=1)
 anomalies = anomalies.drop(features_to_drop, axis=1)
 ```
 
 %% Cell type:markdown id: tags:
 
 #### Split the data
 We split the targets from the features
 
 %% Cell type:code id: tags:
 
 ``` python
 genuine_X = genuine.drop(columns=["Class"]).values
 genuine_y = genuine.Class.values
 
 anomalies_X = anomalies.drop(columns=["Class"]).values
 anomalies_y = anomalies.Class.values
 ```
 
 %% Cell type:markdown id: tags:
 
 Now we use 60% of the *genuine* data as the training set.
 
 %% Cell type:code id: tags:
 
 ``` python
 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 id: tags:
 
 And add the anomalies to the test set
 
 %% Cell type:code id: tags:
 
 ``` python
 # Add anomalies to the test set
 X_test = np.concatenate([X_test, anomalies_X])
 y_test = np.concatenate([y_test, anomalies_y])
 
 print(X_test.shape)
 print(y_test.shape)
 print("Number of anomalies:", y_test.sum())
 ```
 
 %% Cell type:markdown id: tags:
 
 *(Note that because we further split the data we don't have to shuffle it)*
 
 %% Cell type:markdown id: tags:
 
 Now we can further split the test data into a test set and a validation set.
 
 %% Cell type:code id: tags:
 
 ``` python
 X_test, X_val, y_test, y_val = train_test_split(X_test, y_test, test_size=0.5, random_state=42)
 print(X_test.shape)
 print(y_test.shape)
 print("Number of anomalies y_test:", y_test.sum())
 print(X_val.shape)
 print(y_val.shape)
 print("Number of anomalies y_val:", y_val.sum())
 ```
 
 %% Cell type:markdown id: tags:
 
 #### Data Normalization
 Now that we have splitted our data, we can normalize it by applying z-normalization.
 
 > Normalize the data
 
 %% Cell type:code id: tags:
 
 ``` python
 # scaler = ...
 # X_train = ...
 # X_test = ...
 # X_val = ...
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 scaler = StandardScaler()
 X_train = scaler.fit_transform(X_train)
 X_test = scaler.transform(X_test)
 X_val = scaler.transform(X_val)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Baseline model
 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
 
 > Implement the dummy predictor
 
 %% Cell type:code id: tags:
 
 ``` python
 def dummy_predict(X):
     # START YOUR CODE
     pass
     # END YOUR CODE
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 def dummy_predict(X):
     return np.zeros(X.shape[0])
 ```
 
 %% Cell type:markdown id: tags:
 
 > 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 id: tags:
 
 ``` python
 # Calculate the accuracy on the test set
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 y_pred = dummy_predict(X_test)
 accuracy_score(y_pred, y_test)
 ```
 
 %% Cell type:markdown id: tags:
 
 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 id: tags:
 
 Okay, let's be honest there must be a mistake.
 
 > Let's check how many fraudulent and genuine data we have.
 
 %% Cell type:code id: tags:
 
 ``` python
 # Calculate the percentage of fraudulent data
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 print(df.Class.value_counts())
 print("Percentage of fraudulent data:", 100*df.Class.value_counts()[1] /len(df))
 ```
 
 %% Cell type:markdown id: tags:
 
 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 id: tags:
 
 ``` python
 def plot_confusion_matrix(cm):
     fig, (ax1) = plt.subplots(ncols=1, figsize=(5,5))
     sns.heatmap(cm,
                 xticklabels=['Genuine', 'Fraud'],
                 yticklabels=['Genuine', 'Fraud'],
                 annot=True,ax=ax1,
                 linewidths=.2,linecolor="Darkblue", cmap="Blues")
     plt.title('Confusion Matrix', fontsize=14)
     plt.show()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm = confusion_matrix(y_test, y_pred)
 plot_confusion_matrix(cm)
 ```
 
 %% Cell type:markdown id: tags:
 
 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 id: tags:
 
 ### Precision, Recall and F1 Score
 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 id: tags:
 
 ``` python
 precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average="binary")
 print("Precision:", precision)
 print("Recall", recall)
 print("F1", f1)
 ```
 
 %% Cell type:markdown id: tags:
 
 Now it's clear that our dummy anomaly detection model performs really bad! So let's stick to those metrics.
 
 %% Cell type:markdown id: tags:
 
 ### Fit the multivariate gaussian distribution
 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 id: tags:
 
 > 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 id: tags:
 
 ``` python
 def fit(X):
     pass
     # START YOUR CODE HERE
     #return gauss
     # END YOUR CODE HERE
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 def fit(X):
     mu = np.mean(X, axis=0)
     cov = np.cov(X.T)
     gauss = multivariate_normal(mean=mu, cov=cov)
     return gauss
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 gauss = fit(X_train)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Predict anomalies
 Let's use our fitted gaussian distribution to predict anomalies.
 
 > 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 id: tags:
 
 ``` python
 def predict(X, gauss, eps):
     # START YOUR CODE
     pass
     # END YOUR CODE
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 def predict(X, gauss, eps):
     densities = gauss.pdf(X)
     y_pred = densities < eps
     return y_pred
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 eps = 1e-20
 y_pred = predict(X_test, gauss, eps)
 precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average="binary")
 print("Precision:", precision)
 print("Recall", recall)
 print("F1", f1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm = confusion_matrix(y_test, y_pred)
 plot_confusion_matrix(cm)
 ```
 
 %% Cell type:markdown id: tags:
 
 * We have a relatively high recall value, which means that our model does detect most of the anomalies.
 * 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.
 * 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 id: tags:
 
 ### Hyperparameter tuning
 Our threshold `eps` is the hyperparameter we would like to tune.
 
 %% Cell type:markdown id: tags:
 
 But first we take a look at our densities
 
 %% Cell type:code id: tags:
 
 ``` python
 densities = gauss.pdf(X_val)
 print("Min density:", densities.min())
 print("Max density:", densities.max())
 ```
 
 %% Cell type:markdown id: tags:
 
 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`.
 Let's check the values of our log transformed densities.
 
 %% Cell type:code id: tags:
 
 ``` python
 densities = gauss.logpdf(X_val)
 print("Min density:", densities.min())
 print("Max density:", densities.max())
 ```
 
 %% Cell type:markdown id: tags:
 
 These values look much better. However in order to reuse our existing code we need to redefine our `predict` function.
 
 > Reimplement the `predict` function but this time use `logpdf` instead of the `pdf` function.
 
 %% Cell type:code id: tags:
 
 ``` python
 def predict(X, gauss, eps):
     # YOUR CODE
     pass
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 def predict(X, gauss, eps):
     densities = gauss.logpdf(X)
     y_pred = densities < eps
     return y_pred
 ```
 
 %% Cell type:markdown id: tags:
 
 #### Fit our threshold by sorting
 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 id: tags:
 
 ``` python
 expected_anomalies = y_val.sum()
 sorted_densities = np.sort(gauss.logpdf(X_val))
 best_eps = (sorted_densities[expected_anomalies] +
             0.5 * (sorted_densities[expected_anomalies + 1] -
                    sorted_densities[expected_anomalies]))
 
 
 y_pred = predict(X_val, gauss, best_eps)
 f1 = f1_score(y_val, y_pred, average="binary")
 print("Best epsilon {} with f1 score {}".format(best_eps, f1))
 ```
 
 %% Cell type:markdown id: tags:
 
 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 id: tags:
 
 #### Grid Search Approach
 
 %% Cell type:markdown id: tags:
 
 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 id: tags:
 
 ``` python
 densities.min()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 X_anomalies =  anomalies.drop(columns=["Class"]).values
 densities = gauss.logpdf(X_val)
 
 max_eps = densities.max()
 print("Maximum epsilon:", max_eps)
 
 min_eps = -500
 print("Minimum epsilon:", min_eps)
 ```
 
 %% Cell type:markdown id: tags:
 
 This means the optimal value of `epsilon` lies between -500 and -16. Let's try to find it by using a grid search approach.
 > Loop through the list of epsilons and calculate for each epsilon the *precision*, *recall* and *f1* score.
 
 %% Cell type:code id: tags:
 
 ``` python
 f1_scores = []
 recall_scores = []
 precision_scores = []
 
 epsilons = np.linspace(min_eps, max_eps, num=1000)
 
 # START YOUR CODE
 #for eps in tqdm(epsilons):
 
 #best_eps = epsilons[np.argmax(f1_scores)]
 #best_f1 = np.max(f1_scores)
 
 #print("Best epsilon {} with f1 score {}".format(best_eps, best_f1))
 # END YOUR CODE
 ```
 
 %% Cell type:markdown id: tags:
 
 *Click on the dots to display the solution*
 
 %% Cell type:code id: tags:
 
 ``` python
 f1_scores = []
 recall_scores = []
 precision_scores = []
 
 epsilons = np.linspace(min_eps, max_eps, num=1000)
 densities = gauss.logpdf(X_val)
 for eps in tqdm(epsilons):
     y_pred = densities < eps
     precision, recall, f1, _ = precision_recall_fscore_support(y_val, y_pred, average="binary")
     precision_scores.append(precision)
     recall_scores.append(recall)
     f1_scores.append(f1)
 
 best_eps_gridsearch = epsilons[np.argmax(f1_scores)]
 best_f1 = np.max(f1_scores)
 
 print("Best epsilon {} with f1 score {}".format(best_eps_gridsearch, best_f1))
 ```
 
 %% Cell type:markdown id: tags:
 
 ##### Visualization
 Let's visualize the *F1*, *Precision* and *Recall* curves
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, ax = plt.subplots(figsize=(10,8))
 ax.plot(epsilons, f1_scores, label="f1")
 ax.plot(epsilons, precision_scores, label="precision")
 ax.plot(epsilons, recall_scores, label="recall")
 ax.plot(best_eps, best_f1, marker="o")
 ax.set_xlabel("Epsilon")
 ax.set_ylabel("Score")
 ax.legend()
 ax.set_xlim(max_eps, min_eps)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Evaluate the model on the test set
 
 %% Cell type:markdown id: tags:
 
 As we have selected an epsilon, we can check how good our model performs on the test set.
 
 %% Cell type:code id: tags:
 
 ``` python
 y_pred = predict(X_test, gauss, best_eps)
 precision, recall, f1, _ = precision_recall_fscore_support(y_test, y_pred, average="binary")
 print("Precision:", precision)
 print("Recall", recall)
 print("F1", f1)
 ```
 
 %% Cell type:markdown id: tags:
 
 We plot the confusion matrix again
 
 %% Cell type:code id: tags:
 
 ``` python
 cm = confusion_matrix(y_test, y_pred)
 plot_confusion_matrix(cm)
 ```