15.12.2017 08:10
little-brother 10 6000 Источник
Когда наблюдаемых метрик становится слишком много, отслеживание всех графиков самостоятельно становится невозможным. Обычно в этом случае для менее значимых метрик используют проверки на достижение критичных значений. Но даже если значения подобраны хорошо, часть проблем остается незамеченной. Какие это проблемы и как их обнаруживать — под катом.



Disclaimer
Автор хотя и имеет математическое образование, никак не связан ни с Data Mining, ни со статистическим анализом. Данный материал является результатом исследования, проведенного с целью выяснить возможность написания модуля поиска аномалий (пусть даже слабого) для разрабатываемой системы мониторинга.

Что ищем в двух картинках


Примеры аномалий на графиках

Источник Anomaly.io

Конечно в реальности, не всегда все так просто: только на б), д) и е) явная аномалия.

Источник cyberleninka.ru

Текущее положение дел


Коммерческие продукты почти всегда представлены в виде сервиса, использующего как статистику, так и машинное обучение. Вот некоторые из них: AIMS, Anomaly.io (прекрасный блог с примерами), CoScale (возможность интеграции, напр. с Zabbix), DataDog, Grok, Metricly.com и Azure (от Microsoft). У Elastic есть модуль X-Pack на основе машинного обучения.

Open-source продукты, которые можно развернуть у себя:

  • Atlas от Netflix — in-memory база данных для анализа временных рядов. Поиск выполняется методом Хольта-Винтерса.
  • Banshee от Eleme. Для обнаружения используется всего один алгоритм, основанный на правиле трех сигм.
  • Стек KALE от Etsy, состоящий из двух продуктов Skyline, для обнаружения проблем в реальном времени, и Oculus для поиска взаимосвязей. Похоже заброшен.
  • Morgoth — фреймворк для обнаружения аномалий в Kapacitor (модуль уведомлений InfluxDB). Из коробки имеет: правило трех сигм, тест Колмогорова-Смирнова и отклонение Дженсена-Шеннона.
  • Prometheus — набирающая популярность система мониторинга, в которой можно самостоятельно реализовывать некоторые алгоритмы поиска, используя стандартный функционал.
  • RRDTool — база данных, используемая, например, в Cacti и Munin. Умеет делать прогноз методом Хольт-Винтерса. О том, как это можно использовать — в статье Мониторинг прогнозированием, оповещения о потенциальном сбое
  • ~2000 репозитариев на GitHub

На мой взгляд open-source по качеству поиска значительно уступает. Чтобы понять, как работает поиск аномалий и можно ли исправить ситуацию, придется немного окунуться в статистику. Математические детали упрощены и скрыты под спойлерами.

Модель и её компоненты


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

  • Тренд — отражает общее поведение ряда в плане возрастания или убывания значений.
  • Сезонность — периодичные колебания значений, связанные, например, с днем недели или месяца.
  • Случайное значение — то, что останется от ряда после исключения других компонент. Здесь и будем искать аномалии.

Можно включить дополнительные компоненты, например циклическую, как множитель тренда,
аномальную (Catastrophic event) или социальную (праздники). Если тренд или сезонность в данных не просматриваются, то соответсвующие компоненты из модели можно исключить.

В зависимости от того, как связаны между собой компоненты модели, определяют её тип. Так, если все компоненты складываются, чтобы получить наблюдаемый ряд, то говорят, что модель аддиктивна, если умножаются, то мультипликативна, если что-то умножается, а что-то сладывается, то смешаная. Обычно тип модели выбирается исследователем на основе предварительного анализа данных.

Декомпозиция


Выбрав тип модели и набор компонент можно приступать к декомпозиции временного ряда, т.е. его разложению на компоненты.


Источник Anomaly.io

Сперва выделяем тренд, сгладив исходные данные. Метод и степень сглаживания выбираются исследователем.

Методы сглаживания рядов: скользящая средняя, экспоненциальное сглаживание и регрессия
Самый простой способ сгладить временной ряд — использовать вместо исходных значений половину суммы соседних

$s_n = \frac {x_n + x_{n-1}}{2}$


Если использовать не одно, а несколько предшедствующих значений, т.е. среднее арифметическое k-соседних значений, то такое сглаживание называется простым скользящим средним с шириной окна k

$s_n = \frac {x_n + x_{n-1} + ... + x_{n-k}}{k}$


Если для каждого предыдушего значения использовать какой то свой коэффициент, определяющий степень влияния на текущий, то получим взвешенное скользящее среднeе.

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

$s_n = \alpha x_n + (1 - \alpha)s_{n-1}$


