В данной статье я расскажу о нескольких полезных способах улучшения качества предсказаний модели Facebook Prophet. В этой статье я не буду уделять внимание таким подходам, как denoising или различным трансформациям (например, box-cox), а сконцентрируюсь на двух методах, которые позволят получить лучшее качество, не внося изменения в данные.
На Хабре уже есть подробные обзоры Prophet (ссылка), поэтому мы ограничимся кратким введением. Prophet — это аддитивная линейная модель (GAM) для прогнозирования временных рядов, которая объединяет несколько компонентов для выявления различных закономерностей в данных (тренд, сезонность, праздники и т.д.). Даже спустя несколько лет после релиза Prophet остается одним из лучших решений для одномерного прогнозирования (univariate forecasting).
Одной из ключевых особенностей Prophet является то, что он генерирует не только прогнозы, но и признаки как для тренировочных, так и для тестовых данных. Поэтому возникает идея: «Почему бы не использовать эти данные в какой‑нибудь другой модели?». Это похоже на стекинг, где результаты первой модели улучшаются мета-моделью. Такой подход достаточно хорошо известен и иногда используется, например, в Microsoft Azure ML.
Продолжение этого подхода — обучение нескольких Prophet'ов с использованием данных, агрегированных на разных временных шкалах. Например, если у нас есть данные с ежедневной частотой, они могут быть сильно зашумленными, что затрудняет выявление более общих паттернов (недельных или месячных). Чтобы преодолеть эту проблему, можно агрегировать исходные данные на более высоком временном уровне, используя сумму, среднее или медиану, а затем обучить отдельные модели Prophet на этих данных. Таким образом можно сгенерировать дополнительные признаки для обучения мета-модели.
Достаточно прелюдий, перейдем к рассмотрению всех трех подходов:
«Ванильный» Prophet
Prophet + LightGBM
Улучшение предсказаний Prophet с использованием LightGBM и многоуровневой агрегации
А также сравним качество этих подходов на практике.
В качестве примера мы возьмем набор данных о посещениях определенного веб-сайта с 2014 по 2016 годы в разрезе дня (ссылка). Для обучения выберем все данные, кроме данных за последние 30 дней, они будут использоваться для тестирования. В наборе данных есть несколько столбцов, которые могут быть выбраны в качестве целевой переменной, мы остановимся на общем количестве загрузок страницы — Page.Loads.
import pandas as pd
data['ds'] = pd.to_datetime(data['Date'])
data = data.drop('Date', axis=1)
data = data.rename(columns={'Page.Loads': 'y'})
data['y'] = data['y'].str.replace(',', '').astype(float)
data = data[['ds', 'y']]
train_series = data[data.ds < (data.ds.max() - timedelta(days=30))].copy()
test_series = data[data.ds >= (data.ds.max() - timedelta(days=30))].copy()Prophet требует определенных наименований столбцов, где метки времени должны быть в столбце ds, а целевые значения — в столбце y.
Давайте визуализируем наш временной ряд. Мы видим определенные паттерны, а сами данные имеют некритичный уровень шума. Идеально было бы декомпозировать ряд и провести анализ с помощью тестов перед обучением модели, но в данном случае мы не будем проводить никаких манипуляций с данными и попробуем обучить модель сразу.
Vanilla Prophet
Теперь обучим Prophet с использованием стандартных настроек и оценим качество предсказания. При обучении я использую специальный класс suppress_stdout_stderr для отключения всех логов, которые генерирует C-шный бэкэнд. Код этого класса я позаимствовал у одного доброго человека в разделе Issues на официальном GitHub Prophet (ссылка).
import prophet
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
fb = prophet.Prophet()
with suppress_stdout_stderr():
fb.fit(train_series)
predictions = fb.make_future_dataframe(periods=len(test_series), freq='D')
forecast = fb.predict(predictions)
v_fb_df = test_series.copy()
v_fb_df = v_fb_df.merge(forecast[['ds', 'yhat']], on='ds', how='left')np.sqrt(mean_squared_error(v_fb_df['y'], v_fb_df['yhat']))
>>> 512.2325266529127mean_absolute_percentage_error(v_fb_df['y'], v_fb_df['yhat'])
>>> 0.1360530885720578На первый взгляд получаем довольно неплохие результаты: RMSE = 512.23, MAPE = 0.136. Как видно на графике ниже, модель хорошо улавливает паттерны, но немного занижает значения по сравнению с реальными данными, что означает наличие потенциала для улучшения.
Prophet + LGBM
Теперь давайте попробуем улучшить этот результат с помощью LightGBM. Я буду использовать не стандартный LightGBM, а обертку над ним, которая интегрирована в библиотеку Optuna. В отличие от стандартного подбора гиперпараметров в Optuna, интегрированный метод использует пошаговый подход, оптимизируя гиперпараметры по одному за раз с использованием предварительно определенного пространства поиска. Этот подход менее эффективен, но в отличие от Grid Search, где обучение происходит на комбинации всех заданных параметров, пошаговый подход гораздо быстрее, так как количество возможных параметров сокращается с p_1 * p_2 * p_3 до p_1 + p_2 + p_3, где p_n — это заданный набор значений для конкретного параметра.
Давайте создадим небольшую обертку, чтобы упростить использование в будущем.
import optuna
import optuna.integration.lightgbm as lgb
from optuna.integration.lightgbm import Dataset
import pandas as pd
import numpy as np
class OptunaLGBMRegressor:
"""
A wrapper class for the LightGBM Regressor with Optuna for hyperparameters tuning
"""
def __init__(
self,
n_estimators: int,
learning_rate: float = 0.01,
metric: str = 'rmse',
cat_columns: str = 'auto',
seed: int = 42
):
"""
Initializes a new instance of the OptunaLGBMRegressor class
"""
self.params = {
"n_estimators": n_estimators,
"objective": "regression",
"verbosity": -1,
"metric": metric,
"learning_rate": learning_rate,
"boosting_type": 'gbdt',
"random_state": seed
}
self.cat_columns = cat_columns
self.model = None
self.features = None
self.is_fitted_ = False
def _to_datasets(
self, x_train: pd.DataFrame, y_train: np.ndarray, x_val: pd.DataFrame, y_val: np.ndarray
) -> (Dataset, Dataset):
"""
Converts Pandas DataFrames to LightGBM Datasets
"""
self.features = list(x_train.columns)
X_val = x_val[self.features].copy()
dtrain = Dataset(x_train, label=y_train, categorical_feature=self.cat_columns)
dval = Dataset(X_val, label=y_val, categorical_feature=self.cat_columns)
return dtrain, dval
def fit(self, X_train: pd.DataFrame, y_train: np.ndarray, X_val: pd.DataFrame, y_val: np.ndarray) -> None:
dtrain, dval = self._to_datasets(X_train, y_train, X_val, y_val)
self.model = lgb.tuner.train(
self.params,
dtrain,
valid_sets=[dtrain, dval],
verbose_eval=0,
early_stopping_rounds=150
)
self.is_fitted_ = True
def predict(self, X_test: pd.DataFrame) -> np.ndarray:
assert self.is_fitted_, 'Model is not fitted!'
return self.model.predict(X_test[self.features], num_iteration=self.model.best_iteration)
Прежде чем обучать модель на данных из Prophet, необходимо определить правильное разделение на тренировочный и валидационный наборы, поскольку валидационный набор необходим для контроля переобучения модели. Здесь нет одного правильного ответа. С одной стороны, было бы правильно использовать out-of-time выборку, то есть взять данные после определенной даты для валидации (Time Series Split), но в таком случае часть данных может не попасть в тренировочный или валидационный наборы. В случае, когда у нас есть наблюдения, например, только за один год, мы обучаем модель на определенных месяцах, а валидацию проводим на тех, которые модель не видела. Это может привести к тому, что наш подход будет работать не лучшим образом. С другой стороны, мы можем использовать случайное разделение (random split), но тогда мы неправильно предполагаем независимость наблюдений, что для временных рядов, скорее всего, неверно. Я буду использовать второй подход, так как он показывает меньшее переобучение, но я рекомендую принимать решение на основе конкретного случая.
from sklearn.model_selection import train_test_split
gbt_data = train_series.merge(forecast, on='ds', how='left')
train_gbt, val_gbt = train_test_split(gbt_data, test_size=0.15, random_state=42)
lgbm = OptunaLGBMRegressor(n_estimators=300, learning_rate=0.03, metric='mape', seed=42)
lgbm.fit(
train_gbt.drop(['ds', 'y'], 1), train_gbt.y.values,
val_gbt.drop(['ds', 'y'], 1), val_gbt.y.values
)
test_gbt = test_series.merge(forecast, on='ds', how='left')
preds = lgbm.predict(test_gbt.drop(['ds', 'y'], 1))
forecast_df = test_gbt[['ds', 'y', 'yhat']].copy()
forecast_df['gbt_yhat'] = predsnp.sqrt(mean_squared_error(forecast_df['y'], forecast_df['gbt_yhat']))
>>> 461.3712552827671mean_absolute_percentage_error(forecast_df['y'], forecast_df['gbt_yhat'])
>>> 0.11717199190922199Обучаем модель и проверяем результаты. Обучение с использованием LightGBM позволило нам улучшить RMSE до 461.4 и MAPE до 0.117. Как видно, комбинируя широкое охватывание паттернов Prophet с способностью LightGBM обнаруживать сложные взаимосвязи в данных, мы можем получать более точные прогнозы.
Но это еще не предел!
Добавление многоуровневой агрегации
Используя только Prophet и LightGBM, мы в значительной мере полагаемся на стабильность исходного временного ряда. Обычно увеличение гранулярности позволяет улучшить его стабильность. Не секрет, что прогнозирование недельных или месячных данных обычно гораздо проще, чем минутных или дневных временных рядов.
Однако перед тем, как проверять это утверждение, нам нужно написать класс, который будет:
Принимать входные данные и агрегировать их на указанных частотностях с использованием заданных агрегирующих функций.
Обучать объекты класса Prophet для каждой комбинации «частотность — агрегирующая функция».
Прогнозировать временные ряды на указанном горизонте.
Объединять полученные данные в один набор данных.
Удалять неинформативные признаки.
(Дополнительно я добавил поддержку праздников, если кто-то захочет протестировать этот подход на своих данных. Все, что нужно сделать, это передать соответствующий объект (для нужной страны) из библиотеки holidays).
В результате получился не самый компактный, но очень полезный класс, который берет на себя рутину по генерации признаков с использованием Prophet.
import re
import pandas as pd
from typing import Optional
from itertools import product
from tqdm import tqdm
from holidays.holiday_base import HolidayBase
from prophet import Prophet
class ProphetsEnsemble:
"""An ensemble of Prophet models with different aggregation functions and frequencies."""
def __init__(self, freq: str, levels: list, agg_fn: list, holidays_getter: HolidayBase = None):
"""Initializes an ensemble of Prophet models."""
self.freq = freq
self.levels = ['_'.join(x) for x in product(levels, agg_fn)]
self.h_getter = holidays_getter
self.prophets_ = dict()
self.is_fitted_ = False
@staticmethod
def _resample(data: pd.DataFrame, freq: str, how: str) -> pd.DataFrame:
"""Resamples a time series DataFrame."""
if how not in ['median', 'mean', 'sum']:
raise NotImplementedError(f'Unknown function {how}. Only [median, mean, sum] are supported.')
return data.set_index('ds').resample(freq).agg(how).reset_index(drop=False)
@staticmethod
def _merge_key_gen(x, level: str) -> str:
"""Generates a key for merging DataFrames based on the frequency."""
freq = re.sub('[\d]', '', level.split('_')[0])
if freq == 'H':
return f'{x.year}-{x.month}-{x.day}-{x.hour}'
elif freq in ['D', 'M']:
return f'{x.year}-{x.month}-{x.day}' if freq == 'D' else f'{x.year}-{x.month}'
elif freq == 'W':
return f'{x.year}-{x.week}'
raise NotImplementedError(f'Only [H, D, W, M] are supported. {freq} was recieved as input!')
def _get_holidays(self, data: pd.DataFrame) -> Optional[pd.DataFrame]:
"""Extracts holidays from the data."""
if self.h_getter is None:
return None
holidays = data[['ds']].copy()
holidays['holiday'] = holidays['ds'].apply(self.h_getter.get)
return holidays.dropna()
def _fit_level(self, data: pd.DataFrame, level: str) -> None:
"""Fits a Prophet model for a specific aggregation level."""
resampled = self._resample(data, *level.split('_')) if level != self.freq else data.copy()
fb = Prophet(holidays=self._get_holidays(resampled))
with suppress_stdout_stderr():
fb.fit(resampled)
self.prophets_[level] = fb
def _predict_level(self, periods: int, level: str) -> pd.DataFrame:
"""Makes predictions for a specific aggregation level."""
fb = self.prophets_[level]
df = fb.make_future_dataframe(periods=periods, freq=level.split('_')[0])
forecasts = fb.predict(df)
forecasts.columns = [f'{x}_{level}' for x in forecasts.columns]
return forecasts
def _combine_levels(self, base_df: pd.DataFrame, data: pd.DataFrame, level: str) -> pd.DataFrame:
"""Combines predictions from different aggregation levels."""
key = lambda x: self._merge_key_gen(x, level)
return (
base_df.assign(key=base_df['ds'].apply(key))
.merge(data.assign(key=data[f'ds_{level}'].apply(key)), on='key', how='left')
.drop(['key', f'ds_{level}'], axis=1)
)
@staticmethod
def _drop_redundant(data: pd.DataFrame) -> pd.DataFrame:
"""Drops redundant features from the DataFrame."""
redundant = [col for col in data.columns if col != 'ds' and 'yhat' not in col and len(data[col].unique()) == 1]
return data.drop(redundant, axis=1)
def fit(self, data: pd.DataFrame) -> None:
"""Fits the Prophet models for all aggregation levels."""
for level in tqdm([self.freq] + self.levels, 'Fitting prophets...'):
self._fit_level(data, level)
self.is_fitted_ = True
def forecast(self, periods: int) -> pd.DataFrame:
"""Makes forecasts for all aggregation levels and combines them."""
assert self.is_fitted_, 'Model is not fitted'
forecasts = [self._predict_level(periods, level) for level in tqdm([self.freq] + self.levels, 'Forecasting...')]
forecast = forecasts[0].rename(columns={f'ds_{self.freq}': 'ds'})
for level, fore in zip(self.levels, forecasts[1:]):
forecast = self._combine_levels(forecast, fore, level)
return self._drop_redundant(forecast)Теперь попробуем сгенерировать данные и обучить нашу модель LightGBM+Optuna.
pe = ProphetsEnsemble(freq='D', levels=['W', 'M'], agg_fn=['mean', 'median'])
pe.fit(train_series)
pe_forecast = pe.forecast(len(test_series))
gbt_data = train_series.merge(pe_forecast, on='ds', how='left')
train_gbt, val_gbt = train_test_split(gbt_data, test_size=0.15, random_state=42)
lgbm = OptunaLGBMRegressor(n_estimators=300, learning_rate=0.03, metric='mape', seed=42)
lgbm.fit(
train_gbt.drop(['ds', 'y'], 1), train_gbt.y.values,
val_gbt.drop(['ds', 'y'], 1), val_gbt.y.values
)
test_gbt = test_series.merge(pe_forecast, on='ds', how='left')
preds = lgbm.predict(test_gbt.drop(['ds', 'y'], 1))
forecast_df['pe_yhat'] = predsnp.sqrt(mean_squared_error(forecast_df['y'], forecast_df['pe_yhat']))
>>> 429.81662049975967mean_absolute_percentage_error(forecast_df['y'], forecast_df['pe_yhat'])
>>> 0.09649907311310414Как видно, метод сработал и позволил еще повысить качество предсказаний. На графике видно, что итоговая модель практически перестала занижать значения в минимуме, а среди всех подходов именно последний ближе всего к реальным значениям.
Для предсказания всего на один шаг вперёд данный подход может быть немного громоздким или даже избыточным, однако он хорошо себя показывает, когда нужно предсказывать данные на большой горизонт.
Надеюсь, описанные методы помогут вам в решении задач прогнозирования временных рядов. Полный ноутбук доступен по ссылке: GitHub.
Спасибо за внимание!
* Facebook принадлежит компании Meta, признанной экстремистской организацией и запрещенной в РФ