Оценка алгоритмов обнаружения выбросов

Предварительная обработка данных и обучение моделей

Разные модели обнаружения выбросов требуют различной предварительной обработки. При наличии категориальных переменных OrdinalEncoder часто является хорошей стратегией для моделей, основанных на деревьях, таких как IsolationForest, в то время как модели, основанные на ближайших соседей, такие как LocalOutlierFactor, будут подвержены влиянию порядка, индуцированного порядковым кодированием. Чтобы избежать индуцирования порядка, следует использовать OneHotEncoder.

Модели, основанные на ближайших соседей, могут также потребовать масштабирования численных признаков (см., например, Влияние масштабирования на модели k-ближайших соседей). При наличии выбросов хорошим вариантом является использование RobustScaler.

from sklearn.compose import ColumnTransformer
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (
    OneHotEncoder,
    OrdinalEncoder,
    RobustScaler,
)


def make_estimator(name, categorical_columns=None, iforest_kw=None, lof_kw=None):
    """Create an outlier detection estimator based on its name."""
    if name == "LOF":
        outlier_detector = LocalOutlierFactor(**(lof_kw or {}))
        if categorical_columns is None:
            preprocessor = RobustScaler()
        else:
            preprocessor = ColumnTransformer(
                transformers=[("categorical", OneHotEncoder(), categorical_columns)],
                remainder=RobustScaler(),
            )
    else:  # name == "IForest"
        outlier_detector = IsolationForest(**(iforest_kw or {}))
        if categorical_columns is None:
            preprocessor = None
        else:
            ordinal_encoder = OrdinalEncoder(
                handle_unknown="use_encoded_value", unknown_value=-1
            )
            preprocessor = ColumnTransformer(
                transformers=[
                    ("categorical", ordinal_encoder, categorical_columns),
                ],
                remainder="passthrough",
            )

    return make_pipeline(preprocessor, outlier_detector)

Следующая fit_predict функция возвращает средний показатель выбросов для X.

from time import perf_counter


def fit_predict(estimator, X):
    tic = perf_counter()
    if estimator[-1].__class__.__name__ == "LocalOutlierFactor":
        estimator.fit(X)
        y_pred = estimator[-1].negative_outlier_factor_
    else:  # "IsolationForest"
        y_pred = estimator.fit(X).decision_function(X)
    toc = perf_counter()
    print(f"Duration for {model_name}: {toc - tic:.2f} s")
    return y_pred

В остальной части примера мы обрабатываем по одному набору данных в каждом разделе. После загрузки данных целевые переменные изменяются, чтобы состоять из двух классов: 0 — инлиеры, 1 — выбросы. Из-за вычислительных ограничений документации scikit-learn размер выборки некоторых наборов данных уменьшается с помощью стратифицированного train_test_split.

Кроме того, мы устанавливаем n_neighbors для соответствия ожидаемому количеству аномалий expected_n_anomalies = n_samples * expected_anomaly_fraction. Это хорошая эвристика, пока доля выбросов не очень мала, поскольку n_neighbors должно быть по крайней мере больше, чем количество выборок в менее населенном кластере (см. Обнаружение выбросов с помощью Local Outlier Factor (LOF)).

import numpy as np

from sklearn.datasets import fetch_kddcup99
from sklearn.model_selection import train_test_split

X, y = fetch_kddcup99(
    subset="SA", percent10=True, random_state=42, return_X_y=True, as_frame=True
)
y = (y != b"normal.").astype(np.int32)
X, _, y, _ = train_test_split(X, y, train_size=0.1, stratify=y, random_state=42)

n_samples, anomaly_frac = X.shape[0], y.mean()
print(f"{n_samples} datapoints with {y.sum()} anomalies ({anomaly_frac:.02%})")

10065 datapoints with 338 anomalies (3.36%)

y_true = {}
y_pred = {"LOF": {}, "IForest": {}}
model_names = ["LOF", "IForest"]
cat_columns = ["protocol_type", "service", "flag"]