Где ? — коэффициент сглаживания, от 0 до 1. Как легко видеть чем ближе ? к 1, тем больше получаемый ряд будет похож на исходный.

Для определения линейного тренда можно взять методику расчета линейной регрессии $y_t=a+bx_t+\varepsilon_t$ методом наименьших квадратов: $\hat{b}=\frac {\overline{xy}}{\overline{x^2}}$, $\hat{a}=\bar{y}-b \bar {x}$, где $\bar {x}$ и $\bar {y}$ — средние арифметические $x$ и $y$.

Источник Википедия

Для определения сезонной составляющей из исходного ряда вычитаем тренд или делим на него, в зависимости от типа выбранной модели, и еще раз сглаживаем. Затем делим данные по длине сезона (периоду), обычно это неделя, и находим усредненный сезон. Если длина сезона не известна, то можно попытаться найти её:

Дискретным преобразованием Фурье или авто-корреляцией
Честно признаюсь, что не стал разбираться как работает преобразование Фурье. Кому интересно могут заглянуть в следующие статьи: Detect Seasonality using Fourier Transform in R и Простыми словами о преобразовании Фурье. Насколько я понял, исходный ряд/функция представляется в виде бесконечной суммы элементов и берется несколько первых значимых коэффициентов.


Для поиска авто-корреляцией просто сдвигаем функцию вправо и ищем такое положение, чтобы расстояние/площадь между исходной и сдвинутой функцией (выделено красным) было минимально. Очевидно для алгоритма должен быть задан шаг сдвига и максимальный предел, при достижении которого считаем, что поиск периода не удался.

Более подробно о моделях и декомпозиции можно узнать из статей «Extract Seasonal & Trend: using decomposition in R» и «How To Identify Patterns in Time Series Data».

Удалив из исходного ряда тренд и сезонный фактор, получаем случайную компоненту.

Типы аномалий


