Прогнозирование временных рядов с помощью Prophet: ёмкость системы

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

Код для воспроизведения примеров

Все примеры, приведенные в сообщениях из этой серии, можно воспроизвести с помощью кода, который хранится в Github-репозитории ranalytics/intro_to_prophet. Описание используемых в примерах данных по стоимости биткоина представлено в первом сообщении.

Тренд с насыщением

Стандартным подходом для описания роста в системе с ограниченной емкостью является использование логистической функции следующего вида: \[g(t) = \frac{C}{1 + \exp(-k(t - m))}, \] где $C$ - верхний порог ("емкость системы"), $k$ - скорость роста, $t$ - время, а $m$ - параметр, позволяющий "сдвигать" функцию вдоль оси времени. На рис. 1 показаны примеры логистической функции с разными значениями параметров.


logistic_growth <- function(x) {
    C / (1 + exp(-k*(x - m)))
}

par(mfrow = c(1, 3))

# Пример 1:
C <- 10
k <- 0.1
m <- 0
curve(logistic_growth(x), from = -100, to = 100)
abline(v = 0, col = "red", lty = 2)
abline(h = 10, col = "blue", lty = 2)
title("C=10, k=0.1, m=0")

# Пример 2:
C <- 10
k <- 0.1
m <- 20
curve(logistic_growth(x), from = -100, to = 100)
abline(v = 20, col = "red", lty = 2)
abline(h = 10, col = "blue", lty = 2)
title("C=10, k=0.1, m=10")

# Пример 3:
C <- 8
k <- 0.05
m <- 0
curve(logistic_growth(x), from = -100, to = 100,
      ylim = c(0, 10))
abline(v = 0, col = "red", lty = 2)
abline(h = 8, col = "blue", lty = 2)
title("C=8, k=0.05, m=0")

Рис. 1

Есть два важных аспекта, которые не отражены в приведенном выше уравнении логистической функции (Taylor & Letham 2017). Во-первых, емкость многих систем не является постоянной. Например, число людей, имеющих доступ в Интернет (и, соответственно, количество потенциальных пользователей того или иного онлайн-ресурса) со временем возрастает. Поэтому в пакете Prophet постоянная емкость системы $C$ заменяется на изменяющуюся во времени $C(t)$.

Во-вторых, непостоянной обычно бывает и скорость роста $k$. Например, выход новой версии продукта может значительно ускорить рост числа его потребителей. В Prophet для моделирования такого рода изменений в тренде временного ряда вводится понятие точек излома. Предположим, что есть $S$ таких точек на временных отметках $s_j$, $j = 1, 2, \dots, S$. Совокупность всех изменений скорости роста $\delta_j$ можно представить в виде вектора $\boldsymbol \delta \in \mathbb{R}^S$. Тогда скорость роста в любой точке времени $t$ будет равна сумме базовой скорости $k$ и всех изменений, предшествовавших этой точке: $k + \sum_{j:t \gt s_j} \delta_j$. Более наглядно это можно представить с помощью такого вектора $\boldsymbol{a}(t) \in \left\{0, 1\right\}^S$, что \[a_j(t) = \begin{cases} 1 \; \text{если} \; t \ge s_j \\ 0 \; \text{в остальных случаях}. \end{cases}\] Тогда скорость роста в момент времени $t$ равна $k + \boldsymbol{a}(t)^\intercal \boldsymbol{\delta}$. При изменении скорости роста необходимо также изменить параметр сдвига $m$ (см. выше), чтобы обеспечить гладкий стык сегментов кривой тренда на соответствующей временной отметке. Такая поправка рассчитывается следующим образом: \[\gamma_j = \left(s_j - m -\sum_{l \lt j}\gamma_l \right) \left(1 - \frac{k + \sum_{l \lt j} \delta_l}{k + \sum_{l \le j} \delta_l} \right).\] Тогда кусочная логистическая функция (piecewise logistic growth model), которая используется в пакете Prophet для моделирования тренда с насыщением (saturating growth) принимает вид: \[g(t) = \frac{C(t)}{1 + \exp( -(k + \boldsymbol{a}(t)^\intercal \boldsymbol{\delta})(t - (m + \boldsymbol{a}(t)^\intercal \boldsymbol{\gamma})) )}.\] Важным параметром в приведенном уравнении является $C(t)$ - емкость системы в каждый момент времени. Как задать этот параметр для подгонки модели? В большинстве случаев у исследователя уже будет хорошее представление о емкости моделируемой системы (например, компании обычно хорошо знают размер рынка, на котором они работают). Если же такого понимания нет, то придется потратить некоторое время на поиск дополнительной информации в сторонних источниках (например, прогнозы роста численности населения от Всемирного Банка, тренды в Google-запросах и т.п.).