y_true["KDDCup99 - SA"] = y
for model_name in model_names:
    model = make_estimator(
        name=model_name,
        categorical_columns=cat_columns,
        lof_kw={"n_neighbors": int(n_samples * anomaly_frac)},
        iforest_kw={"random_state": 42},
    )
    y_pred[model_name]["KDDCup99 - SA"] = fit_predict(model, X)

Duration for LOF: 1.85 s
Duration for IForest: 0.30 s

from sklearn.datasets import fetch_covtype

X, y = fetch_covtype(return_X_y=True, as_frame=True)
s = (y == 2) + (y == 4)
X = X.loc[s]
y = y.loc[s]
y = (y != 2).astype(np.int32)

X, _, y, _ = train_test_split(X, y, train_size=0.05, stratify=y, random_state=42)
X_forestcover = X  # save X for later use

n_samples, anomaly_frac = X.shape[0], y.mean()
print(f"{n_samples} datapoints with {y.sum()} anomalies ({anomaly_frac:.02%})")

14302 datapoints with 137 anomalies (0.96%)

y_true["forestcover"] = y
for model_name in model_names:
    model = make_estimator(
        name=model_name,
        lof_kw={"n_neighbors": int(n_samples * anomaly_frac)},
        iforest_kw={"random_state": 42},
    )
    y_pred[model_name]["forestcover"] = fit_predict(model, X)

Duration for LOF: 1.84 s
Duration for IForest: 0.23 s

Набор данных Ames Housing

Набор данных Ames housing изначально является регрессионным набором данных, где целевые переменные — это цены продажи домов в Амесе, Айова. Здесь мы преобразуем его в задачу обнаружения выбросов, рассматривая дома с ценой свыше 70 USD/кв.фут. Чтобы упростить задачу, мы удаляем промежуточные цены между 40 и 70 USD/кв.фут.

import matplotlib.pyplot as plt

from sklearn.datasets import fetch_openml

X, y = fetch_openml(name="ames_housing", version=1, return_X_y=True, as_frame=True)
y = y.div(X["Lot_Area"])

# None values in pandas 1.5.1 were mapped to np.nan in pandas 2.0.1
X["Misc_Feature"] = X["Misc_Feature"].cat.add_categories("NoInfo").fillna("NoInfo")
X["Mas_Vnr_Type"] = X["Mas_Vnr_Type"].cat.add_categories("NoInfo").fillna("NoInfo")

X.drop(columns="Lot_Area", inplace=True)
mask = (y < 40) | (y > 70)
X = X.loc[mask]
y = y.loc[mask]
y.hist(bins=20, edgecolor="black")
plt.xlabel("House price in USD/sqft")
_ = plt.title("Distribution of house prices in Ames")

y = (y > 70).astype(np.int32)

n_samples, anomaly_frac = X.shape[0], y.mean()
print(f"{n_samples} datapoints with {y.sum()} anomalies ({anomaly_frac:.02%})")

2714 datapoints with 30 anomalies (1.11%)

Набор данных содержит 46 категориальных признаков. В этом случае проще использовать make_column_selector для их поиска вместо передачи вручную созданного списка.

from sklearn.compose import make_column_selector as selector

categorical_columns_selector = selector(dtype_include="category")
cat_columns = categorical_columns_selector(X)

y_true["ames_housing"] = y
for model_name in model_names:
    model = make_estimator(
        name=model_name,
        categorical_columns=cat_columns,
        lof_kw={"n_neighbors": int(n_samples * anomaly_frac)},
        iforest_kw={"random_state": 42},
    )
    y_pred[model_name]["ames_housing"] = fit_predict(model, X)

Duration for LOF: 0.88 s
Duration for IForest: 0.27 s

X, y = fetch_openml(name="cardiotocography", version=1, return_X_y=True, as_frame=False)
X_cardiotocography = X  # save X for later use
s = y == "3"
y = s.astype(np.int32)

n_samples, anomaly_frac = X.shape[0], y.mean()
print(f"{n_samples} datapoints with {y.sum()} anomalies ({anomaly_frac:.02%})")

