incorrect results of IsolationForest

Question

I inspired by this notebook, and I'm experimenting IsolationForest algorithm for anomaly detection context on the SF version of KDDCUP99 dataset, including 4 attributes. The data is directly fetched from sklearn and after preprocessing (label encoding the categorical feature) passed to the IF algorithm with default setup (except for n_estimator, which I used GridsearchCV to find the optimum number of trees).

The problem is based on the KDDCUP99 dataset documentation, and its User guide and my calculation in the following implementation the anomaly rate is 0.5% and it means if I set contamination=0.005, it should give me the best results but surprisingly, it does not!

The full code is as follows:

from sklearn import datasets
from sklearn import preprocessing
from sklearn.model_selection import train_test_split

from sklearn.ensemble import IsolationForest
from sklearn.metrics import classification_report, ConfusionMatrixDisplay
from sklearn.metrics import confusion_matrix
from sklearn.metrics import recall_score, silhouette_score, roc_curve, roc_auc_score, f1_score, precision_recall_curve, auc
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer
from sklearn.metrics import accuracy_score

import pandas as pd
import numpy as np
import seaborn as sns
import itertools
import matplotlib.pyplot as plt
import datetime

%matplotlib inline


def byte_decoder(val):
    # decodes byte literals to strings
    
    return val.decode('utf-8')
def plot_confusion_matrix(cm, title, classes=["Normal", "Anomaly"],
                          cmap=plt.cm.Blues, save=False, saveas="MyFigure.png"):
    
    # print Confusion matrix with blue gradient colours
    
    cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
    
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=90)
    plt.yticks(tick_marks, classes)

    fmt = '.1%'
    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, format(cm[i, j], fmt),
                horizontalalignment="center",
                color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    
    if save:
        plt.savefig(saveas, dpi=100)

#Load Dataset KDDCUP99 from sklearn
target = 'target'
sf = datasets.fetch_kddcup99(subset='SF', percent10=False)
dfSF=pd.DataFrame(sf.data, 
                  columns=["duration", "service", "src_bytes", "dst_bytes"])
assert len(dfSF)>0, "SF dataset no loaded."

dfSF[target]=sf.target
anomaly_rateSF = 1.0 - len(dfSF.loc[dfSF[target]==b'normal.'])/len(dfSF)

"SF Anomaly Rate is:"+"{:.1%}".format(anomaly_rateSF)
#'SF Anomaly Rate is: 0.5%'

#Data Processing 
toDecodeSF = ['service']
# apply hot encoding to fields of type string
# convert all abnormal target types to a single anomaly class

dfSF['binary_target'] = [1 if x==b'normal.' else -1 for x in dfSF[target]]
    
leSF = preprocessing.LabelEncoder()

for f in toDecodeSF:
    dfSF[f + " (encoded)"] = list(map(byte_decoder, dfSF[f]))
    dfSF[f + " (encoded)"] = leSF.fit_transform(dfSF[f])

for f in toDecodeSF:
  dfSF.drop(f, axis=1, inplace=True)

dfSF.drop(target, axis=1, inplace=True)

#check rate of Anomaly for setting contamination parameter in IF
dfSF["binary_target"].value_counts() / np.sum(dfSF["binary_target"].value_counts())



#train & data split
X_train_sf, X_test_sf, y_train_sf, y_test_sf = train_test_split(dfSF.drop('binary_target', axis=1), dfSF['binary_target'], 
                                                                test_size=0.33, random_state=11, stratify=dfSF['binary_target'])

clfIF = IsolationForest(max_samples="auto", random_state=11, contamination = "auto", n_estimators=100, n_jobs=-1)

clfIF.fit(X_train_sf,y_train_sf)
y_pred_train = clfIF.predict(X_train_sf)

#Results on SF training set:
print(classification_report(y_train_sf, y_pred_train, target_names=["Anomaly", "Normal"]))
print("AUC: ", "{:.1%}".format(roc_auc_score(y_train_sf, y_pred_train)))
#AUC is about 93%.

#Results on SF training set:
clfIF = IsolationForest(max_samples="auto", random_state=11, contamination = "auto", n_estimators=100, n_jobs=-1)

