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

Исследование данных набора Bike Sharing Demand

Сначала загрузим данные из репозитория OpenML.

from sklearn.datasets import fetch_openml

bike_sharing = fetch_openml("Bike_Sharing_Demand", version=2, as_frame=True)
df = bike_sharing.frame

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

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

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(12, 4))
average_week_demand = df.groupby(["weekday", "hour"])["count"].mean()
average_week_demand.plot(ax=ax)
_ = ax.set(
    title="Average hourly bike demand during the week",
    xticks=[i * 24 for i in range(7)],
    xticklabels=["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"],
    xlabel="Time of the week",
    ylabel="Number of bike rentals",
)

Целевая переменная задачи прогнозирования — абсолютное количество аренд велосипедов за каждый час:

df["count"].max()

np.int64(977)

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

y = df["count"] / df["count"].max()

fig, ax = plt.subplots(figsize=(12, 4))
y.hist(bins=30, ax=ax)
_ = ax.set(
    xlabel="Fraction of rented fleet demand",
    ylabel="Number of hours",
)

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

X = df.drop("count", axis="columns")
X

Теперь рассмотрим распределение категориальных переменных, начав с "weather":

X["weather"].value_counts()

weather
clear         11413
misty          4544
rain           1419
heavy_rain        3
Name: count, dtype: int64

Поскольку событий "heavy_rain" всего 3, мы не можем использовать эту категорию для обучения моделей машинного обучения с перекрестной проверкой. Вместо этого мы упростим представление, объединив их с категорией "rain".

X["weather"] = (
    X["weather"]
    .astype(object)
    .replace(to_replace="heavy_rain", value="rain")
    .astype("category")
)

X["weather"].value_counts()

weather
clear    11413
misty     4544
rain      1422
Name: count, dtype: int64

Как и ожидалось, переменная "season" имеет сбалансированное распределение:

X["season"].value_counts()

season
fall      4496
summer    4409
spring    4242
winter    4232
Name: count, dtype: int64

Перекрестная проверка с учетом времени

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

1000 тестовых точек данных должно быть достаточно для оценки эффективности модели. Это чуть меньше полутора месяцев непрерывных тестовых данных:

from sklearn.model_selection import TimeSeriesSplit

ts_cv = TimeSeriesSplit(
    n_splits=5,
    gap=48,
    max_train_size=10000,
    test_size=1000,
)

Давайте вручную проверим различные разбиения, чтобы убедиться, что TimeSeriesSplit работает так, как мы ожидаем. Начнем с первого разбиения:

all_splits = list(ts_cv.split(X, y))
train_0, test_0 = all_splits[0]

X.iloc[test_0]

X.iloc[train_0]

Теперь рассмотрим последнее разбиение:

train_4, test_4 = all_splits[4]

X.iloc[test_4]

X.iloc[train_4]

Все в порядке. Теперь можно приступать к построению прогностических моделей!

Градиентный бустинг

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

Здесь мы используем современный HistGradientBoostingRegressor со встроенной поддержкой категориальных признаков. Поэтому нам нужно лишь установить categorical_features="from_dtype", чтобы признаки с категориальным dtype считались категориальными. Для справки мы извлекаем категориальные признаки из фрейма данных на основе dtype. Для этих признаков внутренние деревья используют специальное правило разбиения.

Числовые переменные не требуют предварительной обработки, и для простоты мы пробуем только гиперпараметры модели по умолчанию:

from sklearn.compose import ColumnTransformer
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import cross_validate
from sklearn.pipeline import make_pipeline

gbrt = HistGradientBoostingRegressor(categorical_features="from_dtype", random_state=42)
categorical_columns = X.columns[X.dtypes == "category"]
print("Categorical features:", categorical_columns.tolist())

Categorical features: ['season', 'holiday', 'workingday', 'weather']

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

import numpy as np