2126 datapoints with 53 anomalies (2.49%)

y_true["cardiotocography"] = y
for model_name in model_names:
    model = make_estimator(
        name=model_name,
        lof_kw={"n_neighbors": int(n_samples * anomaly_frac)},
        iforest_kw={"random_state": 42},
    )
    y_pred[model_name]["cardiotocography"] = fit_predict(model, X)

Duration for LOF: 0.07 s
Duration for IForest: 0.18 s

Построение и интерпретация результатов

Производительность алгоритма связана с тем, насколько хорошо показатель истинно положительных результатов (TPR) при низком значении ложно положительного результата (FPR). Лучшие алгоритмы имеют кривую в верхнем левом углу графика и площадь под кривой (AUC), близкую к 1. Пунктирная линия по диагонали представляет случайную классификацию выбросов и инлиеров.

import math

from sklearn.metrics import RocCurveDisplay

cols = 2
pos_label = 0  # mean 0 belongs to positive class
datasets_names = y_true.keys()
rows = math.ceil(len(datasets_names) / cols)

fig, axs = plt.subplots(nrows=rows, ncols=cols, squeeze=False, figsize=(10, rows * 4))

for ax, dataset_name in zip(axs.ravel(), datasets_names):
    for model_idx, model_name in enumerate(model_names):
        display = RocCurveDisplay.from_predictions(
            y_true[dataset_name],
            y_pred[model_name][dataset_name],
            pos_label=pos_label,
            name=model_name,
            ax=ax,
            plot_chance_level=(model_idx == len(model_names) - 1),
            chance_level_kw={"linestyle": ":"},
        )
    ax.set_title(dataset_name)
_ = plt.tight_layout(pad=2.0)  # spacing between subplots

Мы наблюдаем, что после настройки числа соседей LOF и IForest показывают схожую производительность в терминах ROC AUC для наборов данных forestcover и cardiotocography. Значение для IForest немного лучше для набора данных SA, а LOF значительно лучше работает на наборе данных Ames Housing по сравнению с IForest.

Однако следует помнить, что Isolation Forest, как правило, обучается гораздо быстрее, чем LOF, на наборах данных с большим количеством выборок. LOF требует вычисления парных расстояний для нахождения ближайших соседей, что имеет квадратичную сложность относительно количества наблюдений. Это может сделать этот метод неприемлемым для больших наборов данных.

Изучение влияния параметров

В этом разделе мы исследуем влияние гиперпараметра n_neighbors и выбора масштабирования числовых переменных на модель LOF. Здесь мы используем набор данных Типы лесов в качестве бинарно закодированных категорий, которые вводят естественный масштаб евклидовых расстояний от 0 до 1. Затем мы хотим метод масштабирования, чтобы избежать предоставления преимущества небинарным признакам и достаточно устойчивый к выбросам, чтобы задача их поиска не стала слишком сложной.

X = X_forestcover
y = y_true["forestcover"]

n_samples = X.shape[0]
n_neighbors_list = (n_samples * np.array([0.2, 0.02, 0.01, 0.001])).astype(np.int32)
model = make_pipeline(RobustScaler(), LocalOutlierFactor())

linestyles = ["solid", "dashed", "dashdot", ":", (5, (10, 3))]

fig, ax = plt.subplots()
for model_idx, (linestyle, n_neighbors) in enumerate(zip(linestyles, n_neighbors_list)):
    model.set_params(localoutlierfactor__n_neighbors=n_neighbors)
    model.fit(X)
    y_pred = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_pred,
        pos_label=pos_label,
        name=f"n_neighbors = {n_neighbors}",
        ax=ax,
        plot_chance_level=(model_idx == len(n_neighbors_list) - 1),
        chance_level_kw={"linestyle": (0, (1, 10))},
        linestyle=linestyle,
        linewidth=2,
    )
_ = ax.set_title("RobustScaler with varying n_neighbors\non forestcover dataset")

