Обнаружение аномалий во временных рядах (на примере цен на отели)

Необходимость обнаружения необычных наблюдений (выбросов, или аномалий) во временных рядах часто возникает в таких ситуациях, как мониторинг состояния оборудования, отслеживание неожиданных колебаний на рынке ценных бумаг, учет показателей состояния здоровья пациентов и т.д. Существует большое количество методов и инструментов, позволяющих решать подобные задачи (см., например, Шкодырев и др. 2017). В этом сообщении мы рассмотрим один из таких инструментов, реализованных в R - пакет anomalize.

Данные, рассматриваемые в примерах

Для описания основ работы с пакетом anomalize мы воспользуемся набором данных по ценам на отели, предоставленным компанией Expedia в рамках одного из Kaggle-соревнований. Нам будет достаточно лишь небольшой части из исходной обучающей выборки - данные по суточной стоимости номеров в трех отелях за период, охватывающий примерно 8 месяцев. Этот небольшой набор данных можно найти в Github-репозитории ranalytics/intro_to_anomalize, где хранится также весь код, представленный в этом сообщении.

Структура данных очень проста. Имеем таблицу с тремя столбцами - prop_id (идентификатор отеля), date_time (время, когда пользователь увидел цену на соответствующий гостиничный номер в результате онлайн-поиска) и price_usd  (цена в долларах США; рис. 1):

# ---------------------- Инсталляция необходимых пакетов ---------------------

required_packages <- c("tidyverse", "anomalize", "tibbletime")
existing_packages <- installed.packages()[,"Package"]
to_install <- required_packages[!required_packages %in% existing_packages]
if (length(to_install) > 0) {install.packages(required_packages)}

require(tidyverse)
require(anomalize)
require(tibbletime)

# ------------------- Загрузка данных знакомство с ними ----------------------

dat <- read_csv("data/expedia_hotel_prices.csv")

dat %>% glimpse()

## Observations: 1,647
## Variables: 3
## $ prop_id   <dbl> 83045, 13252, 73738, 73738, 13252...
## $ date_time <dttm> 2012-12-31 08:59:22, 2013-03-20 17:50:44...
## $ price_usd <dbl> 219.00, 160.00, 189.00, 189.00, 186.00....

dat %>% 
    ggplot(., aes(date_time, price_usd)) +
    geom_line() +
    facet_wrap(~prop_id)

Рис. 1

Автоматическое обнаружение аномалий

Для начала попробуем обнаружить необычные наблюдения в одном из трех имеющихся  временных рядов (например, prop_id = 13252). В пакете anomalize для этого к таблице с данными (класса tibble или tibbletime) необходимо последовательно применить команды с участием следующих функций:

• time_decompose() - выполняет разложение временного ряда на отдельные составляющие (сезонную, тренд и случайный шум, или остатки);
• anomalize() - применяет к остаткам один из двух реализованных в пакете методов выявления аномалий (см. ниже);
• time_recompose() - вычисляет верхнюю и нижнюю границы диапазона, в который входят "нормальные" наблюдения.

Получаемый в итоге объект имеет следующую структуру:

result_13252 <- dat %>% 
    filter(prop_id == 13252) %>% 
    time_decompose(price_usd, merge = TRUE) %>%
    anomalize(remainder) %>%
    time_recompose()
## frequency = 2 hours
## trend = 61 hours

result_13252 %>% glimpse()

## Observations: 521
## Variables: 12
## $ prop_id       <dbl> 13252, 13252, 13252, 13252, 13252, 13252, 13...
## $ date_time     <dttm> 2012-11-01 10:16:16, 2012-11-02 17:51:01, 2...
## $ price_usd     <dbl> 184.00, 203.00, 203.00, 186.00, 169.00, 177....
## $ observed      <dbl> 184.00, 203.00, 203.00, 186.00, 169.00, 177....
## $ season        <dbl> 3.345308, -3.345308, 3.345308, -3.345308, 3....
## $ trend         <dbl> 192.4586, 193.1835, 193.9083, 194.6331, 195....
## $ remainder     <dbl> -11.803930, 13.161848, 5.746393, -5.287829, ...
## $ remainder_l1  <dbl> -133.7496, -133.7496, -133.7496, -133.7496, ...
## $ remainder_l2  <dbl> 141.1883, 141.1883, 141.1883, 141.1883, 141....
## $ anomaly       <chr> "No", "No", "No", "No", "No", "No", "No", "N...
## $ recomposed_l1 <dbl> 62.05432, 56.08854, 63.50400, 57.53822, 64.9...
## $ recomposed_l2 <dbl> 336.9923, 331.0265, 338.4420, 332.4762, 339....

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