Если анализировать только случайную компоненту, то многие аномалии можно свести к одному из следующих случаев:

  • Выброс
    Выброс в данных временного ряда Нахождение выбросов в данных — классическая задача, для решения которой уже имеется хороший набор решений:

    Правило трех сигм, межквартильный размах и другие
    Идея подобных тестов — определить насколько далеко располагается отдельное значение от среднего. Если расстояние отличается от «обычного», то значение объявляется выбросом. Время события при этом игнорируется.

    Считаем, что на входе ряд чисел — $x_1$, $x_2$,… $x_n$, всего $n$ штук. $x_i$ — $i$-ое число.

    Стандартные тесты достаточно просты в реализации и требуют лишь вычисление среднего $\bar x = (x_1 + ... + x_n) / n$, стандартного отклонения $S=\sqrt{\frac{1}{n-1}\sum_{i=1}^n\left(x_i-\bar{x}\right)^2}$ и иногда медианы $M$ — среднее значение, если упорядочить все числа по возрастанию и взять то, которое по середине.

    Правило трех сигм
    Если $|x_i - \bar x| > 3 * S$, то $x_i$ считаем выбросом.

    Z-оценка и уточненный метод Iglewicz и Hoaglin
    $x_i$ — выброс, если $|x_i - \bar x|/S$ больше задаваемого порога, обычно равному 3. По сути переписанное правило трех сигм.

    Уточненный метод заключается в следующем: для каждого числа ряда вычисляем $|x_i - M|$ и для получившихся значений находим медиану, обозначаемую $MAD$.
    $x_i$ — выброс, если $0.6745 * |x_i - M|/MAD$ больше порога.

    Межквартильный размах
    Сортируем исходные числа по возрастанию, делим на двое и для каждой части находим медианные значения $Q_1$ и $Q_3$.
    $x_i$ — выброс, если $|x_i - M| > 1.5 * (Q_3 - Q_1)$.

    Тест Граббса
    Находим минимальное $x_{min}$ и максимальное $x_{max}$ значения и для них вычисляем $G_{min} = (\bar x - x_{min})/ S$ и $G_{max} = (x_{max} - \bar x) / S$. Затем выбираем уровень значимости ? (обычно один из 0.01, 0.05 или 0.1), заглядываем в таблицу критичных значений, выбираем значение для n и ?. Если $G_{min}$ или $G_{max}$ больше табличного значения, то считаем соответствующий элемент ряда выбросом.

    Обычно тесты требуют, чтобы исследовалось нормальное распределение, но зачастую это требование игнорируется.
  • Сдвиг
    Сдвиг в данных временного ряда Задача обнаружения сдвига в данных неплохо исследована, поскольку встречается в обработке сигналов. Для её решения можно воспользоваться Twitter Breakout Detection. Отмечу, что это оно работает значительно хуже, если анализируется исходный ряд с компонентами сезонности и тренда.
  • Изменение характера (распределения) значений
    Изменения распределения в данных временного ряда Как и для сдвига, можно использовать тот же пакет BreakoutDetection от Twitter, но с ним не все гладко.

  • Отклонение от «повседневного» (для данных с сезонностью)
    Для обнаружения данной аномалии необходимо сравнить текущий период и несколько предыдущих. Обычно используют

    Метод Хольта-Винтерса
    Метод относится к прогнозированию, поэтому его применение сводится к тому, чтобы сравнить прогнозируемое значение с действительным.

    Основная идея метода в том, что каждая из трех компонент экспоненциально сглаживается, используя отдельный коэффициент сглаживания, поэтому метод зачастую называется тройным экспоненциальным сглаживанием. Формулы рассчета для мультипликативного и аддиктивного сезонов есть в Википедии, а подробности о методе в статье на Хабре.

    Три параметра сглаживания должны выбираться так, чтобы получаемый ряд был «близок» к исходному. На практике такая задача решается перебором, хотя RRDTool требует явного задания этих значений.

    Недостаток метода: требует минимум три сезона данных.

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

    Также Twitter предлагает AnomalyDetection, который на тестовых данных показывает хорошие результаты. К сожалению, как и BreakoutDetection, он уже два года без обновлений.
  • Поведенческие
    Иногда аномалии могут иметь специфичный характер, который в большинстве случаев должен быть проигнорирован, например, одно значение у двух соседних элементов. В таких случаях требуются отдельные алгоритмы.
  • Совместные аномалии
    Отдельно стоит упомянуть о аномалиях, когда значения двух наблюдаемых метрик поотдельности в пределах нормы, но их совместное появление признак проблемы. В общем случае, нахождение таких аномалий может быть решено кластерным анализом, когда пары значений представляют координаты точек на плоскости и все множество событий делится на две/три группы (кластера) — «норма»/«шум» и «аномалия».


    Источник alexanderdyakonov.wordpress.com

    Более слабый метод состоит в том, чтобы отслеживать насколько метрики зависят друг от друга во времени и в случае, когда зависимость теряется, выдавать сообщение об аномалии. Для этого, вероятно, можно использовать один из методов.

    Линейный коэффициент корреляции Пирсона
    Пусть $x_1, x_2, ..., x_n$ и $y_1, y_2, ..., y_n$ два набора чисел и требуется выяснить имеется ли между ними линейная зависимость. Вычисляем для $x$ среднее $\bar x = (x_1 + ... + x_n) / n$ и стандартное отклонение $S_x=\sqrt{\frac{1}{n-1}\sum_{i=1}^n\left(x_i-\bar{x}\right)^2}$. Аналогично для $y$.

    Коэффициент корреляции Пирсона

    $r(x, y) = \frac {\sum_{i=1}^n(x_i - \bar x)(y_i - \bar y)} {S_x S_y}$


    Чем коэффициент по модулю ближе к 1, тем большая зависимость. Если коэффициент ближе к -1, то зависимость обратная, т.е. рост $x$ связан с убыванием $y$.

    Расстояние Дамерау — Левенштейна
    Пусть есть две строки ABC и ADEC. Чтобы получить из первой вторую, необходимо убрать B и добавить D и E. Если каждой операции удаления/добавления символа и перестановке XY в YX задать стоимость, то суммарная стоимость и будет расстояниеи Дамерау — Левенштейна.

    Для определения похожести графиков можно оттолкнуться от алгоритма, использованного в KALE
    Вначале исходный ряд значений, например, ряд вида [960, 350, 350, 432, 390, 76, 105, 715, 715], нормализуется: ищется максимум — ему будет соответствовать 25, и минимум — ему будет соответствовать 0; таким образом, данные пропорционально распределяются в пределе целых чисел от 0 до 25. В итоге мы получаем ряд вида [25, 8, 8, 10, 9, 0, 1, 18, 18]. Затем нормализованный ряд кодируется с помощью 5 слов: sdec (резко вниз), dec (вниз), s (ровно), inc (вверх), sinc (резко вверх). В итоге получается ряд вида [sdec, flat, inc, dec, sdec, inc, sinc, flat].

    Закодировав два временных ряда таким образом, можно узнать расстояние между ними. И если оно меньше какого то порога, то считать, что графики похожи и связь имеется.

Заключение


