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

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

Для разных моделей обнаружения выбросов требуются разные методы предобработки. При наличии категориальных переменных 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_score = estimator[-1].negative_outlier_factor_
    else:  # "IsolationForest"
        y_score = estimator.fit(X).decision_function(X)
    toc = perf_counter()
    print(f"Duration for {model_name}: {toc - tic:.2f} s")
    return y_score

В оставшейся части примера мы обрабатываем по одному набору данных в каждом разделе. После загрузки данных целевые значения преобразуются в два класса: 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_score = {"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_score[model_name]["KDDCup99 - SA"] = fit_predict(model, X)

Duration for LOF: 1.82 s
Duration for IForest: 0.29 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_score[model_name]["forestcover"] = fit_predict(model, X)

Duration for LOF: 1.81 s
Duration for IForest: 0.23 s

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

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

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_score[model_name]["ames_housing"] = fit_predict(model, X)

Duration for LOF: 0.46 s
Duration for IForest: 0.25 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_score[model_name]["cardiotocography"] = fit_predict(model, X)

Duration for LOF: 0.06 s
Duration for IForest: 0.15 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_score[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 на наборах данных Forest covertypes и кардиотокографии. На наборе данных SA оценка IForest немного выше, а на наборе Ames Housing LOF значительно превосходит IForest.

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

Абляционное исследование

В этом разделе мы изучаем влияние гиперпараметра n_neighbors и выбора способа масштабирования числовых переменных на модель LOF. Здесь мы используем набор данных Forest covertypes, поскольку категории, закодированные в двоичном виде, задают естественный масштаб евклидовых расстояний от 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_score = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_score,
        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))},
        curve_kwargs=dict(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_score = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_score,
        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))},
        curve_kwargs=dict(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_score = model[-1].negative_outlier_factor_
    display = RocCurveDisplay.from_predictions(
        y,
        y_score,
        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))},
        curve_kwargs=dict(linestyle=linestyle, linewidth=2),
    )
ax.set_title(
    "Fixed n_neighbors with varying preprocessing\non cardiotocography dataset"
)
plt.show()

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

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)