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

Временная перекрестная проверка

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

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", чтобы признаки с категориальным типом данных рассматривались как категориальные. Для справки, мы извлекаем категориальные признаки из фрейма данных на основе типа данных. Внутренние деревья используют специальное правило разделения деревьев для этих признаков.

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

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=model_prop is not None,
    )
    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}"
    )


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

Мы видим, что мы установили max_iter достаточно большой, чтобы произошла остановка по раннему прекращению.

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

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

Простая линейная регрессия

Как обычно для линейных моделей, категориальные переменные необходимо закодировать с помощью метода one-hot encoding. Для согласованности мы масштабируем числовые признаки в тот же диапазон от 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(),
    ),
    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

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

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

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

Шаги времени как категории

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

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

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(),
    ),
    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

Средняя ошибка этой модели составляет 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))


def cos_transformer(period):
    return FunctionTransformer(lambda x: np.cos(x / period * 2 * np.pi))

Давайте визуализируем эффект расширения признаков на некоторых синтетических данных часов с небольшой экстраполяцией за пределы часа=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-й» час отображается обратно в первый час из-за периодической природы синусоидального/косинусного представления.

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(),
)
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

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

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

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

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,
    )

Давайте снова визуализируем эффект расширения признаков на некоторых синтетических данных часов с небольшой экстраполяцией за пределы часа=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(),
)
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

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

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

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

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

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 encoding, ведут себя аналогично периодическим сплайновым признакам, но имеют более резкие пики: например, они могут лучше моделировать утренний пик в будние дни, поскольку этот пик длится меньше часа. Однако в дальнейшем мы увидим, что то, что может быть преимуществом для линейных моделей, не обязательно является преимуществом для более выразительных моделей.

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

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 encoding и сплайнового кодирования создают больше признаков для временного представления, чем альтернативные методы, что, в свою очередь, даёт последующей модели линейной регрессии большую гибкость (степени свободы), чтобы избежать недообучения.

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

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

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

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

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

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"), ["workingday"]),
        ]
    ),
    PolynomialFeatures(degree=2, interaction_only=True, include_bias=False),
)

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