В случаях, когда рост моделируемой системы неограничен, вместо представленной выше логистической функции для моделирования тренда в пакете Prophet используется кусочно-линейная функция следующего вида (обозначения те же): \[ g(t) = (k + \boldsymbol{a}(t)^\intercal \boldsymbol{\delta})t + (m + \boldsymbol{a}(t)^\intercal \boldsymbol{\gamma}).\] Именно такого рода рост по умолчанию предполагался во всех моделях, которые мы рассматривали в предыдущих сообщениях из этой серии.

Примеры моделей

Для начала построим модель, в которой рост зависимой переменной ограничен по некоторому верхнему порогу. Для этого воспользуемся данными по стоимости биткоина за период с января 2016 г. по середину декабря 2017 г., когда наблюдался практически экспоненциальный рост этой переменной (рис. 2):


library(prophet)
library(ggplot2)
library(dplyr)

# Загрузка ранее сохраненных объектов:
load("./workspaces/intro_to_prophet.RData")

# Ограничим обучающие данные периодом до середины декабря 2017 г.:
train_df_short <- train_df %>% filter(ds <= as.Date("2017-12-15"))

# График роста стоимости биткоина:
ggplot(train_df_short, aes(ds, y)) + geom_point()

Рис. 2

Предположим, что несмотря на свой экспоненциальный характер роста стоимость биткоина в будущем не могла превысить 11.5 (на логарифмической шкале, что соответствует почти 99000$ (!) на исходной шкале). Для введения этого верхнего порога в модель необходимо добавить новый столбец cap в таблицу с обучающими данными:


train_df_short$cap <- 11.5

Теперь подгоним модель по этим обучающим данным, рассчитаем прогноз на следующие 180 дней и изобразим полученный результат (рис. 3):


# Подгонка модели:
M17 <- prophet(train_df_short, growth = "logistic",
               changepoint.range = 0.95,
               changepoint.prior.scale = 0.15)

# Таблица с датами прогнозного периода в 180 дней:
future_df_short <- make_future_dataframe(M17, periods = 180) %>% 
    mutate(cap = 11.5) # верхний порог должен быть добавлен и в эту таблицу

# Расчет прогноза:
forecast_M17 <- predict(M17, future_df_short)

# Визуализация подогнанной модели и полученного прогноза:
plot(M17, forecast_M17)

Рис. 3

Как видим, несмотря на экспоненциальный рост в историческом периоде, предсказанные значения стоимости биткоина постепенно выходят на плато, что обусловлено введенным в модель ограничением на рост.

Многие естественные системы имеют ограничения емкости не только по верхнему, но и нижнему порогу (подобно упомянутой выше численности населения). В Prophet имеется возможность учесть это обстоятельство путем добавления еще одного столбца - floor - в таблицы с обучающими данными и данными для расчета прогноза (верхний порог при этом также должен присутствовать в обеих таблицах). Для демонстрации этой возможности воспользуемся данными по стоимости биткоина в 2018 г., когда наблюдался выраженный тренд на ее снижение (рис. 4):


# Обучающие данные за 2018 г.:
train_df_floor <- train_df %>% 
    filter(ds >= as.Date("2018-01-01") & ds <= as.Date("2018-12-31"))

# Визуализация данных за 2018 г.:
ggplot(train_df_floor, aes(ds, y)) + geom_point()

Рис. 4

Построим модель с верхним пределом cap = 10 и нижним пределом floor = 7 (см. полученный результат на рис. 5):


# Добавляем верхний и нижний пороги стоимости биткоина в обучающие данные:
train_df_both$floor <- 7
train_df_both$cap <- 10.0

# Подгонка модели:
M18 <- prophet(train_df_both, growth = "logistic",
               changepoint.range = 0.85,
               changepoint.prior.scale = 0.15)

# Таблица с датами прогнозного периода в 180 дней:
future_df_both <- make_future_dataframe(M18, periods = 180) %>% 
    mutate(floor = 7, cap = 10)

# Расчет прогноза:
forecast_M18 <- predict(M18, future_df_both)

# Визуализация подогнанной модели и полученного прогноза:
plot(M18, forecast_M18)

Рис. 5

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

***

Этим сообщением я завершаю описание основных возможностей пакета Prophet для прогнозирования временных рядов. Другие статьи из этой серии можно найти по следующим ссылкам: