Anomaly Detection

Anomaly detection describes the process of finding outliers, for some definition of what an outler is. Since the definition of outlier is quite context dependent, there are many anomaly deteciton methods. Anomaly deteciton is generally used in an unsupervised fashion, although we use labeled data for evaluatoin.

This lecture does not focus on specific mehtods of anomaly detection. Instead, we focus on the general principles and pragmatic coding examples uisng scikit-learn. Topics covered:

• Use of PCA for anomaly detection
• Use of validation and test data sets
• Model evaluation and scoring
• ROC curve
• Precision recall curve
• Overfitting
• Example of embarassingly parallel code

%matplotlib inline

import matplotlib as mpl
import matplotlib.pyplot as plt

from scipy.io import loadmat
from scipy.io.arff import loadarff

import numpy as np

import pandas as pd

from collections import OrderedDict

from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import StratifiedKFold, train_test_split, GridSearchCV
from sklearn.metrics import precision_recall_curve, roc_curve, average_precision_score, auc, make_scorer
from sklearn.pipeline import make_pipeline
from sklearn.externals.joblib import Parallel, delayed

import seaborn as sns
sns.set(font_scale=1.5)

Information about data

The Cardiotocogrpahy dataset has the following description (from the linked URL):

The original Cardiotocography (Cardio) dataset from UCI machine learning repository consists of measurements of fetal heart rate (FHR) and uterine contraction (UC) features on cardiotocograms classified by expert obstetricians. This is a classification dataset, where the classes are normal, suspect, and pathologic. For outlier detection, The normal class formed the inliers, while the pathologic (outlier) class is downsampled to 176 points. The suspect class is discarded.

We will try to detect outliers wihtout using any lables.

Load data and inspect

data = loadmat('data/cardio.mat')

data.keys()

dict_keys(['__header__', '__version__', '__globals__', 'X', 'y'])

X = data['X']
y = data['y']

X.shape, y.shape

((1831, 21), (1831, 1))

y = y.squeeze()

Scale data to have zero mean and unit standard deviation

scaler = StandardScaler()

X = scaler.fit_transform(X)

Split into test and trainign data sets

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=20198)

Define score function

def score(a, b):
    """Normalized differnece of squares between a and b."""

    loss = np.sum((a-b)**2, axis=1)
    loss = (loss-loss.min())/(loss.max() - loss.min())
    return loss

Concept of using PCA for outlier detection

We find the robust features of the data set using PCA, expecting these to be mostly contributed by the more abundant inliers. In other words, we expect the inliers to be well summarized by the reduced dimensionality PCA data set, and the outliers to be less well summarized. Now we can reconstruct the original \(n\)-dimensional space from the \(k\)-dimensional reduced space, and compare the reconstructed data vectors with the original ones. If the reconstruction is good, we expect the distance between original and reconstructed vectors to be small, and big otherwise. Hence the normalized distance between original and reconstructed vectors might give some clue as to which data points are poorly captured by the reduced dimension PCA representation, and hence possible outliers.

def fit(n, X):
    pca = PCA(n_components=n)
    X_pca = pca.fit_transform(X)
    X_pca_inverse = pca.inverse_transform(X_pca)
    y_score = score(X, X_pca_inverse)
    return (n, y_score)

Find opitmal parameter (n) by grid search

search = Parallel(n_jobs=-1)(delayed(fit)(n, X_train) for n in range(2,X.shape[1]))

scores =  np.array([(n, average_precision_score(y_train, y_score)) for (n, y_score) in search])

best = max(scores, key=lambda x: x[1])
best

array([2.        , 0.43654539])

Visualize scores

plt.scatter(scores[:, 0], scores[:, 1])
pass

Precision-recall and ROC curves

Confusion matrix

		Truth=P
Predict=P		TP
Predict=N		FN

We will use these definitnios in the construciton of the preciison-recall and ROC curves.

PRC (Precision-recall Curve)

• X axis is the recall

  \[\frac{\text{TP}}{\text{TP + FN}}\]

• Y axis is the preciison

  \[\frac{\text{TP}}{\text{TP + FP}}\]

ROC (Receiver Operating Characteristics)

• X axis is the False Positive Rate (FPR)

  \[\frac{\text{FP}}{\text{FP + TN}}\]

• Y axis is the True Positive Rate (TPR)

  \[\frac{\text{TP}}{\text{TP + FN}}\]

