Гибридные модели можно разделить на два вида : первый - как комбинация прогнозов, полученных разными моделями, второй - на разложении временного ряда на составляющие, каждая из которых моделируется своей моделью. В курсе Kaggle рассматривается модель второго вида и именно ее мы будем рассматривать.
Временной ряд можно представить как комбинацию трех компонент : тренд, сезонность и случайная компонента. В результате получаем две модели временного ряда :
аддитивная модель Yt = Tt + St +et.
мультипликативная Yt = Tt*St*et.
Случайную составляющую опустим и займемся трендом и сезонностью. Алгоритм у нас будет такой :
1. Первой моделью сделаем подгонку тренда, для этого используем линейную регрессию
2. Далее в случае аддитивной модели вычтем из ряда прогноз тренда и полученные остатки запустим во вторую модель для подгонки остатков, а в случае мультипликативной модели разделим значения временного ряда на значения тренда и полученные значения также запустим для подгонки во вторую модель.
3. На третьем шаги совместим два прогноза, для аддитивной модели как сумму, для мультипликативной - как произведение.
Загружаем необходимые библиотеки и делаем необходимые настройки :
from warnings import simplefilter
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from xgboost import XGBRegressor
simplefilter("ignore")
plt.style.use("seaborn-whitegrid")
plt.rc(
"figure",
autolayout=True,
figsize=(11, 4),
titlesize=18,
titleweight='bold',
)
plt.rc(
"axes",
labelweight="bold",
labelsize="large",
titleweight="bold",
titlesize=16,
titlepad=10,
)
plot_params = dict(
color="0.75",
style=".-",
markeredgecolor="0.25",
markerfacecolor="0.25",
)
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from xgboost import XGBRegressor
simplefilter("ignore")
plt.style.use("seaborn-whitegrid")
plt.rc(
"figure",
autolayout=True,
figsize=(11, 4),
titlesize=18,
titleweight='bold',
)
plt.rc(
"axes",
labelweight="bold",
labelsize="large",
titleweight="bold",
titlesize=16,
titlepad=10,
)
plot_params = dict(
color="0.75",
style=".-",
markeredgecolor="0.25",
markerfacecolor="0.25",
)
Для примера рассмотрим набор данных, представляющий месячную выручку восьми магазинов с января 2015 по декабрь 2021.
df = pd.read_excel(
"shops8_rev.xlsx",
parse_dates=['month'])
df = df.set_index("month").to_period()
pd.concat([df.head(), df.tail()])
"shops8_rev.xlsx",
parse_dates=['month'])
df = df.set_index("month").to_period()
pd.concat([df.head(), df.tail()])
| month | Sh01 | Sh02 | Sh03 | Sh04 | Sh05 | Sh06 | Sh07 | Sh08 |
|---|---|---|---|---|---|---|---|---|
| 2015-01 | 157408 | 82190 | 174300 | 77120 | 87860 | 51830 | 81790 | 80500 |
| 2015-02 | 127692 | 83010 | 124520 | 83330 | 79210 | 42700 | 62320 | 76530 |
| 2015-03 | 140148 | 91590 | 182450 | 105470 | 98780 | 54490 | 75460 | 110610 |
| 2015-04 | 167924 | 82680 | 183340 | 111520 | 121480 | 53850 | 88780 | 96930 |
| 2015-05 | 184196 | 98930 | 204770 | 123600 | 128200 | 63730 | 92800 | 116680 |
| 2021-08 | 161308 | 48320 | 100380 | 87310 | 73110 | 83550 | 62500 | 39640 |
| 2021-09 | 165948 | 61810 | 136720 | 105370 | 91580 | 109460 | 83650 | 55120 |
| 2021-10 | 216604 | 66430 | 181200 | 105550 | 92360 | 133850 | 85850 | 63220 |
| 2021-11 | 276204 | 70140 | 197760 | 119400 | 98500 | 154980 | 101530 | 78540 |
| 2021-12 | 236616 | 78480 | 191440 | 150360 | 126040 | 183610 | 124600 | 93280 |
Выделим для работы один магазин
df1=df[['Sh01']]
df1.head()
| month | Sh01 |
|---|---|
| 2015-01 | 157408 |
| 2015-02 | 127692 |
| 2015-03 | 140148 |
| 2015-04 | 167924 |
| 2015-05 | 184196 |
Для получения тренда с помощью линейной регрессии необходимо подготовить независимые переменные, для этого используем функцию DeterministicProcess из модуля statsmodels. Задаем наличие константы в регрессионном уравнении, порядок кривой и контроль за коллениарностью.
dp = DeterministicProcess(
index=df1.index, # даты из обучающих данных
constant=True, # фиктивная функция для смещения (y_intercept)
order=2, # квадратичный тренд
drop=True, # если необходимо, отбросить члены, чтобы избежать коллинеарности
)
in_sample` создает функции для дат, указанных в аргументе `index`
X = dp.in_sample()
X.head()
| month | const | trend | trend_squared |
|---|---|---|---|
| 2015-01 | 1.0 | 1.0 | 1.0 |
| 2015-02 | 1.0 | 2.0 | 4.0 |
| 2015-03 | 1.0 | 3.0 | 9.0 |
| 2015-04 | 1.0 | 4.0 | 16.0 |
| 2015-05 | 1.0 | 5.0 | 25.0 |
Разделим индексы дат на обучающие и тестовые
idx_train, idx_test = train_test_split(
y.index, test_size=12, shuffle=False)
y.index, test_size=12, shuffle=False)
С их помощью разделим данные на обучающую и тестовую частиX_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]
Для создания модели линейной регрессии используем класс LinearRegression из библиотеки Scikit-Learn. Так как раньше задали наличие константы в регрессионном уравнении, для исключения дублирования задаем fit_intercept=False
# Подгонка трендовой модели
model = LinearRegression(fit_intercept=False)
model.fit(X_train, y_train)
# Делаем прогноз
y_fit = pd.DataFrame(
model.predict(X_train),
index=y_train.index,
columns=y_train.columns,
)
y_pred = pd.DataFrame(
model.predict(X_test),
index=y_test.index,
columns=y_test.columns,
)
Выведем график
axs = y_train.plot(color='0.25', figsize=(11, 5))
axs = y_test.plot(color='0.25',ax=axs )
axs = y_fit.plot(color='C0',ax=axs)
axs = y_pred.plot(color='C3',ax=axs)
plt.legend(['train plot', 'test plot','fit plot', 'pred plot'])
Переходим к моделированию сезонности, как остатков модели тренда
Готовим независимые переменные, для этого добавляем номер месяца и убираем целевую переменную - выручку.
X = df
X["Month"] = X.index.month
y = X.pop('Sh01')
X.head()
| month | Month |
|---|---|
| 2015-01 | 1 |
| 2015-02 | 2 |
| 2015-03 | 3 |
| 2015-04 | 4 |
| 2015-05 | 5 |
Делим данные на обучение и тест
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]
Готовим вектор остатков
y_fit = y_fit.stack().squeeze()
y_pred = y_pred.stack().squeeze()
y_resid = y_train - y_fit
Для подгонки остатков используем библиотеку XGBoost, основанную на алгоритме дерева решений и реализующая методы градиентного бустинга класс XGBRegressor
xgb = XGBRegressor()
xgb.fit(X_train, y_resid)
Складываем предсказания тренда и остатков, отражающих сезонность и получаем итоговый прогноз
y_fit_boosted = xgb.predict(X_train) + y_fit
y_pred_boosted = xgb.predict(X_test) + y_pred
Визуализируем полученное
axs = y_train.plot(color='0.25', figsize=(11, 5))
axs = y_test.plot(color='0.25',ax=axs )
axs = y_fit_boosted.unstack().plot(color='C0',ax=axs)
axs = y_pred_boosted.unstack().plot(color='C3',ax=axs)
plt.legend(['train plot', 'test plot','fit plot', 'pred plot'])
Определим функцию для расчета MAPE и выведем метрики качества для обучающей тестовой частей
def MAPE(Y_actual,Y_Predicted):
mape = np.mean(np.abs((Y_actual - Y_Predicted)/Y_actual))*100
return mape
# Рассчитываем метрики качества для обучающей выборки
train_fcst=y_fit_boosted.unstack().loc[:,'Sh01']
mae = metrics.mean_absolute_error(y_train, train_fcst)
mse = metrics.mean_squared_error(y_train, train_fcst)
mape = MAPE(y_train, train_fcst)
rmse = math.sqrt(mse)
r2 = metrics.r2_score(y_train, train_fcst)# Выводим метрики качества для обучающей выборки
print(f'R^2 Score = {r2:.3f}')
print(f'Mean Absolute Percentage Error = {mape:.2f}')
print(f'Mean Absolute Error = {mae:.2f}')
print(f'Mean Squared Error = {mse:.2f}')
print(f'Root Mean Squared Error = {rmse:.2f}')R^2 Score = 0.977 Mean Absolute Percentage Error = 3.68 Mean Absolute Error = 7811.87 Mean Squared Error = 91203651.63 Root Mean Squared Error = 9550.06
# Рассчитываем метрики качества для тестовой выборки
test_fcst=y_pred_boosted.unstack().loc[:,'Sh01']
mae = metrics.mean_absolute_error(y_test, test_fcst)
mse = metrics.mean_squared_error(y_test, test_fcst)
mape = MAPE(y_test, test_fcst)
rmse = math.sqrt(mse)
r2 = metrics.r2_score(y_test, test_fcst)# Выводим метрики для тестовой выборки
print(f'R^2 Score = {r2:.3f}')
print(f'Mean Absolute Percentage Error = {mape:.2f}')
print(f'Mean Absolute Error = {mae:.2f}')
print(f'Mean Squared Error = {mse:.2f}')
print(f'Root Mean Squared Error = {rmse:.2f}')R^2 Score = 0.938 Mean Absolute Percentage Error = 4.86 Mean Absolute Error = 8448.78 Mean Squared Error = 97119241.52 Root Mean Squared Error = 9854.91Переходим к мультипликативной модели сезонности и готовим вектор сезонных
коэффициентовy_sс=y_train/y_fitДля подгонки сезонных коэффициентов также используем библиотеку XGBoost, основанную на алгоритме дерева решений и реализующая методы градиентного бустинга класс XGBRegressorxgb = XGBRegressor() xgb.fit(X_train, y_sk)y_sс_pred_train=xgb.predict(X_train)y_sс_pred_test=xgb.predict(X_test)Умножаем предсказания тренда и сезонных коэффициентов и получаем итоговый прогнозy_fit_mult = y_sc_pred_train*y_fit y_pred_mult = y_sc_pred_test*y_predВыводим графикиaxs = y_train.plot(color='0.25', figsize=(11, 5)) axs = y_test.plot(color='0.25',ax=axs ) axs = pd.DataFrame(y_fit_mult).unstack().plot(color='C0',ax=axs) axs = pd.DataFrame(y_pred_mult).unstack().plot(color='C3',ax=axs) plt.title('multiplicative model') # plt.legend(['train plot', 'test plot','fit plot', 'pred plot'])Выведем метрики качества для обучающей тестовой частей для этого варианта# Рассчитываем метрики качества для обучающей выборкиtrain_fcst=y_fit_mult.unstack().loc[:,'Sh01'] mae = metrics.mean_absolute_error(y_train, train_fcst) mse = metrics.mean_squared_error(y_train, train_fcst) mape = MAPE(y_train, train_fcst) rmse = math.sqrt(mse) r2 = metrics.r2_score(y_train, train_fcst)# Выводим метрики качества для обучающей выборки print(f'R^2 Score = {r2:.3f}') print(f'Mean Absolute Percentage Error = {mape:.2f}') print(f'Mean Absolute Error = {mae:.2f}') print(f'Mean Squared Error = {mse:.2f}') print(f'Root Mean Squared Error = {rmse:.2f}')R^2 Score = 0.987 Mean Absolute Percentage Error = 2.79 Mean Absolute Error = 6150.18 Mean Squared Error = 49410004.75 Root Mean Squared Error = 7029.23# Рассчитываем метрики качества для тестовой выборкиtest_fcst=y_pred_mult.unstack().loc[:,'Sh01'] mae = metrics.mean_absolute_error(y_test, test_fcst) mse = metrics.mean_squared_error(y_test, test_fcst) mape = MAPE(y_test, test_fcst) rmse = math.sqrt(mse) r2 = metrics.r2_score(y_test, test_fcst) # Выводим метрики для тестовой выборкиprint(f'R^2 Score = {r2:.3f}') print(f'Mean Absolute Percentage Error = {mape:.2f}') print(f'Mean Absolute Error = {mae:.2f}') print(f'Mean Squared Error = {mse:.2f}') print(f'Root Mean Squared Error = {rmse:.2f}')R^2 Score = 0.909 Mean Absolute Percentage Error = 4.32 Mean Absolute Error = 8813.82 Mean Squared Error = 142721216.77 Root Mean Squared Error = 11946.60Как видно, обе модели показали примерно равное качество.И в заключении покажем, как можно сделать для первого варианта расчет сразудля нескольких магазин
Загружаем данные
df = pd.read_excel(
"shops8_rev.xlsx",
parse_dates=['month'],index_col='month',).to_period()
df = pd.concat({'Sales': df}, names=[None, 'Shops'], axis=1)
pd.concat([df.head(), df.tail()])
| Sales | ||||||||
|---|---|---|---|---|---|---|---|---|
| Shops | Sh01 | Sh02 | Sh03 | Sh04 | Sh05 | Sh06 | Sh07 | Sh08 |
| month | ||||||||
| 2015-01 | 157408 | 82190 | 174300 | 77120 | 87860 | 51830 | 81790 | 80500 |
| 2015-02 | 127692 | 83010 | 124520 | 83330 | 79210 | 42700 | 62320 | 76530 |
| 2015-03 | 140148 | 91590 | 182450 | 105470 | 98780 | 54490 | 75460 | 110610 |
| 2015-04 | 167924 | 82680 | 183340 | 111520 | 121480 | 53850 | 88780 | 96930 |
| 2015-05 | 184196 | 98930 | 204770 | 123600 | 128200 | 63730 | 92800 | 116680 |
y = df.copy()
dp = DeterministicProcess(
index=df.index, # даты из обучающих данных
constant=True, # фиктивная функция для смещения (y_intercept)
order=2, # квадратичный тренд
drop=True, # если необходимо, отбросить члены, чтобы избежать коллинеарности
)
in_sample` создает функции для дат, указанных в аргументе `index`
X = dp.in_sample()
X.head()
| month | const | trend | trend_squared |
|---|---|---|---|
| 2015-01 | 1.0 | 1.0 | 1.0 |
| 2015-02 | 1.0 | 2.0 | 4.0 |
| 2015-03 | 1.0 | 3.0 | 9.0 |
| 2015-04 | 1.0 | 4.0 | 16.0 |
| 2015-05 | 1.0 | 5.0 | 25.0 |
Делим данные на обучение и тест
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]
# Подгонка трендовой модели
model = LinearRegression(fit_intercept=False)
model.fit(X_train, y_train)
# Делаем прогноз
y_fit = pd.DataFrame(
model.predict(X_train),
index=y_train.index,
columns=y_train.columns,
)
y_pred = pd.DataFrame(
model.predict(X_test),
index=y_test.index,
columns=y_test.columns,
)
Выводим график
axs = y_train.plot(color='0.25', figsize=(11, 15),subplots=True, sharex=True,
title=['Sh01','Sh02','Sh03','Sh04','Sh05','Sh06','Sh07','Sh08',], )
axs = y_test.plot(color='0.25', subplots=True, sharex=True, ax=axs)
axs = y_fit.plot(color='C0', subplots=True, sharex=True, ax=axs)
axs = y_pred.plot(color='C3', subplots=True, sharex=True, ax=axs)
for ax in axs: ax.legend([])
_ = plt.suptitle("Trends")
Проводим преобразование данных, деление на обучающую и тестовую выборки, рассчитываем остатки модели, обучаем по ним модель градиентного бустинга и в конце сложив предсказания тренда и остатков, отражающих сезонность и получаем итоговые прогнозы с их последующей визуализацией.
X = df.stack()
display(X.head())
y = X.pop('Sales')
X = X.reset_index('Shops')
X['Shops'] , _= X['Shops'].factorize()
X["Month"] = X.index.month
X.head()
| month |
|
Month | |
|---|---|---|---|
| 2015-01 | 0 | 1 | |
| 2015-01 | 1 | 1 | |
| 2015-01 | 2 | 1 | |
| 2015-01 | 3 | 1 | |
| 2015-01 | 4 | 1 |
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]
y_fit = y_fit.stack().squeeze()
y_pred = y_pred.stack().squeeze()
y_resid = y_train - y_fit
xgb = XGBRegressor()
xgb.fit(X_train, y_resid)
y_fit_boosted = xgb.predict(X_train) + y_fit
y_pred_boosted = xgb.predict(X_test) + y_pred
axs = y_train.unstack(['Shops']).plot(
color='0.25', figsize=(11, 20), subplots=True, sharex=True,
title=['Sh01','Sh02','Sh03','Sh04','Sh05','Sh06','Sh07','Sh08',],
)
axs = y_test.unstack(['Shops']).plot(
color='0.25', subplots=True, sharex=True, ax=axs,
)
axs = y_fit_boosted.unstack(['Shops']).plot(
color='C0', subplots=True, sharex=True, ax=axs,
)
axs = y_pred_boosted.unstack(['Shops']).plot(
color='C3', subplots=True, sharex=True, ax=axs,
)
for ax in axs: ax.legend([])
Комментариев нет:
Отправить комментарий