Разумеется, многие алгоритмы нахождения аномалий уже реализованы на языке R, предназначенном для статистической обработки данных, в виде пакетов: tsoutliers, strucchange, Twitter Anomaly Detection и других. Подробнее о R в статьях А вы уже применяете R в бизнесе? и Мой опыт введения в R. Казалось бы, подключай пакеты и используй. Однако есть проблема — задание параметров статистических проверок, которые в отличии от критических значений далеко не очевидны для большинства и не имеют универсальных значений. Выходом из данной ситуации может быть их подбор перебором (ресурсоёмко), с редким периодичным уточнением, независимо для каждой метрики. С другой стороны, большая часть аномалий, не связанных с сезонностью, хорошо определяется визуально, что наталкивает на мысль использовать нейронную сеть на отрендеренные графики.

Приложение


Ниже привожу собственные алгоритмы, которые работают сопоставимо с Twitter Breakout по результатам, и несколько быстрее по скорости при реализации на Java Script.

Алгоритм кусочно-линейной аппроксимации временного ряда
Аппроксимация
  1. Если ряд сильно зашумленный, то усредняем, напр. по 5 элементам.
  2. В результат включаем первую и последние точки ряда.
  3. Находим максимально удаленную точку ряда от текущей ломаной и добавляем ее в набор.
  4. Повторяем пока среднее отклонение от ломаной до исходного ряда не будет меньше среднего отклонения в исходном ряде или пока не достигнуто предельное число вершин ломаной (в этом случае, вероятно, апроксимация не удалась).

Алгоритм нахождения сдвига в данных
Обнаружение сдвига
  1. Апроксимируем исходный ряд ломаной
  2. Для каждого отрезка ломаной, кроме первого и последнего:
    • Находим его высоту $h$, как разность y-координат начала и конца. Если высота меньше игнорируемого интервала, то такой отрезок игнорируется
    • Оба соседних отрезка $L$ и $R$ делим на двое, каждый участок $L_1$, $L_2$, $R_1$, $R_2$ аппроксимируем своей прямой и находим среднее расстояние от прямой до ряда — $dL_1$, $dL_2$, $dR_1$, $dR_2$.
    • Если $|dL_1 - dL_2|$ и $|dR_1 - dR_2|$ значительно меньше чем $h$, то считаем, что сдвиг обнаружен

Алгоритм нахождения изменений в распределении
Обнаружение изменения распределения
  1. Исходный ряд разбивается на отрезки, длинной, зависящей от числа данных, но не менее трех.
  2. В каждом отрезке ищется минимум и максимум. Заменяя каждый отрезок его центром, образуется два ряда — минимумы и максимумы. Далее ряды обрабатываются отдельно.
  3. Ряд линейно аппроксимируется и для каждой его вершины, кроме первой и последней, выполняется сравнение данных исходного ряда, лежащих слева и справа до соседних вершин ломаной, тестом Колмогорова-Смирнова. Если разница обнаружена, то точка добавляется в результат.

Тест Колмогорова-Смирнова
Пусть $x_1, x_2, ..., x_{n_1}$ и $y_1, y_2, ..., y_{n_2}$ два набора чисел и требуется оценить существенность различий между ними.

Вначале значения обоих рядов разбиваются на несколько (около десятка) категорий. Далее для каждой категории вычисляется число $f_x$, вошедших в него значений из ряда $x$, и делится на длину ряда $n_1$. Аналогично для ряда $y$. Для каждой категории находим $|f_x - f_y|$ и затем общий максимум $d_{max}$ по всем категориям. Проверяемое значение критерия вычисляется по формуле $\lambda = d_{max} \frac {n_1 n_2} {n_1 + n_2}$.

Выбирается уровень значимости $\alpha$ (один из 0.01, 0.05, 0.1) и по нему определяется критичное значение по таблице. Если $\lambda$ больше критичного значения, то считается, что группы различаются существенно.