cyclic_spline_interactions_pipeline = make_pipeline(
    FeatureUnion(
        [
            ("marginal", cyclic_spline_transformer),
            ("interactions", hour_workday_interaction),
        ]
    ),
    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

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

Предыдущий анализ выявил необходимость моделирования взаимодействий между "workingday" и "hours". Другим примером такого нелинейного взаимодействия, которое мы хотели бы смоделировать, может быть влияние дождя, которое может быть неодинаковым в рабочие дни и в выходные и праздничные дни, например.

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

В качестве альтернативы мы можем использовать метод Nyström для вычисления приближенного полиномиального ядра. Давайте попробуем последний:

from sklearn.kernel_approximation import Nystroem

cyclic_spline_poly_pipeline = make_pipeline(
    cyclic_spline_transformer,
    Nystroem(kernel="poly", degree=2, n_components=300, random_state=0),
    RidgeCV(alphas=alphas),
)
evaluate(cyclic_spline_poly_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.053 +/- 0.002
Root Mean Squared Error: 0.076 +/- 0.004

Мы наблюдаем, что эта модель может почти соперничать с производительностью деревьев градиентного бустинга со средней ошибкой около 5% от максимального спроса.

Обратите внимание, что, хотя последним шагом этого конвейера является линейная регрессионная модель, промежуточные шаги, такие как извлечение признаков сплайнов и приближение ядра Nyström, являются высоко нелинейными. В результате составной конвейер является гораздо более выразительным, чем простая линейная регрессионная модель с исходными признаками.

Для полноты мы также оцениваем комбинацию one-hot кодирования и приближения ядра:

one_hot_poly_pipeline = make_pipeline(
    ColumnTransformer(
        transformers=[
            ("categorical", one_hot_encoder, categorical_columns),
            ("one_hot_time", one_hot_encoder, ["hour", "weekday", "month"]),
        ],
        remainder="passthrough",
    ),
    Nystroem(kernel="poly", degree=2, n_components=300, random_state=0),
    RidgeCV(alphas=alphas),
)
evaluate(one_hot_poly_pipeline, X, y, cv=ts_cv)

Mean Absolute Error:     0.082 +/- 0.006
Root Mean Squared Error: 0.111 +/- 0.011

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

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

gbrt.fit(X.iloc[train_0], y.iloc[train_0])
gbrt_predictions = gbrt.predict(X.iloc[test_0])

one_hot_poly_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
one_hot_poly_predictions = one_hot_poly_pipeline.predict(X.iloc[test_0])

cyclic_spline_poly_pipeline.fit(X.iloc[train_0], y.iloc[train_0])
cyclic_spline_poly_predictions = cyclic_spline_poly_pipeline.predict(X.iloc[test_0])

Опять же, мы сосредоточимся на последних 4 днях тестового набора:

last_hours = slice(-96, None)
fig, ax = plt.subplots(figsize=(12, 4))
fig.suptitle("Predictions by non-linear regression models")
ax.plot(
    y.iloc[test_0].values[last_hours],
    "x-",
    alpha=0.2,
    label="Actual demand",
    color="black",
)
ax.plot(
    gbrt_predictions[last_hours],
    "x-",
    label="Gradient Boosted Trees",
)
ax.plot(
    one_hot_poly_predictions[last_hours],
    "x-",
    label="One-hot + polynomial kernel",
)
ax.plot(
    cyclic_spline_poly_predictions[last_hours],
    "x-",
    label="Splines + polynomial kernel",
)
_ = ax.legend()

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

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

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

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

Наконец, давайте более количественно оценим ошибки предсказания этих трех моделей с помощью диаграмм рассеяния «фактический против предсказанного» спроса:

from sklearn.metrics import PredictionErrorDisplay

fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(13, 7), sharex=True, sharey="row")
fig.suptitle("Non-linear regression models", y=1.0)
predictions = [
    one_hot_poly_predictions,
    cyclic_spline_poly_predictions,
    gbrt_predictions,
]
labels = [
    "One hot +\npolynomial kernel",
    "Splines +\npolynomial kernel",
    "Gradient Boosted\nTrees",
]
plot_kinds = ["actual_vs_predicted", "residual_vs_predicted"]
for axis_idx, kind in enumerate(plot_kinds):
    for ax, pred, label in zip(axes[axis_idx], predictions, labels):
        disp = PredictionErrorDisplay.from_predictions(
            y_true=y.iloc[test_0],
            y_pred=pred,
            kind=kind,
            scatter_kwargs={"alpha": 0.3},
            ax=ax,
        )
        ax.set_xticks(np.linspace(0, 1, num=5))
        if axis_idx == 0:
            ax.set_yticks(np.linspace(0, 1, num=5))
            ax.legend(
                ["Best model", label],
                loc="upper center",
                bbox_to_anchor=(0.5, 1.3),
                ncol=2,
            )
        ax.set_aspect("equal", adjustable="box")
plt.show()

Это визуализация подтверждает выводы, сделанные на предыдущем графике.

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

Заключительные замечания

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

Регрессор Nystroem + RidgeCV также мог бы быть заменён на MLPRegressor с одним или двумя скрытыми слоями, и мы получили бы довольно похожие результаты.

Используемый в данном случае набор данных сэмплируется по часам. Однако циклические сплайн-признаки могли бы очень эффективно моделировать время в течение дня или недели с более мелким временным разрешением (например, с измерениями каждые минуту вместо каждого часа), не добавляя больше признаков. One-hot кодирование временных представлений не предлагало бы такой гибкости.

Наконец, в этом блокноте мы использовали RidgeCV, поскольку он очень эффективен с точки зрения вычислений. Однако он моделирует целевую переменную как гауссовскую случайную переменную с постоянной дисперсией. Для задач регрессии с положительными значениями, вероятно, было бы более уместно использовать распределение Пуассона или Гамма. Это можно было бы сделать, используя GridSearchCV(TweedieRegressor(power=2), param_grid({"alpha": alphas})) вместо RidgeCV.

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

Download Jupyter notebook: plot_cyclical_feature_engineering.ipynb

Download Python source code: plot_cyclical_feature_engineering.py

Download zipped: plot_cyclical_feature_engineering.zip

Примеры по теме

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

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

Интерполяция многочленами и сплайнами

Графики частичной зависимости и индивидуального условного ожидания