Мы наблюдаем, что количество соседей оказывает большое влияние на производительность модели. Если у вас есть (по крайней мере некоторые) истинные метки, важно настроить n_neighbors соответственно. Удобный способ сделать это — исследовать значения для n_neighbors порядка величины ожидаемой аномалии.

from sklearn.preprocessing import MinMaxScaler, SplineTransformer, StandardScaler

preprocessor_list = [
    None,
    RobustScaler(),
    StandardScaler(),
    MinMaxScaler(),
    SplineTransformer(),
]
expected_anomaly_fraction = 0.02
lof = LocalOutlierFactor(n_neighbors=int(n_samples * expected_anomaly_fraction))

fig, ax = plt.subplots()
for model_idx, (linestyle, preprocessor) in enumerate(
    zip(linestyles, preprocessor_list)
):
    model = make_pipeline(preprocessor, lof)
    model.fit(X)
    y_pred = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_pred,
        pos_label=pos_label,
        name=str(preprocessor).split("(")[0],
        ax=ax,
        plot_chance_level=(model_idx == len(preprocessor_list) - 1),
        chance_level_kw={"linestyle": (0, (1, 10))},
        linestyle=linestyle,
        linewidth=2,
    )
_ = ax.set_title("Fixed n_neighbors with varying preprocessing\non forestcover dataset")

С одной стороны, RobustScaler масштабирует каждый признак независимо, используя межквартильный размах (IQR) по умолчанию, который представляет собой размах между 25-м и 75-м процентилями данных. Он центрирует данные, вычитая медиану, а затем масштабирует их, деля на IQR. IQR устойчив к выбросам: медиана и межквартильный размах менее чувствительны к экстремальным значениям, чем размах, среднее значение и стандартное отклонение. Кроме того, RobustScaler не сжимает маргинальные выбросы, в отличие от StandardScaler.

С другой стороны, MinMaxScaler масштабирует каждый признак индивидуально так, что его диапазон отображается в диапазоне от нуля до единицы. Если в данных есть выбросы, они могут сместить его к минимальному или максимальному значениям, что приводит к совершенно другой распределению данных с большими маргинальными выбросами: все значения, не являющиеся выбросами, в результате могут быть сжаты почти вместе.

Мы также оценили отсутствие предобработки (передав None в конвейер), StandardScaler и SplineTransformer. Для получения более подробной информации обратитесь к соответствующей документации.

Обратите внимание, что оптимальная предобработка зависит от набора данных, как показано ниже:

X = X_cardiotocography
y = y_true["cardiotocography"]

n_samples, expected_anomaly_fraction = X.shape[0], 0.025
lof = LocalOutlierFactor(n_neighbors=int(n_samples * expected_anomaly_fraction))

fig, ax = plt.subplots()
for model_idx, (linestyle, preprocessor) in enumerate(
    zip(linestyles, preprocessor_list)
):
    model = make_pipeline(preprocessor, lof)
    model.fit(X)
    y_pred = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_pred,
        pos_label=pos_label,
        name=str(preprocessor).split("(")[0],
        ax=ax,
        plot_chance_level=(model_idx == len(preprocessor_list) - 1),
        chance_level_kw={"linestyle": (0, (1, 10))},
        linestyle=linestyle,
        linewidth=2,
    )
ax.set_title(
    "Fixed n_neighbors with varying preprocessing\non cardiotocography dataset"
)
plt.show()

Общее время выполнения скрипта: (0 минут 51.064 секунды)

Download Jupyter notebook: plot_outlier_detection_bench.ipynb

Download Python source code: plot_outlier_detection_bench.py

Download zipped: plot_outlier_detection_bench.zip

Связанные примеры

Обнаружение выбросов с помощью Local Outlier Factor (LOF)

Сравнение алгоритмов обнаружения аномалий для выявления выбросов на наборах данных игрушечного типа

Сравнение эффекта различных масштабируемых преобразований данных с выбросами

Обнаружение новых выбросов с помощью Local Outlier Factor (LOF)