Отстающие признаки для прогнозирования временных рядов

Анализ набора данных о спросе на велосипеды

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

Адрес URL файла parquet можно найти в описании 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://openml1.win.tue.nl/datasets/0004/44063/dataset_44063.pq",
    sha256="d120af76829af0d256338dc6dd4be5df4fd1f35bf3a283cab66a51c1c6abd06a",
)
bike_sharing_data_file

PosixPath('/home/circleci/scikit_learn_data/openml1.win.tue.nl/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.3889873516666431

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

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

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.44300751539296973

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

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.44300752, 0.27772182, 0.3697178 ])

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

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

CV MAPE: 0.363 ± 0.068

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

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\), как это делают методы наименьших квадратов и Пуассона, можно попытаться оценить квантили условного распределения.

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

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

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 минут 9,647 секунд)

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

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

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

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

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

Квантильная регрессия