Лаговые признаки для прогнозирования временных рядов

Анализ набора данных Bike Sharing Demand

Для начала загрузим данные из репозитория OpenML в виде исходного файла Parquet, чтобы показать, как работать с произвольным файлом Parquet, а не скрывать этот шаг за удобным инструментом, таким как sklearn.datasets.fetch_openml.

URL-адрес файла Parquet можно найти в описании набора данных Bike Sharing Demand в формате JSON с идентификатором 44063 на openml.org (https://openml.org/search?type=data&status=active&id=44063).

Для проверки целостности загруженного файла также указан его хеш sha256.

import numpy as np
import polars as pl

from sklearn.datasets import fetch_file

pl.Config.set_fmt_str_lengths(20)

bike_sharing_data_file = fetch_file(
    "https://data.openml.org/datasets/0004/44063/dataset_44063.pq",
    sha256="d120af76829af0d256338dc6dd4be5df4fd1f35bf3a283cab66a51c1c6abd06a",
)
bike_sharing_data_file

PosixPath('/home/circleci/scikit_learn_data/data.openml.org/datasets_0004_44063/dataset_44063.pq')

Загрузим файл Parquet с помощью Polars для создания признаков. Polars автоматически кэширует общие подвыражения, которые используются в нескольких выражениях (например, pl.col("count").shift(1) ниже). Дополнительные сведения см. на странице https://docs.pola.rs/user-guide/lazy/optimizations/.

df = pl.read_parquet(bike_sharing_data_file)

Теперь рассмотрим статистическую сводку набора данных, чтобы лучше понять данные, с которыми мы работаем.

import polars.selectors as cs

summary = df.select(cs.numeric()).describe()
summary

Посмотрим на количество наблюдений для сезонов "fall", "spring", "summer" и "winter" в наборе данных, чтобы убедиться, что они сбалансированы.

import matplotlib.pyplot as plt

df["season"].value_counts()

Создание лаговых признаков с помощью Polars

Рассмотрим задачу прогнозирования спроса на следующий час на основе спроса в прошлом. Поскольку спрос — непрерывная переменная, интуитивно можно было бы использовать любую модель регрессии. Однако у нас нет обычного набора данных (X_train, y_train). Вместо этого у нас есть только данные о спросе y_train, последовательно упорядоченные по времени.

lagged_df = df.select(
    "count",
    *[pl.col("count").shift(i).alias(f"lagged_count_{i}h") for i in [1, 2, 3]],
    lagged_count_1d=pl.col("count").shift(24),
    lagged_count_1d_1h=pl.col("count").shift(24 + 1),
    lagged_count_7d=pl.col("count").shift(7 * 24),
    lagged_count_7d_1h=pl.col("count").shift(7 * 24 + 1),
    lagged_mean_24h=pl.col("count").shift(1).rolling_mean(24),
    lagged_max_24h=pl.col("count").shift(1).rolling_max(24),
    lagged_min_24h=pl.col("count").shift(1).rolling_min(24),
    lagged_mean_7d=pl.col("count").shift(1).rolling_mean(7 * 24),
    lagged_max_7d=pl.col("count").shift(1).rolling_max(7 * 24),
    lagged_min_7d=pl.col("count").shift(1).rolling_min(7 * 24),
)
lagged_df.tail(10)

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

lagged_df.head(10)

Теперь можно выделить лаговые признаки в матрицу X, а целевую переменную (прогнозируемые значения) — в массив с таким же размером первого измерения y.

lagged_df = lagged_df.drop_nulls()
X = lagged_df.drop("count")
y = lagged_df["count"]
print("X shape: {}\ny shape: {}".format(X.shape, y.shape))

X shape: (17210, 13)
y shape: (17210,)

Наивная оценка регрессии для прогнозирования спроса на велосипеды в следующий час

Случайным образом разделим наш табличный набор данных на обучающую и тестовую выборки, чтобы обучить модель регрессионного дерева с градиентным бустингом (GBRT) и оценить её с помощью средней абсолютной процентной ошибки (MAPE). Если наша модель предназначена для прогнозирования (то есть для предсказания будущих данных на основе прошлых), не следует использовать обучающие данные, относящиеся к более позднему периоду, чем тестовые данные. В машинном обучении для временных рядов предположение «i.i.d» (независимость и одинаковое распределение) не выполняется, поскольку точки данных не являются независимыми и связаны временными отношениями.

from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

model = HistGradientBoostingRegressor().fit(X_train, y_train)

Рассмотрим качество работы модели.

from sklearn.metrics import mean_absolute_percentage_error

y_pred = model.predict(X_test)
mean_absolute_percentage_error(y_test, y_pred)

0.3951983212888564

Корректная оценка прогнозирования на следующий час

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

from sklearn.model_selection import TimeSeriesSplit

ts_cv = TimeSeriesSplit(
    n_splits=3,  # to keep the notebook fast enough on common laptops
    gap=48,  # 2 days data gap between train and test
    max_train_size=10000,  # keep train sets of comparable sizes
    test_size=3000,  # for 2 or 3 digits of precision in scores
)
all_splits = list(ts_cv.split(X, y))

Обучим модель и оценим её качество на основе MAPE.

train_idx, test_idx = all_splits[0]
X_train, X_test = X[train_idx, :], X[test_idx, :]
y_train, y_test = y[train_idx], y[test_idx]

model = HistGradientBoostingRegressor().fit(X_train, y_train)
y_pred = model.predict(X_test)
mean_absolute_percentage_error(y_test, y_pred)

0.45822977660508035

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

from sklearn.model_selection import cross_val_score

cv_mape_scores = -cross_val_score(
    model, X, y, cv=ts_cv, scoring="neg_mean_absolute_percentage_error"
)
cv_mape_scores

array([0.45822978, 0.2825124 , 0.36374485])

Разброс значений между разбиениями довольно велик! В реальных условиях рекомендуется использовать больше разбиений, чтобы точнее оценить вариативность. Далее будем приводить средние значения метрик перекрёстной проверки и их стандартное отклонение.

print(f"CV MAPE: {cv_mape_scores.mean():.3f} ± {cv_mape_scores.std():.3f}")

CV MAPE: 0.368 ± 0.072

Можно вычислить несколько комбинаций метрик оценки и функций потерь; они приведены ниже.

from collections import defaultdict

from sklearn.metrics import (
    make_scorer,
    mean_absolute_error,
    mean_pinball_loss,
    root_mean_squared_error,
)
from sklearn.model_selection import cross_validate


def consolidate_scores(cv_results, scores, metric):
    if metric == "MAPE":
        scores[metric].append(f"{value.mean():.2f} ± {value.std():.2f}")
    else:
        scores[metric].append(f"{value.mean():.1f} ± {value.std():.1f}")

    return scores


scoring = {
    "MAPE": make_scorer(mean_absolute_percentage_error),
    "RMSE": make_scorer(root_mean_squared_error),
    "MAE": make_scorer(mean_absolute_error),
    "pinball_loss_05": make_scorer(mean_pinball_loss, alpha=0.05),
    "pinball_loss_50": make_scorer(mean_pinball_loss, alpha=0.50),
    "pinball_loss_95": make_scorer(mean_pinball_loss, alpha=0.95),
}
loss_functions = ["squared_error", "poisson", "absolute_error"]
scores = defaultdict(list)
for loss_func in loss_functions:
    model = HistGradientBoostingRegressor(loss=loss_func)
    cv_results = cross_validate(
        model,
        X,
        y,
        cv=ts_cv,
        scoring=scoring,
        n_jobs=2,
    )
    time = cv_results["fit_time"]
    scores["loss"].append(loss_func)
    scores["fit_time"].append(f"{time.mean():.2f} ± {time.std():.2f} s")

    for key, value in cv_results.items():
        if key.startswith("test_"):
            metric = key.split("test_")[1]
            scores = consolidate_scores(cv_results, scores, metric)

Моделирование неопределённости прогнозов с помощью квантильной регрессии

Вместо моделирования ожидаемого значения распределения \(Y|X\), как это делают функции потерь наименьших квадратов и Пуассона, можно попробовать оценить квантили условного распределения.

Для заданной точки данных \(x_i\) предполагается, что \(Y|X=x_i\) является случайной величиной, поскольку мы ожидаем, что количество аренд невозможно предсказать по признакам со стопроцентной точностью. На него могут влиять другие переменные, которые не учитываются должным образом существующими лаговыми признаками. Например, по данным об аренде велосипедов за предыдущие часы нельзя в полной мере предсказать, будет ли дождь в следующий час. Это и называется алеаторной неопределённостью.

Квантильная регрессия позволяет точнее описать это распределение, не делая строгих предположений о его форме.

quantile_list = [0.05, 0.5, 0.95]

for quantile in quantile_list:
    model = HistGradientBoostingRegressor(loss="quantile", quantile=quantile)
    cv_results = cross_validate(
        model,
        X,
        y,
        cv=ts_cv,
        scoring=scoring,
        n_jobs=2,
    )
    time = cv_results["fit_time"]
    scores["fit_time"].append(f"{time.mean():.2f} ± {time.std():.2f} s")

    scores["loss"].append(f"quantile {int(quantile * 100)}")
    for key, value in cv_results.items():
        if key.startswith("test_"):
            metric = key.split("test_")[1]
            scores = consolidate_scores(cv_results, scores, metric)

scores_df = pl.DataFrame(scores)
scores_df

Рассмотрим функции потерь, которые минимизируют каждую метрику.

def min_arg(col):
    col_split = pl.col(col).str.split(" ")
    return pl.arg_sort_by(
        col_split.list.get(0).cast(pl.Float64),
        col_split.list.get(2).cast(pl.Float64),
    ).first()


scores_df.select(
    pl.col("loss").get(min_arg(col_name)).alias(col_name)
    for col_name in scores_df.columns
    if col_name != "loss"
)

Несмотря на то что распределения оценок перекрываются из-за дисперсии в наборе данных, среднее значение RMSE действительно ниже, когда loss="squared_error", тогда как среднее значение MAPE ниже, когда loss="absolute_error", как и ожидалось. То же самое справедливо для средней функции потерь Pinball при квантилях 5 и 95. Оценка для функции потерь при квантиле 50 перекрывается с оценками, полученными при минимизации других функций потерь; то же верно и для MAE.

Качественный анализ прогнозов

Теперь можно визуализировать качество модели для 5-го процентиля, медианы и 95-го процентиля:

all_splits = list(ts_cv.split(X, y))
train_idx, test_idx = all_splits[0]

X_train, X_test = X[train_idx, :], X[test_idx, :]
y_train, y_test = y[train_idx], y[test_idx]

max_iter = 50
gbrt_mean_poisson = HistGradientBoostingRegressor(loss="poisson", max_iter=max_iter)
gbrt_mean_poisson.fit(X_train, y_train)
mean_predictions = gbrt_mean_poisson.predict(X_test)

gbrt_median = HistGradientBoostingRegressor(
    loss="quantile", quantile=0.5, max_iter=max_iter
)
gbrt_median.fit(X_train, y_train)
median_predictions = gbrt_median.predict(X_test)

gbrt_percentile_5 = HistGradientBoostingRegressor(
    loss="quantile", quantile=0.05, max_iter=max_iter
)
gbrt_percentile_5.fit(X_train, y_train)
percentile_5_predictions = gbrt_percentile_5.predict(X_test)

gbrt_percentile_95 = HistGradientBoostingRegressor(
    loss="quantile", quantile=0.95, max_iter=max_iter
)
gbrt_percentile_95.fit(X_train, y_train)
percentile_95_predictions = gbrt_percentile_95.predict(X_test)

Теперь рассмотрим прогнозы, полученные с помощью регрессионных моделей:

last_hours = slice(-96, None)
fig, ax = plt.subplots(figsize=(15, 7))
plt.title("Predictions by regression models")
ax.plot(
    y_test[last_hours],
    "x-",
    alpha=0.2,
    label="Actual demand",
    color="black",
)
ax.plot(
    median_predictions[last_hours],
    "^-",
    label="GBRT median",
)
ax.plot(
    mean_predictions[last_hours],
    "x-",
    label="GBRT mean (Poisson)",
)
ax.fill_between(
    np.arange(96),
    percentile_5_predictions[last_hours],
    percentile_95_predictions[last_hours],
    alpha=0.3,
    label="GBRT 90% interval",
)
_ = ax.legend()

Здесь интересно отметить, что ширина синей области между оценками 5-го и 95-го процентилей меняется в зависимости от времени суток:

• Ночью синяя полоса гораздо уже: пара моделей с высокой уверенностью предсказывает небольшое количество аренд велосипедов. Более того, это, по-видимому, верные прогнозы, поскольку фактический спрос остаётся в пределах этой синей полосы.
• Днём синяя полоса гораздо шире: неопределённость растёт, вероятно, из-за изменчивости погоды, которая может оказывать очень большое влияние, особенно в выходные.
• Также видно, что в будние дни характер поездок на работу по-прежнему заметен в оценках 5-го и 95-го процентилей.
• Наконец, ожидается, что в 10% случаев фактический спрос не будет находиться между оценками 5-го и 95-го процентилей. На этом тестовом интервале фактический спрос, по-видимому, выше, особенно в часы пик. Это может указывать на то, что оценка 95-го процентиля недооценивает пиковый спрос. Это можно количественно подтвердить, вычислив эмпирические значения покрытия, как описано в разделе калибровка доверительных интервалов.

Сравним качество нелинейных регрессионных моделей с лучшими моделями:

from sklearn.metrics import PredictionErrorDisplay

fig, axes = plt.subplots(ncols=3, figsize=(15, 6), sharey=True)
fig.suptitle("Non-linear regression models")
predictions = [
    median_predictions,
    percentile_5_predictions,
    percentile_95_predictions,
]
labels = [
    "Median",
    "5th percentile",
    "95th percentile",
]
for ax, pred, label in zip(axes, predictions, labels):
    PredictionErrorDisplay.from_predictions(
        y_true=y_test,
        y_pred=pred,
        kind="residual_vs_predicted",
        scatter_kwargs={"alpha": 0.3},
        ax=ax,
    )
    ax.set(xlabel="Predicted demand", ylabel="True demand")
    ax.legend(["Best model", label])

plt.show()

Заключение

В этом примере мы изучили прогнозирование временных рядов с использованием лаговых признаков. Мы сравнили наивную регрессию (с использованием стандартизированного разбиения данных train_test_split) с корректной стратегией оценки временных рядов с помощью TimeSeriesSplit. Мы увидели, что модель, обученная с помощью train_test_split со значением параметра по умолчанию shuffle, установленным в True, показала чрезмерно оптимистичную среднюю абсолютную процентную ошибку (MAPE). Результаты, полученные при разбиении с учётом времени, лучше отражают качество нашей модели регрессии временных рядов. Мы также проанализировали неопределённость прогнозов нашей модели с помощью квантильной регрессии. Прогнозы на основе 5-го и 95-го процентилей, полученные с помощью loss="quantile", дают количественную оценку неопределённости прогнозов нашей модели регрессии временных рядов. Оценить неопределённость также можно с помощью MAPIE, в котором реализованы методы, основанные на последних исследованиях в области конформного предсказания и позволяющие одновременно оценивать алеаторную и эпистемическую неопределённость. Кроме того, функциональность sktime позволяет расширять оценщики scikit-learn за счёт рекурсивного прогнозирования временных рядов, обеспечивающего динамические прогнозы будущих значений.

Общее время выполнения скрипта: (0 минут 8.792 секунды)

Download Jupyter notebook: plot_time_series_lagged_features.ipynb

Download Python source code: plot_time_series_lagged_features.py

Download zipped: plot_time_series_lagged_features.zip

Связанные примеры

Баланс между сложностью модели и оценкой при кросс-валидации

Создание признаков, связанных со временем

Интервалы прогнозирования для регрессии с градиентным бустингом

Признаки деревьев с гистограммным градиентным бустингом