def evaluate(model, X, y, cv, model_prop=None, model_step=None):
    cv_results = cross_validate(
        model,
        X,
        y,
        cv=cv,
        scoring=["neg_mean_absolute_error", "neg_root_mean_squared_error"],
        return_estimator=True,
    )
    if model_prop is not None:
        if model_step is not None:
            values = [
                getattr(m[model_step], model_prop) for m in cv_results["estimator"]
            ]
        else:
            values = [getattr(m, model_prop) for m in cv_results["estimator"]]
        print(f"Mean model.{model_prop} = {np.mean(values)}")
    mae = -cv_results["test_neg_mean_absolute_error"]
    rmse = -cv_results["test_neg_root_mean_squared_error"]
    print(
        f"Mean Absolute Error:     {mae.mean():.3f} +/- {mae.std():.3f}\n"
        f"Root Mean Squared Error: {rmse.mean():.3f} +/- {rmse.std():.3f}"
    )
    # To display the fitted estimator diagrams in the notebook.
    return cv_results["estimator"][0]


evaluate(gbrt, X, y, cv=ts_cv, model_prop="n_iter_")

Mean model.n_iter_ = 100.0
Mean Absolute Error:     0.044 +/- 0.003
Root Mean Squared Error: 0.068 +/- 0.005

HistGradientBoostingRegressor(random_state=42)

Мы видим, что задали max_iter достаточно большим, чтобы сработала ранняя остановка.

Средняя ошибка этой модели составляет около 4–5% от максимального спроса. Это очень хороший результат для первой попытки без настройки гиперпараметров! Нам нужно было лишь явно указать категориальные переменные. Обратите внимание, что временные признаки передаются без изменений, то есть без предварительной обработки. Однако для моделей на основе деревьев это не представляет большой проблемы: они могут обучаться немонотонной зависимости между порядковыми входными признаками и целевой переменной.

Как мы увидим далее, для моделей линейной регрессии это не так.

Наивная линейная регрессия

Как обычно для линейных моделей, категориальные переменные необходимо преобразовать в one-hot-представление. Для единообразия масштабируем числовые признаки к одному диапазону 0–1 с помощью MinMaxScaler, хотя в данном случае это не сильно влияет на результаты, поскольку признаки уже имеют сопоставимые масштабы:

from sklearn.linear_model import RidgeCV
from sklearn.preprocessing import MinMaxScaler, OneHotEncoder

one_hot_encoder = OneHotEncoder(handle_unknown="ignore", sparse_output=False)
alphas = np.logspace(-6, 6, 25)
naive_linear_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
        ],
        remainder=MinMaxScaler(),
        verbose_feature_names_out=False,
    ),
    RidgeCV(alphas=alphas),
)


evaluate(
    naive_linear_pipeline, X, y, cv=ts_cv, model_prop="alpha_", model_step="ridgecv"
)

Mean model.alpha_ = 2.7298221281347037
Mean Absolute Error:     0.142 +/- 0.014
Root Mean Squared Error: 0.184 +/- 0.020

Pipeline(steps=[('columntransformer',
                 ColumnTransformer(remainder=MinMaxScaler(),
                                   transformers=[('categorical',
                                                  OneHotEncoder(handle_unknown='ignore',
                                                                sparse_output=False),
                                                  Index(['season', 'holiday', 'workingday', 'weather'], dtype='str'))],
                                   verbose_feature_names_out=False)),
                ('ridgecv',
                 RidgeCV(alphas=array([1.00000000e-06, 3.16227766e-06, 1.00000000e-05, 3.16227766e-05,
       1.00000000e-04, 3.16227766e-04, 1.00000000e-03, 3.16227766e-03,
       1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01,
       1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01,
       1.00000000e+02, 3.16227766e+02, 1.00000000e+03, 3.16227766e+03,
       1.00000000e+04, 3.16227766e+04, 1.00000000e+05, 3.16227766e+05,
       1.00000000e+06])))])

Радует, что выбранное alpha_ находится в указанном нами диапазоне.

Результаты неудовлетворительны: средняя ошибка составляет около 14% от максимального спроса. Это более чем втрое превышает среднюю ошибку модели градиентного бустинга. Можно предположить, что наивное исходное кодирование периодических временных признаков (простое масштабирование min-max) мешает модели линейной регрессии правильно использовать временную информацию: линейная регрессия не моделирует автоматически немонотонные зависимости между входными признаками и целевой переменной. Нелинейные члены необходимо создавать для входных данных.

