Практическое руководство и платформа для использования любой модели scikit-learn для прогнозирования временных рядов в Python.
Введение
Существует много так называемых традиционных моделей для прогнозирования временных рядов, таких как семейство моделей SARIMAX, экспоненциальное сглаживание или BATS и TBATS.
Однако очень редко мы упоминаем наиболее распространенные модели машинного обучения для регрессии, такие как деревья решений, случайные леса, повышение градиента или даже регрессор опорных векторов. Мы видим, что эти модели широко применяются в типичных задачах регрессии, но не для прогнозирования временных рядов.
Отсюда и причина написания этой статьи! Здесь мы разрабатываем структуру, чтобы сформулировать проблему временных рядов как задачу контролируемого обучения, что позволяет нам использовать любую модель из нашей любимой библиотеки: scikit-learn!
К концу этой статьи у вас будут инструменты и знания для применения любой модели машинного обучения для прогнозирования временных рядов вместе со статистическими моделями, упомянутыми выше.
Давайте начнем!
Полный исходный код доступен на GitHub.
Изучите новейшие методы анализа временных рядов с помощью моей бесплатной шпаргалка по временным рядам на Python! Получите реализацию статистических методов и методов глубокого обучения на Python и TensorFlow!
Подготовка набора данных
Во-первых, мы импортируем все библиотеки, необходимые для завершения нашего руководства.
import numpy as np import pandas as pd import statsmodels.api as sm import matplotlib.pyplot as plt
Здесь мы используем библиотеку statsmodels для импорта набора данных, который представляет собой еженедельную концентрацию CO2 с 1958 по 2001 год.
data = sm.datasets.co2.load_pandas().data
Хороший первый шаг — визуализировать наши данные с помощью следующего блока кода.
fig, ax = plt.subplots(figsize=(16, 11))
ax.plot(data['co2'])
ax.set_xlabel('Time')
ax.set_ylabel('CO2 concentration (ppmw)')
fig.autofmt_xdate()
plt.tight_layout()
На изображении выше мы видим четкую положительную тенденцию в данных, так как концентрация со временем увеличивается. Мы также наблюдаем годовой сезонный характер. Это связано со сменой сезонов, когда зимой концентрация CO2 выше, чем летом. Наконец, мы видим некоторые недостающие данные в начале набора данных.
Давайте затем обработаем эти недостающие данные с помощью интерполяции. Мы просто выполним линейную интерполяцию между двумя известными точками, чтобы заполнить недостающие значения.
data = data.interpolate()
Теперь, когда у нас нет недостающих данных, мы готовы приступить к моделированию!
Моделирование с помощью scikit-learn
Как вы увидите, самая большая проблема в прогнозировании временных рядов с помощью scikit-learn заключается в правильной постановке задачи.
Есть 3 различных способа, которыми мы можем сформулировать задачу прогнозирования временных рядов как задачу обучения с учителем:
- Предсказать следующий временной шаг, используя предыдущее наблюдение
- Предсказать следующий временной шаг, используя последовательность прошлых наблюдений
- Предсказать последовательность будущих временных шагов, используя последовательность прошлых наблюдений
Разберем каждую ситуацию подробнее!
Предсказать следующий временной шаг, используя предыдущее наблюдение
Это самая основная установка. Модель выводит прогноз для следующего временного шага, учитывая только предыдущее наблюдение, как показано на рисунке ниже.
Это простой вариант использования с небольшим количеством практических применений, поскольку модель, скорее всего, не собирается учиться чему-либо только на основе предыдущего наблюдения. Тем не менее, это служит хорошей отправной точкой, чтобы помочь нам понять более сложные сценарии позже.
Итак, на данный момент наш набор данных выглядит так:
Что не очень полезно. Правильно, все значения находятся в одном столбце, но нам нужно отформатировать набор данных таким образом, чтобы текущее наблюдение было функцией для прогнозирования следующего наблюдения (цели).
Таким образом, мы добавляем второй столбец, который просто сдвигает столбец co2 таким образом, что значение в 1958–03–29 теперь является предиктором значения в 1958–04–05.
df = data.copy() df['y'] = df['co2'].shift(-1)
Как вы можете, наш набор данных теперь отформатирован таким образом, что каждое текущее наблюдение является предиктором для следующего наблюдения! Обратите внимание, что у нас отсутствует значение конца нашего набора данных. Это нормально для последнего известного наблюдения. Мы удалим эту строку на следующем шаге.
А пока давайте разделим наш набор данных на набор для обучения и тестирования, чтобы запустить наши модели и оценить их. Здесь мы используем данные за последние два года для обучающей выборки. Поскольку у нас еженедельные данные, а в году 52 недели, это означает, что последние 104 выборки сохраняются для тестового набора.
train = df[:-104] test = df[-104:] test = test.drop(test.tail(1).index) # Drop last row
Базовая модель
Конечно, нам нужна базовая модель, чтобы определить, лучше ли использовать модели машинного обучения. Здесь мы наивно предсказываем, что следующее наблюдение будет иметь то же значение, что и текущее наблюдение.
Другими словами, мы просто устанавливаем столбец co2 в качестве наших базовых прогнозов.
test = test.copy() test['baseline_pred'] = test['co2']
Большой! Сделав этот шаг, давайте перейдем к более сложным моделям.
Дерево решений
Здесь давайте применим регрессор дерева решений. Эту модель можно заменить любой моделью из библиотеки scikit-learn! Обратите внимание, что мы используем случайное состояние для обеспечения воспроизводимости.
from sklearn.tree import DecisionTreeRegressor X_train = train['co2'].values.reshape(-1,1) y_train = train['y'].values.reshape(-1,1) X_test = test['co2'].values.reshape(-1,1) # Initialize the model dt_reg = DecisionTreeRegressor(random_state=42) # Fit the model dt_reg.fit(X=X_train, y=y_train) # Make predictions dt_pred = dt_reg.predict(X_test) # Assign predictions to a new column in test test['dt_pred'] = dt_pred
Обратите внимание, что мы сохраняем часть моделирования простой. Мы не выполняем никакой перекрестной проверки или настройки гиперпараметров, хотя эти методы могут применяться здесь обычно, как и в любой другой задаче регрессии.
Не стесняйтесь применять эти методы и посмотрите, сможете ли вы улучшить производительность.
Усиление градиента
Чтобы попробовать разные модели, давайте теперь применим усиление градиента.
from sklearn.ensemble import GradientBoostingRegressor gbr = GradientBoostingRegressor(random_state=42) gbr.fit(X_train, y=y_train.ravel()) gbr_pred = gbr.predict(X_test) test['gbr_pred'] = gbr_pred
Большой! Теперь у нас есть прогнозы на основе двух моделей машинного обучения и базового уровня. Пришло время оценить эффективность каждого метода.
Оценка
Здесь мы используем среднюю абсолютную процентную ошибку (MAPE). Это особенно информативная метрика ошибок, так как она возвращает процент, который легко интерпретировать. Обязательно применяйте его только тогда, когда у вас нет значений, близких к 0, как в данном случае.
К сожалению, MAPE еще не реализован в scikit-learn, поэтому мы должны определить функцию вручную.
def mape(y_true, y_pred):
return round(np.mean(np.abs((y_true - y_pred) / y_true)) * 100, 2)
Затем мы можем оценить каждую модель и создать следующую гистограмму:
baseline_mape = mape(test['y'], test['baseline_pred'])
dt_mape = mape(test['y'], test['dt_pred'])
gbr_mape = mape(test['co2'], test['gbr_pred'])
# Generate bar plot
fig, ax = plt.subplots(figsize=(7, 5))
x = ['Baseline', 'Decision Tree', 'Gradient Boosting']
y = [baseline_mape, dt_mape, gbr_mape]
ax.bar(x, y, width=0.4)
ax.set_xlabel('Regressor models')
ax.set_ylabel('MAPE (%)')
ax.set_ylim(0, 0.3)
for index, value in enumerate(y):
plt.text(x=index, y=value + 0.02, s=str(value), ha='center')
plt.tight_layout()
Глядя на рисунок выше, мы видим, что базовый вариант имеет наилучшую производительность, так как у него самый низкий показатель MAPE. Это имеет смысл, так как концентрация CO2, кажется, не меняется резко от одной недели к другой.
В данном случае использование моделей машинного обучения не дало нам никакой дополнительной ценности. Опять же, это может быть связано с тем, что модель учится только на одном наблюдении, чтобы сделать прогноз. Возможно, было бы лучше дать ему последовательность в качестве входных данных, чтобы предсказать следующий временной шаг.
Это подводит нас к следующему сценарию!
Предсказать следующий временной шаг, используя последовательность прошлых наблюдений
Как мы видели в предыдущем примере, использование одного наблюдения для предсказания следующего временного шага не очень хорошо работает. Теперь давайте попробуем использовать последовательность в качестве входных данных для модели и предсказать следующий временной шаг, как показано ниже.
Опять же, мы должны отформатировать наш набор данных таким образом, чтобы у нас была последовательность прошлых наблюдений, действующих как предикторы для следующего временного шага.
Мы можем легко написать функцию, которая добавляет сдвинутые столбцы, чтобы получить желаемую длину ввода.
def window_input(window_length: int, data: pd.DataFrame) -> pd.DataFrame:
df = data.copy()
i = 1
while i < window_length:
df[f'x_{i}'] = df['co2'].shift(-i)
i = i + 1
if i == window_length:
df['y'] = df['co2'].shift(-i)
# Drop rows where there is a NaN
df = df.dropna(axis=0)
return df
Давайте используем эту функцию для ввода 5 наблюдений, чтобы предсказать следующий временной шаг.
new_df = window_input(5, data)
Глядя на рисунок выше, мы видим, что наш набор данных устроен таким образом, что у нас есть пять наблюдений для прогнозирования следующего временного шага, хранящихся в столбце y.
Это по существу заботится о самой сложной части! Теперь нужно просто применить разные модели и посмотреть, какая из них работает лучше всего.
Прежде чем перейти к этому шагу, давайте сначала разделим наши данные на обучающий и тестовый наборы.
from sklearn.model_selection import train_test_split X = new_df[['co2', 'x_1', 'x_2', 'x_3', 'x_4']].values y = new_df['y'].values X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42, shuffle=False)
Базовая модель
Опять же, давайте определим базовую модель для этой ситуации. Здесь мы просто предскажем среднее значение входной последовательности
baseline_pred = []
for row in X_test:
baseline_pred.append(np.mean(row))
Дерево решений
Опять же, давайте применим регрессор дерева решений. Это так же просто, как и предыдущая реализация.
dt_reg_5 = DecisionTreeRegressor(random_state=42) dt_reg_5.fit(X_train, y_train) dt_reg_5_pred = dt_reg_5.predict(X_test)
Усиление градиента
Ради последовательности, давайте также попробуем повышение градиента.
gbr_5 = GradientBoostingRegressor(random_state=42) gbr_5.fit(X_train, y_train.ravel()) gbr_5_pred = gbr_5.predict(X_test)
Оценка
Теперь мы можем оценить производительность каждой модели. Опять же, мы используем MAPE и отображаем результаты в виде гистограммы. По сути, это тот же код, что и раньше.
baseline_mape = mape(y_test, baseline_pred)
dt_5_mape = mape(y_test, dt_reg_5_pred)
gbr_5_mape = mape(y_test, gbr_5_pred)
# Generate the bar plot
fig, ax = plt.subplots()
x = ['Baseline', 'Decision Tree', 'Gradient Boosting']
y = [baseline_mape, dt_5_mape, gbr_5_mape]
ax.bar(x, y, width=0.4)
ax.set_xlabel('Regressor models')
ax.set_ylabel('MAPE (%)')
ax.set_ylim(0, 2.5)
for index, value in enumerate(y):
plt.text(x=index, y=value + 0.1, s=str(value), ha='center')
plt.tight_layout()
Удивительно, но на рисунке выше мы видим, что модели машинного обучения не превосходят базовый уровень.
Хотя настройка гиперпараметров, вероятно, может повысить производительность моделей машинного обучения, я подозреваю, что окно ввода слишком короткое. Ввода в пять недель недостаточно, чтобы модель уловила тренд и сезонные компоненты, поэтому нам может потребоваться более длительный временной интервал.
Мы попробуем именно это в следующем сценарии!
Предсказать последовательность будущих временных шагов, используя последовательность прошлых наблюдений
Последний сценарий использует последовательность наблюдений для прогнозирования последовательности будущих временных шагов, как показано ниже.
Здесь наша модель требуется для вывода последовательности прогнозов. Это можно рассматривать как проблему регрессии с несколькими выходами.
Итак, первый шаг — правильно отформатировать наш набор данных. Мы разрабатываем другую функцию, которая использует метод shift для форматирования набора данных как задачи регрессии с несколькими выходами.
def window_input_output(input_length: int, output_length: int, data: pd.DataFrame) -> pd.DataFrame:
df = data.copy()
i = 1
while i < input_length:
df[f'x_{i}'] = df['co2'].shift(-i)
i = i + 1
j = 0
while j < output_length:
df[f'y_{j}'] = df['co2'].shift(-output_length-j)
j = j + 1
df = df.dropna(axis=0)
return df
Здесь мы будем использовать последовательность из 26 наблюдений, чтобы предсказать следующие 26 временных шагов. Другими словами, мы вводим полгода, чтобы предсказать следующую половину. Обратите внимание, что входная и выходная последовательности не обязательно должны иметь одинаковую длину. Это произвольное решение с моей стороны.
seq_df = window_input_output(26, 26, data)
Как видите, теперь у нас есть набор данных, в котором 26 наблюдений используются в качестве предикторов для следующих 26 временных шагов.
Прежде чем перейти к моделированию, мы снова разделим данные на обучающий и тестовый наборы. Здесь мы резервируем последние две строки для тестового набора, так как это дает нам 52 тестовых образца.
X_cols = [col for col in seq_df.columns if col.startswith('x')]
X_cols.insert(0, 'co2')
y_cols = [col for col in seq_df.columns if col.startswith('y')]
X_train = seq_df[X_cols][:-2].values
y_train = seq_df[y_cols][:-2].values
X_test = seq_df[X_cols][-2:].values
y_test = seq_df[y_cols][-2:].values
Базовая модель
В качестве базовой модели мы просто повторим входную последовательность. Другими словами, мы берем входную последовательность и выводим ту же последовательность, что и наши базовые прогнозы.
Это очень тривиальное предсказание, поэтому мы реализуем его, когда будем готовы оценивать модели.
Дерево решений
Опять же, давайте попробуем применить дерево решений. Обратите внимание, что дерево решений может давать прогнозы с несколькими выходами, поэтому здесь нам не нужно выполнять какую-либо дополнительную работу.
dt_seq = DecisionTreeRegressor(random_state=42) dt_seq.fit(X_train, y_train) dt_seq_preds = dt_seq.predict(X_test)
Усиление градиента
Теперь повышение градиента требует немного дополнительной работы. Он не может обрабатывать цель с несколькими выходами. Попытка немедленно подобрать модель повышения градиента приведет к ошибке.
Здесь мы должны обернуть модель таким образом, чтобы ее прогноз использовался в качестве входных данных для подачи следующего прогноза. Это достигается с помощью оболочки RegressorChain из scikit-learn.
from sklearn.multioutput import RegressorChain gbr_seq = GradientBoostingRegressor(random_state=42) chained_gbr = RegressorChain(gbr_seq) chained_gbr.fit(X_train, y_train) gbr_seq_preds = chained_gbr.predict(X_test)
Это позволило обучить модель и делать прогнозы без каких-либо ошибок. Как я уже упоминал, за кулисами модель предсказывает следующий временной шаг и использует этот прогноз для создания следующего прогноза. Недостаток этого метода в том, что если первый прогноз плохой, то и остальная часть последовательности, скорее всего, будет плохой.
Теперь оценим каждую модель.
Оценка
Опять же, мы используем MAPE для оценки наших методов прогнозирования.
mape_dt_seq = mape(dt_seq_preds.reshape(1, -1), y_test.reshape(1, -1))
mape_gbr_seq = mape(gbr_seq_preds.reshape(1, -1), y_test.reshape(1, -1))
mape_baseline = mape(X_test.reshape(1, -1), y_test.reshape(1, -1))
# Generate the bar plot
fig, ax = plt.subplots()
x = ['Baseline', 'Decision Tree', 'Gradient Boosting']
y = [mape_baseline, mape_dt_seq, mape_gbr_seq]
ax.bar(x, y, width=0.4)
ax.set_xlabel('Regressor models')
ax.set_ylabel('MAPE (%)')
ax.set_ylim(0, 1)
for index, value in enumerate(y):
plt.text(x=index, y=value + 0.05, s=str(value), ha='center')
plt.tight_layout()
На рисунке выше видно, что нам наконец удалось обучить модели машинного обучения, которые превосходят базовые! Здесь модель дерева решений является лучшей моделью, поскольку она обеспечивает самый низкий показатель MAPE.
Опционально визуализируем прогнозы за последний год.
fig, ax = plt.subplots(figsize=(16, 11))
ax.plot(np.arange(0, 26, 1), X_test[1], 'b-', label='input')
ax.plot(np.arange(26, 52, 1), y_test[1], marker='.', color='blue', label='Actual')
ax.plot(np.arange(26, 52, 1), X_test[1], marker='o', color='red', label='Baseline')
ax.plot(np.arange(26, 52, 1), dt_seq_preds[1], marker='^', color='green', label='Decision Tree')
ax.plot(np.arange(26, 52, 1), gbr_seq_preds[1], marker='P', color='black', label='Gradient Boosting')
ax.set_xlabel('Timesteps')
ax.set_ylabel('CO2 concentration (ppmv)')
plt.xticks(np.arange(1, 104, 52), np.arange(2000, 2002, 1))
plt.legend(loc=2)
fig.autofmt_xdate()
plt.tight_layout()
Кроме того, поддерживая MAPE, мы видим, что дерево решений и повышение градиента более точно соответствуют фактическим значениям, чем базовые прогнозы.
В этой статье мы увидели, как сформулировать задачу прогнозирования временных рядов как проблему регрессии, которую можно решить с помощью регрессионных моделей scikit-learn.
Мы исследовали следующие сценарии:
- Предсказать следующий временной шаг, используя предыдущее наблюдение
- Предсказать следующий временной шаг, используя последовательность прошлых наблюдений
- Предсказать последовательность будущих временных шагов, используя последовательность прошлых наблюдений
Теперь у нас есть структура, позволяющая представить любую задачу прогнозирования временных рядов как задачу контролируемого обучения, к которой можно применить любую модель регрессора из scikit-learn.
Я надеюсь, что вы нашли эту статью полезной!
Обязательно загрузите мою бесплатную шпаргалка по прогнозированию временных рядов на Python, охватывающую как статистические модели, так и модели глубокого обучения!
Ваше здоровье! 🍺
