# CADE — интересный способ поиска аномалий в многомерных данных

### Введение

Один из способов искать аномалии в наборе данных — использовать соответствующую данным плотность случайной величины как оценку «аномальности» точек. Интуиция проста — чем меньше значение плотности распределения в точке, тем более эта точка «аномальна». И наоборот, области с высокими значениями плотности распределения соответствуют наборам более «нормальных» точек. То есть зная плотность распределения, мы, в одной из возможных постановок задачи, можем находить и аномалии. И именно в обладании (а реальной жизни в _приближении_) плотности вероятности и возникает трудность, особенно когда мы работаем с данными большой размерности, и стандартные способы приближения плотности становятся вычислительно дорогими. На помощь приходят более «хитрые» способы приблизить плотность вероятности. Существует целый класс таких методов, и один из них произвёл на меня особое впечатление оригинальностью своей мысли. Он показался мне достаточно интересным для того, чтобы более популярно изложить его на Хабре. Этот метод называется **CADE (Classifier Adjusted Density Estimation)**, и впервые он был изложен еще в 2001 году в статье \[1], хотя его реальный потенциал был исследован позже.

CADE — это метод приближения плотности распределения, который хорошо справляется с большими размерностями и неинформативными признаками. Что важно, для корректной работы он требует предположения о независимости признаков. Однако по утверждениям [этой статьи](https://seppe.net/aa/papers/CADE.pdf) \[2], метод ранжирует точки по аномальности на том же уровне, а при некоторых условиях даже лучше, чем _LOC (Local Outlier Factor)_ — один из популярных density‑based методов для поиска аномалий. Суть работы CADE заключается в использовании любого любимого вами алгоритма классификации, обученного различать известные данные от искусственно сгенерированных точек. И как это нам поможет приблизить плотность распределения? Итак, давайте рассмотрим.

### Как работает CADE

Если вкратце, то CADE предлагает следующее:

* Смешать исходные данные с данными из произвольного (но известного заранее) распределения;
* Построить классификатор, способный различать данные из двух распределений, то есть решить задачу двухклассовой классификации;
* Использовать предсказания классификатора в аналитически выведенной формуле, чтобы восстановить плотность распределения изначальных данных.

Как это работает? Давайте разбираться.

Я буду использовать нотации из вышеупомянутой [статьи](https://seppe.net/aa/papers/CADE.pdf) \[2], на которую я и опирался, чтобы вам было проще её перечитать, если у вас возникнет такое желание.

Пусть _T_ — это имеющиеся данные (с _небольшой_ долей аномалий). _P(X|T)_ — их плотность распределения, которую мы хотим приблизить. Первое, что нам надо сделать — это выбрать то самое «произвольное» распределение, из которого мы можем сгенерировать данные. Обозначим его как _A,_ а соответствующую плотность распределения — _P(X|A)_. Авторы статьи утверждают, что в качестве такого распределения хорошо подходит равномерное или нормальное распределения, хотя ограничений на выбор распределения нет. Желательно только, чтобы оно полностью покрывало данные, то есть была ненулевая вероятность выбрать точку, близкую к изначальным данным. Что касается размера выборки, то авторы использовали размер, равный размеру выборки _T_, однако я поделюсь своими соображениями на этот счёт позже. Точки из _A_ иногда называют _искусственными аномалиями_.

Следующий шаг — выбрать классификатор. В качестве него мы можем использовать всё, что способно рассчитать вероятность принадлежность точки к изначальному классу _T_, а именно _P(C = T|X)_. Классификатором может быть Random Forest Classifier, KNN, Decision Tree Classifier, подходящая по архитектуре нейронная сеть и так далее.

Естественно, для каждой точки сумма вероятностей быть в классе _T_ и классе _A_ равна 1. Другими словами:

![P(C=T|X)=1-P(C=A|X)](https://habrastorage.org/getpro/habr/upload_files/246/21b/5bd/24621b5bdd143f2b11b20a8326baa853.svg)

Теперь приступим к выводу формулы, которая позволит нам использовать предсказание классификатора для приближения изначальной плотности распределения _P(X|T)_. Для этого мы применим теорему Байеса к _P(C=T|X)_:

![P(C=T|X)=\frac{P(X|T)P(C=T)}{P(X)}=\frac{P(X|T)P(C=T)}{P(X|T)P(C=T)+P(X|A)P(C=A)}](https://habrastorage.org/getpro/habr/upload_files/ea5/af4/092/ea5af40928ed68d419a52dbd7280953c.svg)

Далее, решив это уравнение относительно _P(X|T)_, получим следующую формулу:

![P(X|T)=\frac{P(C=A)}{P(C=T)}P(X|A)\frac{P(C=T|X)}{1-P(C=T|X)}](https://habrastorage.org/getpro/habr/upload_files/510/ccb/f73/510ccbf73eb4cdcd2b640173d8bee2f9.svg)

Это и есть приближение искомой плотности распределения по CADE! Оно состоит из трёх основных множителей.

1. Это отношение, которое может быть приближено как отношение количества элементов в сгенерированном и изначальном множестве данных:

![\frac{P(C=A)}{P(C=T)}](https://habrastorage.org/getpro/habr/upload_files/c5a/e83/360/c5ae83360114140f99662bec90435487.svg)

2. Известная плотность распределения у сгенерированных данных:

![P(X|A)](https://habrastorage.org/getpro/habr/upload_files/b8e/734/314/b8e7343149a8d11fcdd1db35a03bc76d.svg)

3. Отношение, вычисляемое по предсказаниям классификатора:

![\frac{P(C=T|X)}{1-P(C=T|X)}](https://habrastorage.org/getpro/habr/upload_files/667/3d7/a47/6673d7a4704895a15d3326e7ebcbae06.svg)

Если _P(X|A)_ — константа, то есть _A_ было выбрано из равномерного распределения, то ранжирование точек по значениям плотности распределения будет совпадать с ранжированием по значениям классификатора (умножение на положительную константу не повлияет на порядок).

Если же мы угадали и выбрали _P(X|A)_ близким к _P(X|T),_ то классификатор не сможет различить распределения и для всех точек предскажет значение 0.5. Тогда ранжирование логичным образом будет совпадать с ранжированием по _P(X|A)._

Покончим с формулами! Давайте посмотрим на небольшую анимацию о процессе работы CADE без вычислений и финальной формулы.

<figure><img src="https://habrastorage.org/getpro/habr/upload_files/1dc/8f2/c88/1dc8f2c88cbd220cb2101ea5856abb9b.gif" alt="Только области с большим количеством реальных данных определяются как области с &#x22;нормальными&#x22; значениями" height="975" width="1776"><figcaption><p>Только области с большим количеством реальных данных определяются как области с "нормальными" значениями</p></figcaption></figure>

Как я упомянул выше, исследователи из \[2] показали, что метод действительно хорошо работает с данными высокой размерности и при наличии «шумных» атрибутов. А также использование множителя _P(X|A)_, что является отличительной особенностью CADE, улучшает работу метода при наличии коррелирующих признаков. Если у вас будет желание, я рекомендую ознакомиться с результатами из этого исследования. В нем вы также увидите значения метрики ROC AUC для разных наборов данных и классификаторов.

### Работаем с CADE на практике

На момент написания статьи, я не обнаружил реализации CADE на Python _(не путать с другим CADE — Compass‑aligned Distributional Embeddings, реализация для которого_ [_есть в PyPI_](https://pypi.org/project/cade/)_)._ Поэтому в этой статье, а также на Github ([cade‑outliers](https://github.com/AntonSarr/cade-outliers)), я оставлю пример кода, который реализует поиск аномалий с помощью этого метода.

Предложения по улучшению кода приветствуются!

Класс `CADEOutliers` отвечает за то, чтобы приближать плотность распределения исходных данных и возвращать точки, ранжированные по значениям плотности.

```python
import numpy as np


class CADEOutliers():
    """Класс для обнаружения выбросов с использованием CADE (Classifier Adjusted Density Estimation)"""
    
    _supported_A_dist = ["uniform"]
    
    def __init__(self, classifier, A_dist, A_size):
        """Создать объект обнаружения выбросов CADE.
        
        Параметры
        ----------
        classifier : model
            Модель классификатора sklearn, поддерживающая методы fit и predict_proba.
        
        A_dist : string
            Распределение для генерации искусственных аномалий.
        
        A_size : float | int
            Если float, то количество искусственных аномалий равно |X| * A_size, где X - анализируемые данные. 
            Если int, то количество искусственных аномалий равно A_size.
        """
        # Установка классификатора
        self.classifier = classifier
        
        # Установка типа распределения
        if A_dist not in self._supported_A_dist:
            raise ValueError(f"{A_dist} распределение пока не поддерживается для генерации искусственных аномалий. Выберите из следующего списка: {self._supported_A_dist}")
        self.A_dist = A_dist
        
        # Установка количества искусственных аномалий для генерации
        if isinstance(A_size, int) and A_size >= 1:
            self.sample_size = A_size
        elif isinstance(A_size, float) and A_size >= 0:
            self.sample_size = None
        else:
            raise ValueError("A_size должен быть либо int[1, inf) либо float[0, inf)")
        self.A_size = A_size
        
    def _generate_A(self, dist, attrs):
        """Генерация искусственных аномалий.
        
        Параметры
        ----------
        dist : string
            Распределение для генерации искусственных аномалий.
        attrs : dict
            Параметры для генерации искусственных аномалий. 
            Если dist равен "uniform", attrs = {'low': array, 'high': array, 'size': (sample_size, n_dim)} с минимальными и максимальными порогами для каждого измерения.
        """
        if dist == 'uniform':
            generator = np.random.uniform
        else:
            raise ValueError(f"{dist} распределение пока не поддерживается для генерации искусственных аномалий. Выберите из следующего списка: {self._supported_A_dist}")
        
        A = generator(**attrs)
        
        return A
    
    def fit(self, X):
        """Оценка распределения X.
        
        Параметры
        ----------
        X : np.ndarray
            Массив с образцами
        """
        if self.sample_size is None:
            sample_size = int(len(X) * self.A_size) + 1
        else:
            sample_size = int(self.sample_size)
        
        if self.A_dist == 'uniform':
            n_dim = X.shape[1]
            attrs = {
                'low': X.min(axis=0),
                'high': X.max(axis=0),
                'size': (sample_size, n_dim)
            }
            # вычисление плотности вероятности для равномерного распределения
            P_A = 1
            for s in X.max(axis=0) - X.min(axis=0):
                if s == 0:
                    s = 1
                P_A /= s
            self.P_A = lambda x: np.array([P_A] * x.shape[0])
        else:
            attrs = None
            
        A = self._generate_A(self.A_dist, attrs)
        
        combined_data = np.vstack([X, A])
        target = np.hstack([np.zeros(X.shape[0]), np.ones(A.shape[0])])
        
        self.classifier.fit(combined_data, target)
        self.A_X_ratio = A.shape[0] / X.shape[0]
        
        
    def get_ranking(self, X):
        """Возвращает ранжирование для каждого образца в X.
        Более высокие значения указывают на большую уверенность в том, что точка является выбросом.
        
        Параметры
        ----------
        X : np.ndarray
            Массив с образцами
        """
        predictions = self.classifier.predict_proba(X)[:, 0]
        
        anomaly_score = self.P_A(X) * self.A_X_ratio * (predictions / (1 - predictions + 1e-5))
        
        return anomaly_score
```

Генерирование данных и применение к ним CADE:

```python
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

from sklearn.metrics import roc_auc_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

from my_module import CADEOutliers  # Убедитесь, что CADEOutliers импортируется из правильного модуля


# Генерация данных
# X состоит из смеси 3 нормальных распределений
# 3-е распределение - это небольшая часть аномалий
X = np.vstack([2 * np.random.randn(5000, 2), 7 + np.random.randn(3000, 2), 25 + np.random.randn(100, 2)])
y = np.hstack([np.zeros(8000), np.ones(100)])

print('Размерность X: {}'.format(X.shape))

# Разделение данных на обучающую и тестовую части
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)

# Создание объекта CADE для обнаружения выбросов с использованием равномерного распределения искусственных аномалий размером 100% от заданного набора данных
cade = CADEOutliers(
    classifier=RandomForestClassifier(max_depth=3), 
    A_dist='uniform', 
    A_size=1.
)

# Обучение оценки плотности вероятности X с использованием X_train
cade.fit(X_train)

# Получение ранжирования для X_test. Меньшие значения указывают на более высокое сходство с аномалиями
ranking_by_dens = cade.get_ranking(X_test)

print('ROC AUC ранжирования аномалий с использованием CADE на тестовых данных: {}'.format(roc_auc_score(y_test, -ranking_by_dens)))

fig, ax = plt.subplots(1, 1, figsize=(7, 7))
sns.scatterplot(x=X_test[:, 0], y=X_test[:, 1], hue=ranking_by_dens, ax=ax, palette='magma')
plt.title('Точки с ранжированием, полученным с помощью CADE (низкие значения указывают на аномалии)')
plt.show()
```

Несколько наблюдений, которые могут оказаться полезными для вас:

1. Переобученный классификатор более склонен считать любые увиденные в процессе обучения точки «нормальными» (даже те, что кажутся аномалиями), а любые области без исходных данных при обучении — аномальными (даже те, что лежат рядом с увиденными точками). Другими словами, переобучение классификатора ведет к более детальному описанию плотности вероятности и большей чувствительностью к обнаружению аномалий.
2. Бóльшая плотность заполнения пространства искусственными аномалиями также ведет к более детальному описанию плотности распределения, т. е. к тому же эффекту, что и переобученный классификатор.
3. Из статьи \[2] — в качестве классификатора неплохо подходит Random Forest Classifier и KNN, а в качестве распределения для искусственных аномалий — равномерное.

***

Я буду рад, если CADE показался интересным и вам. Особенно, если вам захочется применить его на практике и поделиться своими наблюдениями о его работе. Делитесь своим мнением в комментариях!

### Источники

\[1] T. Hastie, R. Tibshirani, and J. H. Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. New York: Springer-Verlag, 2001.

\[2] [https://seppe.net/aa/papers/CADE.pdf](https://seppe.net/aa/papers/CADE.pdf)