clfIF.fit(X_train_sf, y_train_sf)
y_pred_test = clfIF.predict(X_test_sf)
#AUC is about 93%.


scoring = {'AUC': 'roc_auc', 'Recall': make_scorer(recall_score, pos_label=-1)}
gs = GridSearchCV(IsolationForest(max_samples="auto", random_state=11, n_jobs=-1),
                  param_grid={'n_estimators': range(1, 110, 5)},
                  scoring=scoring, refit='Recall', return_train_score=True, cv=3, verbose=1, n_jobs=-1)
gs.fit(X_train_sf, y_train_sf)
results = gs.cv_results_
gs.best_estimator_
max_recall_n_estimators = pd.DataFrame(results).iloc[np.argmax(pd.DataFrame(results)["mean_test_Recall"])]["param_n_estimators"]
print(max_recall_n_estimators)



fig, ax = plt.subplots(figsize=(8, 20))
# fig.tight_layout()

for idx, cont in zip(range(1, 4), [0.005, 0.1, 0.2]):
    iso_for_ = IsolationForest(random_state=11, 
                             n_estimators=max_recall_n_estimators,  
                             max_samples="auto",  
                             contamination=cont, 
                             n_jobs=-1)
    iso_for_.fit(X_train_sf, y_train_sf)
  
    y_pred_ = iso_for_.predict(X_test_sf)

    ax = plt.subplot(1, 3, idx)
    cm = confusion_matrix(y_test_sf, y_pred_)
    cm = np.round(cm / cm.sum(axis=1)[:, np.newaxis], decimals=2)
    cmd = ConfusionMatrixDisplay(cm, display_labels=["Anomaly", "Normal"])
    cmd.plot(ax=ax, xticks_rotation=45, cmap=plt.cm.Blues, values_format=".2%")
    cmd.ax_.set_title(f"IF (contamination = {cont}) \nconfusion matrix - SF")
    cmd.im_.colorbar.remove()
    cmd.ax_.set_xlabel('')
    cmd.ax_.set_ylabel('')

fig.text(0, 0.5, "True Label", rotation=90, va='center')
fig.text(0.4, 0.4, 'Predicted Label', ha='left')
plt.subplots_adjust(wspace=0.9)
plt.suptitle("Confusion matrix of sklearn IF for various contaminations (SF)", y=0.63)
plt.show()

This is the results of the confusion matrix:

Probably I am missing something here, and any help will be highly appreciated.

Edit: We know that IF algorithm label anomaly labels based on outliers score results and try to plot them on train and test set for better realization as follows:

outlier_scores_train = iso_for_sf.decision_function(X_train_sf)
outlier_scores_test  = iso_for_sf.decision_function(X_test_sf)

train_df = pd.concat([X_train_sf.reset_index(drop=True), pd.DataFrame(outlier_scores_train, columns=["outlier score"])], axis=1)
test_df = pd.concat([X_test_sf.reset_index(drop=True), pd.DataFrame(outlier_scores_test, columns=["outlier score"])], axis=1)

sns.distplot(outlier_scores_train, kde=True, norm_hist=True, color="b", label="Trainset outlier scores")
sns.distplot(outlier_scores_test,  kde=True, norm_hist=True, color="r", label="Testset outlier scores", hist_kws={'alpha': 0.3})

plt.xlabel("outlier score")
plt.ylabel("frequency")
plt.title("Outlier scores Distribution")

plt.legend()
plt.show()

Answer 1

Contamination in isolation forest is basically a hyperparameter. It is used to pick the threshold for the score to distinguish anomalies from non-anomalies. It should be close to the proportion of anomalies assuming that your training data is perfectly representative to the test (real!) data and that the algorithm was able to learn the representation of the data reasonably well for this task. As you can see, those are quite a big assumptions to be met.

This may be more obvious with simpler algorithm like EllipticEnvelope, that fits a multivariate Gaussian to your data and cuts-off the values from low density regions based on the value of the contamination hyperparameter. It would work perfectly if you set the value of contamination to the proportion of anomalous samples if your data is distributed as a multivariate Gaussian, but gives no such guarantees for other kinds of data. Isolation forest is a bit more complicated, but same problems apply.

TL;DR Just treat it as a hyperparameter.