Комментарии (10)


  1. JekaMas
    15.12.2017 20:02
    #10569360

    Вариант от facebook, prophet не пробовали?


    1. little-brother Автор
      16.12.2017 12:20
      #10570168

      Prophet я видел, но не стал в обзор, поскольку как я понял, у них больше прогнозирование. В целом, если хочется учитывать праздничную составляющую из коробки, то да, можно им воспользоваться. Тут вопрос насколько это будет быстро.


  1. IvanKhozyainov
    16.12.2017 04:40
    #10569938

    Большая, интересная тема. Спасибо за статью. Да, в блоге anomaly.io неплохо написано про статистические методы. Кстати, если кому будет интересно посмотреть видео по похожей теме (хотя там больше про machine learning) — есть канал на ютуб www.youtube.com/user/numenta/videos, grokstream.com использует эти разработки в своем решении.

    У меня пару вопросов на понимание.

    1. Про совместные аномалии
    Хороший термин, раньше не встречал русского обозначения, но все-таки надо как-то обзывать эти многомерные аномалии, почему бы и не так. Он где-то использовался в русскоязычной литературе?
    Из возможных интуитивных подходов — использование сворачивания этой 2d-функции в одномерную и искать уже аномалии в значениях f(x,y). Например, как это делается в регрессионном анализе при поиске зависимостей.

    2. «Алгоритм нахождения сдвига в данных»
    Нахождение среднего тут неявно заложено в алгоритме аппроксимации. А как аппроксимировали? Своим методом, который «Алгоритм кусочно-линейной аппроксимации временного ряда»? Или метод наименьших квадратов?
    В этом смысле очень похоже на алгоритм, который используется в R-пакете твиттера, но подход с другой стороны. В математическом смысле суть похожа: «находим локальное среднее значение, сравниваем с соседними»
    А часто бывает когда результат апроксимации совпадает с «сырой» исходной ломаной?
    Мне правда это интересно с той точки зрения, насколько могут быть похожи или различаться решения одной и той же задачи. Может быть решение вообще только единственное верное в идеальном случае. Но это вопрос больше философский.

    P.S. Вашу статью, кстати, вполне можно как шпаргалку использовать для первичной навигации по теме :)
    С наступающим!


    1. little-brother Автор
      16.12.2017 13:13
      #10570212

      1. В русскоязычных материалах термин «совместные» вроде не встречал (да и в целом не так её и много просмотрел). Да, наверно можно и свертку сделать, у меня просто проблема в нехватке знаний :(

      2. Аппроксимировал своим методом. Аппроксимация мои методом близка к исходному ряду, если будет аппроксимироваться не достаточно сглаженный ряд. Возможно можно использовать какой то параметр гладкости внутри метода, чтобы так не промахиваться.
      Метод наверняка так себе, поскольку мне оказалось проще придумать своё, нежели использовать готовое. Через язык статей, тому кто в этой теме не постоянно, прорываться очень сложно.


      1. IvanKhozyainov
        17.12.2017 01:56
        #10570992

        1. Честно говоря, русскоязычной литературы по этой теме не так много находил — в основном все-таки англоязычная. Поэтому подумалось «а вдруг».
        2. Интересно, что в математике в этом смысле есть сходство с программированием. Готовые исследованные решения типа методов Эйлера как фреймворки, а самостоятельно придуманные методы — это что-то вроде своих самописных решений в программировании. Преимущества и недостатки очень похожи.


  1. A1essandro
    16.12.2017 12:20
    #10570170

    … что наталкивает на мысль использовать нейронную сеть на отрендеренные графики.

    По идее можно не рендерить, а просто использовать тот же персептрон. Входные данные просто нужно будет преобразовать: привести количество точек по x к количеству нейронов входного слоя, а значения y нормировать от [0, 1] либо [-1, 1].


    1. little-brother Автор
      16.12.2017 12:41
      #10570184

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


      1. A1essandro
        17.12.2017 12:56
        #10571222

        Еще про совместные аномалии и нейронные сети — обратите внимание на Нейронную Сеть Кохонена. Очень простая по структуре (два слоя). Легко реализуема. Обучение без учителя для этой сети тоже легче реализовать чем тот же самый алгоритм обратного распространения ошибки. Дает неплохие результаты по классификации. Правда иногда нужно поиграться с начальными весами.


  1. openbsod
    16.12.2017 12:20
    #10570172

    Благодарю вас, очень увлекательная статья. Если позволите, пара вопросов: на каком объеме данных (примерно) вы проводите вычисления, с какой частотой поступают метрики на анализ, к-примеру, «150k значений/sec, 300GB данных в час»? Спасибо.


    1. little-brother Автор
      16.12.2017 13:27
      #10570230

      Объема данных у меня никакого нет — несколько простых железок, опрашиваемых раз в две минуты и всё. Разрабатываемая система мониторинга уровня Cacti, а не Zabbix/Prometheus/Influx.

      Задачу выделения аномалий хорошо можно распараллелить, если использовать численные/статистические методы, но чтобы анализировать скажем 150К/с это уже надо делать балансировщики и сопутствующую архитектуру. Но даже в этом случае придется переписать все расчеты на низкоуровневые языки. R, Python и, используемый мной, Java Script — слишком медленные для такого.