• Результатом применения time_decompose() к исходным данным стал расчет четырех столбцов: observed (исходные наблюдения зависимой переменной), season (сезонная составляющая временного ряда), trend (тренд) и remainder (остатки). Параметр merge = TRUE привел к тому, что эти новые столбцы были добавлены  к исходным данным.
• Функция anomalize() была применена к столбцу remainder и рассчитала три новых столбца - remainder_l1 и remainder_l2 (нижняя и верхняя границы, за пределами которых наблюдения считаются аномалиями), а также индикаторную переменную anomaly, которая принимает значение "No", если наблюдение не является выбросом, и "Yes" если да.
• Функция time_recompose() восстанавила исходный временной ряд с использованием значений из season, trend, remainder_l1 и remainder_l2 и попутно рассчитала два новых столбца - recomposed_l1 и recomposed_l2, которые соответствуют верхней и нижней границам диапазона "нормальных" наблюдений. 

Суть всех этих вычислений станет понятнее, если мы представим результаты графически. Функция plot_anomaly_decomposition() изображает компоненты временного ряда, рассчитываемые функцией time_decompose(). На получаемом графике аномалии показаны в виде обведенных кружком красных точек (рис. 2):

result_13252 %>% plot_anomaly_decomposition()

Рис. 2

Еще одна функция - plot_anomalies() - позволяет отдельно изобразить исходный временной ряд и обнаруженные в нем необычные наблюдения:

result_13252 %>% plot_anomalies()

Рис. 3

Поскольку plot_anomalies() возвращает графический объект класса ggplot, мы при желании можем легко добавить к нему другие ggplot-элементы обычным в таких случаях способом. Добавим, например, линию, последовательно соединяющую все точки:

result_13252 %>% plot_anomalies() %>% geom_line()

Рис. 4

Ручная настройка параметров

При создании объекта result_13252 все вычисления были выполнены с использованием настроек, заданных по умолчанию. Часто для получения приемлемого результата этого будет достаточно. Однако, в зависимости от свойств анализируемого временного ряда, иногда потребуется более тонкая настройка соответствующих параметров.

Функция time_decompose()

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

• frequency - определяет "частоту" сезонной составляющей, т.е. число наблюдений, входящих в один сезонный цикл (например, для дневных данных этот аргумент по умолчанию будет равен "1 week", т.е. одной неделе);
• trend - определяет число наблюдений в отдельных отрезках временного ряда, используемых для расчета тренда (например, для дневных данных по умолчанию trend = "3 months", т.е. трем месяцам);
• method - метод, используемый для разложения ряда на компоненты. Этот аргумент принимает два возможных значения - "stl" (по умолчанию) и "twitter". STL ("Seasonal and Trend decomposition using Loess") - это классический метод разложения временных рядов, в основе которого лежит сглаживание данных по методу нелинейной локальной регрессии Loess. Метод "twitter" - более новый и, как следует из его названия, был предложен исследователями из компании Twitter (Vallis et al. 2014). Согласно этому методу, тренд оценивается путем вычисления медиан в пределах отдельных промежутков, на которые разбивает исходный ряд (см. выше - их длина определяется параметром trend).

По умолчанию аргументы frequency и trend принимают значения "auto" (автоматический выбор). В этом случае для выбора конкретных подходящих ситуации значений frequency и trend функция time_decompose() обращается к другой функции - get_time_scale_template(), которая возвращает "справочник" типичных значений:

get_time_scale_template()

##  A tibble: 8 x 3
##  time_scale frequency trend   
##  <chr>      <chr>     <chr>   
##  second     1 hour    12 hours
##  minute     1 day     14 days 
##  hour       1 day     1 month 
##  day        1 week    3 months
##  week       1 quarter 1 year  
##  month      1 year    5 years 
##  quarter    1 year    10 years
##  year       5 years   30 years