To plot the PRC and ROC, we calcualte the necessary statistics for a range of cutoff values or thresholds - for example, the cutoff may be whether the score is lower or higher than 0.2, or 0.3 etc.

Imbalanced data sets

In real world settlings, imbalanced data sets, where members of one class are abundant while the other are infrequent, are common. Imbalanced data sets are nearly always due to the preponderance of true negatives (TN), for example, when screening natural compounds, nearly all will be expected to have no or minimal anti-cancer activity.

Notice that TN is NOT used in the PRC, which is hence robust to data sets imbalanced becasue most of the data come from the TN class. In contrast, the use of the FPR makes the ROC sensitive to such data sets.

Test on hold-out data

n = int(best[0])
n, y_score = fit(n, X_test)

average_precision_score(y_test, y_score)

0.4636398109078834

Visualize scores for inliers and outliers

df = pd.DataFrame(dict(score=y_score, y=y_test))

g = sns.swarmplot(data=df, x='y', y='score')
g.set_xticklabels(['Inlier', 'Outliner'])
g.set_xlabel('')
pass

Evaluate model

precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fpr, tpr, thresholds = roc_curve(y_test, y_score)

Plot ROC curve

While this looks reasoably good, it is misleading because the data set is unbalanced. Why is ROC misleading when there is unbalanced data?

plt.step(fpr, tpr, where='post')
plt.title('ROC curve')
plt.xlabel('FPR')
plt.ylabel('TPR')
pass

Plot precision-recall curve

Precision-recall is more useful for evaluating models for unbalanced data sets. Why is precision-recall robust to unbalanced data?

plt.step(recall, precision, where='post')
plt.title('Precision-Recall curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
pass

Case study

This illustrates danger of over-fitting.

Dataset information

Please download the Credit Card Fraud Detection data set and put it in the data folder. You will need to have a Kaggle account to do so.

df = pd.read_csv('data/creditcard.csv')

X = df.loc[:, df.columns.drop('Class')]
y = df.Class

y.value_counts()

0    284315
1       492
Name: Class, dtype: int64

X = scaler.fit_transform(X)

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=2018)

X_train.shape, y_train.shape

((213605, 30), (213605,))

search = Parallel(n_jobs=-1)(delayed(fit)(n, X_train) for n in range(2,X.shape[1]))

scores =  np.array([(n, average_precision_score(y_train, y_score)) for (n, y_score) in search])

This is an amazingly good average precision score

best = max(scores, key=lambda x: x[1])
best

array([27.        ,  0.69112164])

n = int(best[0])
n, y_score = fit(n, X_train)

precision, recall, thresholds = precision_recall_curve(y_train, y_score)
fpr, tpr, thresholds = roc_curve(y_train, y_score)

The precision-recall curve is also amazingly good

plt.step(recall, precision, where='post')
plt.title('Precision-Recall curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
pass

We remember that we should test on the hold-out set

n = int(best[0])
n, y_score = fit(n, X_test)

The results are markedly less impressive

average_precision_score(y_test, y_score)

0.006965112161324013

precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fpr, tpr, thresholds = roc_curve(y_test, y_score)

plt.step(recall, precision, where='post')
plt.title('Precision-Recall curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
pass

Why?

There are problay two reasons for the surprisingly bad perfromance.

• The first is taht the easy fraud cases, just by chance, ended up mostly in the trainign set
• Then grid-search found the a fortuitously effective value of \(n\)
• When we evaluate on the test set, this fails miserably because
  • The hard fraud cases are over-represented in the test cases
  • The value of \(n\) is overfitted on the training data and does not generalize

Let’s check what happeens with a differnet random number seed

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=2019)

X_train.shape, y_train.shape

((213605, 30), (213605,))

search = Parallel(n_jobs=-1)(delayed(fit)(n, X_train) for n in range(2,X.shape[1]))

scores =  np.array([(n, average_precision_score(y_train, y_score)) for (n, y_score) in search])

The average precison score is much worse

This strongly suggests that the previous result was a fluke of the distribution of the easy and hard fraud cases.

best = max(scores, key=lambda x: x[1])
best

array([10.        ,  0.19941431])

n = int(best[0])
n, y_score = fit(n, X_train)

precision, recall, thresholds = precision_recall_curve(y_train, y_score)
fpr, tpr, thresholds = roc_curve(y_train, y_score)

plt.step(recall, precision, where='post')
plt.title('Precision-Recall curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
pass