Автор материала: Дмитрий Сергеев, Zeptolab. Материал распространяется на условиях лицензии Creative Commons CC BY-NC-SA 4.0. Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала.
Заполните код в клетках и выберите ответы в веб-форме.
Импортируем необходимые библиотеки
In [83]:
import warnings
warnings.filterwarnings('ignore')
import random
random.seed(42)
import seaborn as sns
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import TimeSeriesSplit, cross_val_score
import matplotlib.pyplot as plt
%matplotlib inline
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
Для работы возьмем реальные часовые данные по суммарному просмотру рекламы в одном из приложений.
In [84]:
df = pd.read_csv('../../mlcourse_open//data/ads_hour.csv',index_col=['Date'],
parse_dates=['Date'])
In [85]:
with plt.style.context('bmh'):
plt.figure(figsize=(15, 8))
plt.title('Ads watched (hour ticks)')
plt.plot(df.ads);
Тут есть и тренд, и различная сезонность, и немножко выбросов.
В этом домашнем задании мы не будем пытаться подобрать наилучший ARIMA-процесс, порождающий наш ряд, и не затронем библиотеку Prophet от Facebook. Работа с этими методами очень подробна расписана во второй и третьей части ноутбуков по временным рядам в репозитории открытого курса по машинному обучению.
А здесь мы сосредоточимся на построении и отборе признаков для временных рядов.
Возьмем две полезные функции из лекции – для расчета среднего значения по таргету и для подготовки обучающего и тестового датасетов:
In [86]:
def code_mean(data, cat_feature, real_feature):
"""
Возвращает словарь, где ключами являются уникальные категории признака cat_feature,
а значениями - средние по real_feature
"""
return dict(data.groupby(cat_feature)[real_feature].mean())
def prepareData(data, lag_start=5, lag_end=14, test_size=0.15, target_encoding=False):
"""
series: pd.DataFrame
dataframe with timeseries
lag_start: int
initial step back in time to slice target variable
example - lag_start = 1 means that the model
will see yesterday's values to predict today
lag_end: int
final step back in time to slice target variable
example - lag_end = 4 means that the model
will see up to 4 days back in time to predict today
test_size: float
size of the test dataset after train/test split as percentage of dataset
target_encoding: boolean
if True - add target averages to the dataset
"""
data = pd.DataFrame(data.copy())
data.columns = ["y"]
# считаем индекс в датафрейме, после которого начинается тестовыый отрезок
test_index = int(len(data) * (1 - test_size))
# добавляем лаги исходного ряда в качестве признаков
for i in range(lag_start, lag_end):
data["lag_{}".format(i)] = data.y.shift(i)
data.index = data.index.to_datetime()
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5,6])*1
if target_encoding:
# считаем средние только по тренировочной части, чтобы избежать лика
data['weekday_average'] = list(map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday))
data["hour_average"] = list(map(code_mean(data[:test_index], 'hour', "y").get, data.hour))
# выкидываем закодированные средними признаки
data.drop(["hour", "weekday"], axis=1, inplace=True)
data = data.dropna()
data = data.reset_index(drop=True)
# разбиваем весь датасет на тренировочную и тестовую выборку
X_train = data.loc[:test_index].drop(["y"], axis=1)
y_train = data.loc[:test_index]["y"]
X_test = data.loc[test_index:].drop(["y"], axis=1)
y_test = data.loc[test_index:]["y"]
return X_train, X_test, y_train, y_test
Воспользуйтесь предложенными функциями и приготовьте из исходного датасета df необходимые для обучения датафреймы. Для тестирования отложите 30% данных, в качестве начального лага возьмите значения временного ряда 12 часов назад, а в качестве последнего лага - значения ряда 48 часов назад. Таким образом модель будет способна строить предсказания на 12 часов вперёд, имея фактические наблюдения за предыдущие полтора дня.
Также отмасштабируйте признаки, для этого воспользуйстесь StandardScaler и создайте переменные X_train_scaled и X_test_scaled. Не забудьте, что обучать scaler нужно только на тренировочной выборке, чтобы информация о средних отклонениях не просочилась из будущего в прошлое.
In [87]:
scaler = StandardScaler()
X_train, X_test, y_train, y_test = prepareData(df, test_size=0.3, lag_start=12, lag_end=48)
In [88]:
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)
Теперь обучите простую линейную регрессию на полученных данных:
In [89]:
lin_reg = LinearRegression()
In [90]:
lin_reg.fit(X_train, y_train)
Out[90]:
Оценим качество работы линейной регрессии на тренировочной части при помощи кросс-валидации. Для этого необходимо сначала создать объект-генератор разбиений для временного ряда при помощи TimeSeriesSplit. Число фолдов задайте равным 5. Затем воспользуйтесь методом cross_val_score, передав туда в качестве параметра cv созданный генератор разбиений. Метрикой качества будет выступать neg_mean_absolute_error.
Не забудьте взять среднее от полученного значения и домножить его на -1.
In [91]:
tscv = TimeSeriesSplit(n_splits=5)
scores = cross_val_score(estimator=lin_reg, cv=tscv, X=X_train, y=y_train, scoring='neg_mean_absolute_error')
-np.average(scores)
Out[91]:
Вопрос 1. Чему равно среднее значение MAE на кросс-валидации?
Теперь посмотрим на результаты работы линейной регрессии на отложенной выборке. В этом нам поможет функция plotModelResults. В ней Вам нужно дописать часть, ответственную за расчет доверительных интервалов. Для расчета отклонения воспользуйтесь только что реализованным способом оценки качества на кросс-валидации – оттуда понадобится среднее значение MAE и стандартное отклонение этой величины. Доверительные интервалы строятся по формуле $\hat{y} \pm (\text{error} + scale * \text{standard deviation})$
In [92]:
def plotModelResults(model, X_train, X_test, plot_intervals=False):
"""
Строит график прогнозных и фактических значений, а также доверительных интервалов прогноза
"""
# получаем предсказания по модели
prediction = model.predict(X_test)
plt.figure(figsize=(20, 7))
plt.plot(prediction, "g", label="prediction", linewidth=2.0)
plt.plot(y_test.values, label="actual", linewidth=2.0)
if plot_intervals:
cv = TimeSeriesSplit(n_splits=5)
scores = cross_val_score(estimator=model,
cv=tscv,
X=X_train,
y=y_train,
scoring='neg_mean_absolute_error')
mae = -np.mean(scores)
deviation = np.std(scores)
scale = 1.96
lower = prediction - (mae + scale*deviation)
upper = prediction + (mae + scale*deviation)
plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
plt.plot(upper, "r--", alpha=0.5)
pass
mae = mean_absolute_error(prediction, y_test)
mape = mean_absolute_percentage_error(prediction, y_test)
plt.title("MAE {}, MAPE {}%".format(round(mae), round(mape, 2)))
plt.legend(loc="best")
plt.grid(True);
Для визуализации коэффициентов модели воспользуйтесь следующими функциями:
In [93]:
def getCoefficients(model):
"""
Возвращает отсортированные по абсолютному значению коэффициенты модели
"""
coefs = pd.DataFrame(model.coef_, X_train.columns)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
return coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)
def plotCoefficients(model):
"""
Отрисовывает отсортированные по абсолютному значению коэффициенты модели
"""
coefs = getCoefficients(model)
plt.figure(figsize=(20, 7))
coefs.coef.plot(kind='bar')
plt.grid(True, axis='y')
plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed');
In [94]:
plotModelResults(lin_reg, X_train, X_test, plot_intervals=True)
plotCoefficients(lin_reg)
Замечательно, интервалы построены, прогнозы достаточно точные, и всё бы хорошо, если бы не одно но. В модели сейчас находится множество признаков и, возможно, какие-то из них лишние, а если их убрать, качество модели серьезно не пострадает. Чтобы убедиться, что у нас есть лишние признаки, постройте тепловую карту корреляций по X_train, поспользовавшись функцией heatmap из библиотеки seaborn:
In [95]:
plt.figure(figsize=(12, 10))
sns.heatmap(X_train);
Действительно, налицо множество корреляций между нашими признаками, и даже видна "сезонность" в этих корреляциях на каждом 24-м лаге. Попробуем регуляризовать нашу модель и убрать несколько лишних переменных.
Обучите на полученных отмасштабированных данных Лассо-регрессию на кросс-валидации (LassoCV), снова передав в качестве параметра cv созданный ранее генератор разбиений.
Постройте график прогнозов с интервалами и убедитесь, что ошибка на тестовой выборке практически не изменилась.
In [96]:
lasso = LassoCV(cv=tscv)
lasso.fit(X_train_scaled, y_train);
In [97]:
plotModelResults(lasso, X_train_scaled, X_test_scaled, plot_intervals=True)
Замечательно, качество нас всё ещё устраивает, интервалы по-прежднему адекватные, посмотрим, сильно ли мы упростили модель. Воспользуйтесь предложенной функцией getCoefficients и скажите, сколько признаков вошли в модель с зануленными весами.
In [98]:
getCoefficients(lasso)
Out[98]:
Вопрос 2. Сколько коэффициентов линейной модели равны нулю (с точностью до 6-го знака)?
Итак, признаковое пространство хоть немного, но удалось сократить. Что если пойти дальше и сжать наши линейно зависимые признаки в более компактное представление? Для этого воспользуемся методом главных компонент.
In [99]:
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
def plotPCA():
"""
График накопленного процента объясненной дисперсии по компонентам
"""
features = range(pca.n_components_)
variance = np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)
plt.figure(figsize=(20, 10))
plt.bar(features, variance)
# дополнительно отметим уровень, при котором объяснены 95% дисперсии
plt.hlines(y = 95, xmin=0, xmax=len(features), linestyles='dashed', colors='red')
plt.xlabel('PCA components')
plt.ylabel('variance')
plt.xticks(features)
plt.show()
In [107]:
pca = PCA(0.95)
pca.fit(X_train_scaled)
variances = np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)
plt.plot(variances)
variances
Out[107]:
In [64]:
pca.n_components_
Out[64]:
Вопрос 3. Какое минимальное число компонент необходимо взять, чтобы объяснить минимум 95% всей дисперсии тренировочного датасета?
Снова создайте объект pca, указав при этом оптимальное число компонент (объясняющее минимум 95% изменчивости). Затем создайте две новых переменных – pca_features_train и pca_features_test, присвоив им, соответственно, преобразованные при помощи PCA отмасштабированные датафреймы.
In [80]:
pca = PCA(n_components=6)
pca.fit(X_train_scaled)
pca_features_train = pca.transform(X_train_scaled)
pca_features_test = pca.transform(X_test_scaled)
In [81]:
pca_features_train.shape, pca_features_test.shape
Out[81]:
На полученных данных снова обучите простую линейную регрессию и постройте график прогноза.
In [82]:
lin_reg = LinearRegression()
lin_reg.fit(pca_features_train, y_train)
plotModelResults(lin_reg, pca_features_train, pca_features_test, plot_intervals=True)
Вопрос 4. Какая средняя абсолютная ошибка получилась у линейной регрессии, обученной на нескольких главных компонентах?