Видим, например, что в случае с данными, где расстояние между отдельными наблюдениями выражается в секундах, параметры frequency и trend будут выражаться в часах и по умолчанию составят "1 hour" (1 час) и "12 hours" (12 часов) соответственно. Следует отметить однако, что в зависимости от расстояния между наблюдениями и длины самого временного ряда, конкретные выбранные программой значения frequency и trend необязательно окажутся равными приведенным выше значениям из "справочника". Именно это произошло выше при создании объекта result_13252 (frequency = 2 hours, trend = 61 hours). Помимо строковых значений ("1 hour", "2 days", "2 weeks" и т.д.), обоим параметрам можно присваивать и числовые значения (см. ниже) - часто это оказывается более удобным подходом. Отметим также, что "шаблонные" значения параметров frequency и trend можно изменить "глобально" с помощью функции set_time_scale_template() (пример не приводится).

Для отмены автоматического выбора значений frequency и trend достаточно указать желаемые значения при вызове функции time_decompose(). В приведенном ниже примере мы присвоим этим параметрам значения 2 и 6 (часов) соответственно (рис. 5):

dat %>% 
    filter(prop_id == 13252) %>% 
    time_decompose(price_usd, 
                   frequency = 2,
                   trend = 6) %>%
    anomalize(remainder) %>%
    time_recompose() %>% 
    plot_anomaly_decomposition()

Рис. 5

Как видно из рис. 5, уменьшение параметров frequency и trend привело к классификации большего количества наблюдений как аномалий. Это обычный результат в таких случаях, обусловленный переобучением Loess-модели тренда (обратите внимание на то, каким негладким стал тренд на рис. 5 по сравнению с рис. 2). В отличие от Loess, второй метод разложения временных рядов - "twitter" - более устойчив к локальным колебаниям и поэтому при тех же значениях frequency и trend обычно классифицирует меньшее количество наблюдений как аномалии. Это легко проверить (сравните рис. 5 и рис. 6):

dat %>% 
    filter(prop_id == 13252) %>% 
    time_decompose(price_usd, 
                   frequency = 2,
                   trend = 6, 
                   method = "twitter") %>%
    anomalize(remainder) %>%
    time_recompose() %>% 
    plot_anomaly_decomposition()

Рис. 6

Функция anomalize()

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

• method - принимает два возможных значения в соответствии с названием метода обнаружения аномалий: "iqr" (по умолчанию) и "gesd". Метод IQR (от "interquartile range", интерквартильный размах) работает быстрее, но не так точен, как GESD ("Generalized Extreme Studentized Deviate test").
• alpha - контролирует ширину диапазона "нормальных" значений (по умолчанию равен 0.05). Чем меньше этот параметр, тем шире диапазон значений, рассматриваемых как "нормальные", и тем меньше наблюдений будут классифицированы как аномальные.
• max_anoms - максимальная доля наблюдений, которые позволено классифицировать как аномалии (0.2 по умолчанию, т.е. 20%).

Рассмотрим влияние параметра alpha, снизив его значение с принятого по умолчанию 0.05 до 0.025. При этом оставим параметры frequency и trend функции time_decompose() такими же, как на рис. 5:

dat %>% 
    filter(prop_id == 13252) %>% 
    time_decompose(price_usd, 
                   frequency = 2,
                   trend = 6, 
                   method = "twitter") %>%
    anomalize(remainder, alpha = 0.025) %>%
    time_recompose() %>% 
    plot_anomaly_decomposition()

Рис. 7

Как видно из рис. 7, несмотря на все еще имеющее место переобучение Loess-модели тренда, к классу аномальных было отнесено намного меньшее количество наблюдений. Аналогичного эффекта можно было бы добиться также путем снижения параметра max_anoms функции anomalize() (пример не приводится).

Одновременный анализ нескольких рядов

Отличительной особенностью пакета anomalize является возможность параллельно анализировать несколько временных рядов на наличие аномалий. Для этого достаточно в таблице с данными иметь группирующую переменную, которая идентифицирует отдельные ряды. В нашем примере с ценами на гостиничные номера такой одновременный анализ выглядел бы следующим образом (параметрам всех функций здесь присвоены их автоматически выбранные значения; обратите внимание на аргумент time_recompose = TRUE функции plot_anomalies() - он позволяет включить изображение диапазона, в который входят "нормальные" наблюдения):

dat %>% 
    time_decompose(price_usd) %>%
    anomalize(remainder) %>%
    time_recompose() %>% 
    plot_anomalies(time_recomposed = TRUE)

Рис. 8

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

Заключение

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

Пакет anomalize позволяет обнаруживать лишь единичные точки-выбросы. Другим распространенным типом аномалий является сдвиг среднего уровня временного ряда. В следующем сообщении я опишу один из эффективных методов для выявления подобных сдвигов.