Например, исходное числовое кодирование признака "hour" не позволяет линейной модели распознать, что увеличение часа утром с 6 до 8 должно сильно увеличить количество аренд велосипедов, тогда как увеличение на сопоставимую величину вечером с 18 до 20 должно сильно уменьшить прогнозируемое количество аренд велосипедов.

Временные шаги как категории

Поскольку временные признаки закодированы дискретным образом с использованием целых чисел (24 уникальных значения в признаке «часы»), мы можем рассматривать их как категориальные переменные, используя кодирование one-hot, и тем самым игнорировать любые предположения, связанные с упорядочиванием значений часа.

Использование кодирования one-hot для временных признаков предоставляет линейной модели гораздо больше возможностей, поскольку мы добавляем по одному дополнительному признаку для каждого дискретного временного уровня.

one_hot_linear_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
            ("one_hot_time", one_hot_encoder, ["hour", "weekday", "month"]),
        ],
        remainder=MinMaxScaler(),
        verbose_feature_names_out=False,
    ),
    RidgeCV(alphas=alphas),
)

evaluate(one_hot_linear_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.099 +/- 0.011
Root Mean Squared Error: 0.131 +/- 0.011

Pipeline(steps=[('columntransformer',
                 ColumnTransformer(remainder=MinMaxScaler(),
                                   transformers=[('categorical',
                                                  OneHotEncoder(handle_unknown='ignore',
                                                                sparse_output=False),
                                                  Index(['season', 'holiday', 'workingday', 'weather'], dtype='str')),
                                                 ('one_hot_time',
                                                  OneHotEncoder(handle_unknown='ignore',
                                                                sparse_output=False),
                                                  ['hour', 'weekday',
                                                   'month'])],
                                   verbose_featur...
                 RidgeCV(alphas=array([1.00000000e-06, 3.16227766e-06, 1.00000000e-05, 3.16227766e-05,
       1.00000000e-04, 3.16227766e-04, 1.00000000e-03, 3.16227766e-03,
       1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01,
       1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01,
       1.00000000e+02, 3.16227766e+02, 1.00000000e+03, 3.16227766e+03,
       1.00000000e+04, 3.16227766e+04, 1.00000000e+05, 3.16227766e+05,
       1.00000000e+06])))])

['year', 'temp', 'feel_temp', 'humidity', 'windspeed']

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

Однако при этом появляется очень большое количество новых признаков. Если бы время суток было представлено в минутах с начала дня, а не в часах, кодирование one-hot создало бы 1440 признаков вместо 24. Это могло бы привести к значительному переобучению. Чтобы избежать этого, можно использовать sklearn.preprocessing.KBinsDiscretizer, чтобы заново разбить на интервалы уровни детализированных порядковых или числовых переменных, сохранив при этом преимущества кодирования one-hot, обеспечивающего немонотонную выразительность.

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

Тригонометрические признаки

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

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

from sklearn.preprocessing import FunctionTransformer


def sin_transformer(period):
    return FunctionTransformer(
        lambda x: np.sin(x / period * 2 * np.pi), feature_names_out="one-to-one"
    )


def cos_transformer(period):
    return FunctionTransformer(
        lambda x: np.cos(x / period * 2 * np.pi), feature_names_out="one-to-one"
    )

Визуализируем влияние такого расширения признаков на синтетических данных о часах, немного экстраполировав их за пределы hour=23:

import pandas as pd

hour_df = pd.DataFrame(
    np.arange(26).reshape(-1, 1),
    columns=["hour"],
)
hour_df["hour_sin"] = sin_transformer(24).fit_transform(hour_df)["hour"]
hour_df["hour_cos"] = cos_transformer(24).fit_transform(hour_df)["hour"]
hour_df.plot(x="hour")
_ = plt.title("Trigonometric encoding for the 'hour' feature")

Используем двумерную диаграмму рассеяния, закодировав часы цветом, чтобы нагляднее увидеть, как это представление отображает 24 часа суток в двумерное пространство, напоминающее своего рода аналоговые часы на 24 часа. Обратите внимание, что «25-й» час отображается как 1-й час из-за периодической природы представления с помощью синуса и косинуса.

fig, ax = plt.subplots(figsize=(7, 5))
sp = ax.scatter(hour_df["hour_sin"], hour_df["hour_cos"], c=hour_df["hour"])
ax.set(
    xlabel="sin(hour)",
    ylabel="cos(hour)",
)
_ = fig.colorbar(sp)

Теперь можно построить конвейер извлечения признаков, используя эту стратегию:

cyclic_cossin_transformer = ColumnTransformer(
    transformers=[
        ("categorical", one_hot_encoder, categorical_columns),
        ("month_sin", sin_transformer(12), ["month"]),
        ("month_cos", cos_transformer(12), ["month"]),
        ("weekday_sin", sin_transformer(7), ["weekday"]),
        ("weekday_cos", cos_transformer(7), ["weekday"]),
        ("hour_sin", sin_transformer(24), ["hour"]),
        ("hour_cos", cos_transformer(24), ["hour"]),
    ],
    remainder=MinMaxScaler(),
    verbose_feature_names_out=True,
)
cyclic_cossin_linear_pipeline = make_pipeline(
    cyclic_cossin_transformer,
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_cossin_linear_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.125 +/- 0.014
Root Mean Squared Error: 0.166 +/- 0.020

Pipeline(steps=[('columntransformer',
                 ColumnTransformer(remainder=MinMaxScaler(),
                                   transformers=[('categorical',
                                                  OneHotEncoder(handle_unknown='ignore',
                                                                sparse_output=False),
                                                  Index(['season', 'holiday', 'workingday', 'weather'], dtype='str')),
                                                 ('month_sin',
                                                  FunctionTransformer(feature_names_out='one-to-one',
                                                                      func=<function sin_transformer.<locals>.<lambda> at 0x784...
                 RidgeCV(alphas=array([1.00000000e-06, 3.16227766e-06, 1.00000000e-05, 3.16227766e-05,
       1.00000000e-04, 3.16227766e-04, 1.00000000e-03, 3.16227766e-03,
       1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01,
       1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01,
       1.00000000e+02, 3.16227766e+02, 1.00000000e+03, 3.16227766e+03,
       1.00000000e+04, 3.16227766e+04, 1.00000000e+05, 3.16227766e+05,
       1.00000000e+06])))])

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

Периодические сплайновые признаки

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

from sklearn.preprocessing import SplineTransformer


def periodic_spline_transformer(period, n_splines=None, degree=3):
    if n_splines is None:
        n_splines = period
    n_knots = n_splines + 1  # periodic and include_bias is True
    return SplineTransformer(
        degree=degree,
        n_knots=n_knots,
        knots=np.linspace(0, period, n_knots).reshape(n_knots, 1),
        extrapolation="periodic",
        include_bias=True,
    )

Теперь визуализируем влияние такого расширения признаков на синтетические данные о часах, выполнив небольшую экстраполяцию за пределы hour=23:

hour_df = pd.DataFrame(
    np.linspace(0, 26, 1000).reshape(-1, 1),
    columns=["hour"],
)
splines = periodic_spline_transformer(24, n_splines=12).fit_transform(hour_df)
splines_df = pd.DataFrame(
    splines,
    columns=[f"spline_{i}" for i in range(splines.shape[1])],
)
pd.concat([hour_df, splines_df], axis="columns").plot(x="hour", cmap=plt.cm.tab20b)
_ = plt.title("Periodic spline-based encoding for the 'hour' feature")

Благодаря использованию параметра extrapolation="periodic" мы видим, что кодирование признаков остается гладким при экстраполяции за полночь.

Теперь можно построить прогнозирующий конвейер, используя эту альтернативную стратегию создания периодических признаков.

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

cyclic_spline_transformer = ColumnTransformer(
    transformers=[
        ("categorical", one_hot_encoder, categorical_columns),
        ("cyclic_month", periodic_spline_transformer(12, n_splines=6), ["month"]),
        ("cyclic_weekday", periodic_spline_transformer(7, n_splines=3), ["weekday"]),
        ("cyclic_hour", periodic_spline_transformer(24, n_splines=12), ["hour"]),
    ],
    remainder=MinMaxScaler(),
    verbose_feature_names_out=False,
)
cyclic_spline_linear_pipeline = make_pipeline(
    cyclic_spline_transformer,
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_linear_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.097 +/- 0.011
Root Mean Squared Error: 0.132 +/- 0.013

Pipeline(steps=[('columntransformer',
                 ColumnTransformer(remainder=MinMaxScaler(),
                                   transformers=[('categorical',
                                                  OneHotEncoder(handle_unknown='ignore',
                                                                sparse_output=False),
                                                  Index(['season', 'holiday', 'workingday', 'weather'], dtype='str')),
                                                 ('cyclic_month',
                                                  SplineTransformer(extrapolation='periodic',
                                                                    knots=array([[ 0.],
       [ 2.],
       [ 4.],
       [ 6.],
       [ 8.],
       [10.],
       [12.]]),
                                                                    n_knots...
                 RidgeCV(alphas=array([1.00000000e-06, 3.16227766e-06, 1.00000000e-05, 3.16227766e-05,
       1.00000000e-04, 3.16227766e-04, 1.00000000e-03, 3.16227766e-03,
       1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01,
       1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01,
       1.00000000e+02, 3.16227766e+02, 1.00000000e+03, 3.16227766e+03,
       1.00000000e+04, 3.16227766e+04, 1.00000000e+05, 3.16227766e+05,
       1.00000000e+06])))])

Сплайновые признаки позволяют линейной модели эффективно использовать периодические признаки, связанные со временем, и уменьшить ошибку примерно с 14% до 10% от максимального спроса. Это сопоставимо с результатом, который мы получили для признаков, закодированных методом one-hot.

Качественный анализ влияния признаков на предсказания линейной модели

Здесь мы хотим визуализировать влияние выбранных способов конструирования признаков на форму предсказаний, зависящую от времени.

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

naive_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
naive_linear_predictions = naive_linear_pipeline.predict(X.iloc[test_0])

one_hot_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
one_hot_linear_predictions = one_hot_linear_pipeline.predict(X.iloc[test_0])

cyclic_cossin_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_cossin_linear_predictions = cyclic_cossin_linear_pipeline.predict(X.iloc[test_0])

cyclic_spline_linear_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_spline_linear_predictions = cyclic_spline_linear_pipeline.predict(X.iloc[test_0])

Чтобы получить качественное представление о результатах, рассмотрим предсказания за последние 96 часов (4 дня) тестовой выборки:

last_hours = slice(-96, None)
fig, ax = plt.subplots(figsize=(12, 4))
fig.suptitle("Predictions by linear models")
ax.plot(
    y.iloc[test_0].values[last_hours],
    "x-",
    alpha=0.2,
    label="Actual demand",
    color="black",
)
ax.plot(naive_linear_predictions[last_hours], "x-", label="Ordinal time features")
ax.plot(
    cyclic_cossin_linear_predictions[last_hours],
    "x-",
    label="Trigonometric time features",
)
ax.plot(
    cyclic_spline_linear_predictions[last_hours],
    "x-",
    label="Spline-based time features",
)
ax.plot(
    one_hot_linear_predictions[last_hours],
    "x-",
    label="One-hot time features",
)
_ = ax.legend()

По приведённому выше графику можно сделать следующие выводы:

• Исходные порядковые признаки, связанные со временем создают проблемы, поскольку не учитывают естественную периодичность: в конце каждого дня мы наблюдаем резкий скачок предсказаний, когда значение признака часа меняется с 23 на 0. Подобные артефакты можно ожидать и в конце каждой недели или каждого года.
• Как и ожидалось, тригонометрические признаки (синус и косинус) не создают разрывов в полночь, однако модели линейной регрессии не удаётся эффективно использовать эти признаки для правильного моделирования внутридневных изменений. Эту проблему можно потенциально решить, используя тригонометрические признаки для высших гармоник или дополнительные тригонометрические признаки для естественного периода с разными фазами.
• Признаки на основе периодических сплайнов одновременно решают обе эти проблемы: они повышают выразительность линейной модели, позволяя ей сосредоточиться на определённых часах благодаря использованию 12 сплайнов. Кроме того, параметр extrapolation="periodic" обеспечивает плавное представление между hour=23 и hour=0.
• Признаки, закодированные методом one-hot, ведут себя подобно признакам на основе периодических сплайнов, но дают более резкие изменения: например, они лучше моделируют утренний пик в будние дни, поскольку этот пик длится меньше часа. Однако далее мы увидим, что то, что может быть преимуществом для линейных моделей, не обязательно является таковым для более выразительных моделей.

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

naive_linear_pipeline[:-1].transform(X).shape

(17379, 19)

one_hot_linear_pipeline[:-1].transform(X).shape

(17379, 59)

cyclic_cossin_linear_pipeline[:-1].transform(X).shape

(17379, 22)

cyclic_spline_linear_pipeline[:-1].transform(X).shape

(17379, 37)

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

Наконец, мы видим, что ни одна из линейных моделей не может аппроксимировать фактический спрос на аренду велосипедов, особенно в периоды пиков, которые могут быть очень резкими в часы пик в рабочие дни, но гораздо более сглаженными в выходные: наиболее точные линейные модели на основе сплайнов или кодирования one-hot склонны прогнозировать пики поездок на работу на велосипеде даже в выходные и недооценивать связанные с поездками на работу события в будние дни.

Эти систематические ошибки предсказания указывают на недообучение и могут объясняться отсутствием взаимодействий между признаками, например между «workingday» и признаками, производными от «hours». Эта проблема будет рассмотрена в следующем разделе.

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

Линейные модели автоматически не учитывают эффекты взаимодействия между входными признаками. Не помогает и то, что некоторые признаки нелинейны по отдельности, как в случае признаков, построенных с помощью SplineTransformer (или кодирования one-hot, или разбиения на интервалы).

Однако можно использовать класс PolynomialFeatures для часов, закодированных сплайнами с крупным шагом, чтобы явно моделировать взаимодействие «workingday»/»hours», не вводя слишком много новых переменных:

from sklearn.pipeline import FeatureUnion
from sklearn.preprocessing import PolynomialFeatures

hour_workday_interaction = make_pipeline(
    ColumnTransformer(
        [
            ("cyclic_hour", periodic_spline_transformer(24, n_splines=8), ["hour"]),
            (
                "workingday",
                FunctionTransformer(
                    lambda x: x == "True", feature_names_out="one-to-one"
                ),
                ["workingday"],
            ),
        ],
        verbose_feature_names_out=False,
    ),
    PolynomialFeatures(degree=2, interaction_only=True, include_bias=False),
)

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

cyclic_spline_interactions_pipeline = make_pipeline(
    FeatureUnion(
        [
            ("marginal", cyclic_spline_transformer),
            ("interactions", hour_workday_interaction),
        ],
        verbose_feature_names_out=True,
    ).set_output(transform="pandas"),
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_interactions_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.078 +/- 0.009
Root Mean Squared Error: 0.104 +/- 0.009

Pipeline(steps=[('featureunion',
                 FeatureUnion(transformer_list=[('marginal',
                                                 ColumnTransformer(remainder=MinMaxScaler(),
                                                                   transformers=[('categorical',
                                                                                  OneHotEncoder(handle_unknown='ignore',
                                                                                                sparse_output=False),
                                                                                  Index(['season', 'holiday', 'workingday', 'weather'], dtype='str')),
                                                                                 ('cyclic_month',
                                                                                  SplineTransformer(extrapolation='periodic',
                                                                                                    knots=array([[ 0.],
       [ 2....
                 RidgeCV(alphas=array([1.00000000e-06, 3.16227766e-06, 1.00000000e-05, 3.16227766e-05,
       1.00000000e-04, 3.16227766e-04, 1.00000000e-03, 3.16227766e-03,
       1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01,
       1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01,
       1.00000000e+02, 3.16227766e+02, 1.00000000e+03, 3.16227766e+03,
       1.00000000e+04, 3.16227766e+04, 1.00000000e+05, 3.16227766e+05,
       1.00000000e+06])))])