В этом задании мы на примерах увидим, как переобучаются линейные модели, разберем, почему так происходит, и выясним, как диагностировать и контролировать переобучение.
Во всех ячейках, где написан комментарий с инструкциями, нужно написать код, выполняющий эти инструкции. Остальные ячейки с кодом (без комментариев) нужно просто выполнить. Кроме того, в задании требуется отвечать на вопросы; ответы нужно вписывать после выделенного слова "Ответ:".
Напоминаем, что посмотреть справку любого метода или функции (узнать, какие у нее аргументы и что она делает) можно с помощью комбинации Shift+Tab. Нажатие Tab после имени объекта и точки позволяет посмотреть, какие методы и переменные есть у этого объекта.
In [3]:
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
Мы будем работать с датасетом "bikes_rent.csv", в котором по дням записаны календарная информация и погодные условия, характеризующие автоматизированные пункты проката велосипедов, а также число прокатов в этот день. Последнее мы будем предсказывать; таким образом, мы будем решать задачу регрессии.
Загрузите датасет с помощью функции pandas.read_csv в переменную df. Выведите первые 5 строчек, чтобы убедиться в корректном считывании данных:
In [4]:
# (0 баллов)
# Считайте данные и выведите первые 5 строк
df = pd.read_csv('bikes_rent.csv')
df.head(5)
Out[4]:
Для каждого дня проката известны следующие признаки (как они были указаны в источнике данных):
Итак, у нас есть вещественные, бинарные и номинальные (порядковые) признаки, и со всеми из них можно работать как с вещественными. С номинальныеми признаками тоже можно работать как с вещественными, потому что на них задан порядок. Давайте посмотрим на графиках, как целевой признак зависит от остальных
In [3]:
fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(15, 10))
for idx, feature in enumerate(df.columns[:-1]):
df.plot(feature, "cnt", subplots=True, kind="scatter", ax=axes[idx / 4, idx % 4])
Блок 1. Ответьте на вопросы (каждый 0.5 балла):
Давайте более строго оценим уровень линейной зависимости между признаками и целевой переменной. Хорошей мерой линейной зависимости между двумя векторами является корреляция Пирсона. В pandas ее можно посчитать с помощью двух методов датафрейма: corr и corrwith. Метод df.corr вычисляет матрицу корреляций всех признаков из датафрейма. Методу df.corrwith нужно подать еще один датафрейм в качестве аргумента, и тогда он посчитает попарные корреляции между признаками из df и этого датафрейма.
In [5]:
# Код 1.1 (0.5 балла)
# Посчитайте корреляции всех признаков, кроме последнего, с последним с помощью метода corrwith:
df[df.columns[:-1]].corrwith(df["cnt"])
Out[5]:
В выборке есть признаки, коррелирующие с целевым, а значит, задачу можно решать линейными методами.
По графикам видно, что некоторые признаки похожи друг на друга. Поэтому давайте также посчитаем корреляции между вещественными признаками.
In [5]:
# Код 1.2 (0.5 балла)
# Посчитайте попарные корреляции между признаками temp, atemp, hum, windspeed(mph), windspeed(ms) и cnt
# с помощью метода corr:
df[['temp','atemp', 'hum','windspeed(mph)','windspeed(ms)','cnt']].corr()
Out[5]:
На диагоналях, как и полагается, стоят единицы. Однако в матрице имеются еще две пары сильно коррелирующих столбцов: temp и atemp (коррелируют по своей природе) и два windspeed (потому что это просто перевод одних единиц в другие). Далее мы увидим, что этот факт негативно сказывается на обучении линейной модели.
Напоследок посмотрим средние признаков (метод mean), чтобы оценить масштаб признаков и доли 1 у бинарных признаков.
In [6]:
# Код 1.3 (0.5 балла)
# Выведите средние признаков
print df.mean()
Признаки имеют разный масштаб, значит для дальнейшей работы нам лучше нормировать матрицу объекты-признаки.
Итак, в наших данных один признак дублирует другой, и есть еще два очень похожих. Конечно, мы могли бы сразу удалить дубликаты, но давайте посмотрим, как бы происходило обучение модели, если бы мы не заметили эту проблему.
Для начала проведем масштабирование, или стандартизацию признаков: из каждого признака вычтем его среднее и поделим на стандартное отклонение. Это можно сделать с помощью метода scale.
Кроме того, нужно перемешать выборку, это потребуется для кросс-валидации.
In [9]:
from sklearn.preprocessing import scale
from sklearn.utils import shuffle
In [10]:
df_shuffled = shuffle(df, random_state=123)
X = scale(df_shuffled[df_shuffled.columns[:-1]])
y = df_shuffled["cnt"]
Давайте обучим линейную регрессию на наших данных и посмотрим на веса признаков.
In [6]:
from sklearn.linear_model import LinearRegression
In [11]:
# Код 2.1 (1 балл)
# Создайте объект линейного регрессора, обучите его на всех данных и выведите веса модели
# (веса хранятся в переменной coef_ класса регрессора).
# Можно выводить пары (название признака, вес), воспользовавшись функцией zip, встроенной в язык python
# Названия признаков хранятся в переменной df.columns
alg_lnr = LinearRegression()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)
Мы видим, что веса при линейно-зависимых признаках по модулю значительно больше, чем при других признаках.
Чтобы понять, почему так произошло, вспомним аналитическую формулу, по которой вычисляются веса линейной модели в методе наименьших квадратов:
$w = (X^TX)^{-1} X^T y$.
Если в X есть коллинеарные (линейно-зависимые) столбцы, матрица $X^TX$ становится вырожденной, и формула перестает быть корректной. Чем более зависимы признаки, тем меньше определитель этой матрицы и тем хуже аппроксимация $Xw \approx y$. Такая ситуацию называют проблемой мультиколлинеарности, вы обсуждали ее на лекции.
С парой temp-atemp чуть менее коррелирующих переменных такого не произошло, однако на практике всегда стоит внимательно следить за коэффициентами при похожих признаках.
Решение проблемы мультиколлинеарности состоит в регуляризации линейной модели. К оптимизируемому функционалу прибавляют L1 или L2 норму весов, умноженную на коэффициент регуляризации $\alpha$. В первом случае метод называется Lasso, а во втором --- Ridge. Подробнее об этом также рассказано в лекции.
Обучите регрессоры Ridge и Lasso с параметрами по умолчанию и убедитесь, что проблема с весами решилась.
In [12]:
from sklearn.linear_model import Lasso, Ridge
In [12]:
# Код 2.2 (0.5 балла)
# Обучите линейную модель с L1-регуляризацией
alg_lnr = Lasso()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)
In [13]:
# Код 2.3 (0.5 балла)
# Обучите линейную модель с L2-регуляризацией
alg_lnr = Ridge()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)
Блок 2. Поясните, каким образом введение регуляризации решает проблему с весами и мультиколлинеарностью.
Ваш ответ (1 балл): регуляризация изменяет функцию ошибку, таким образом модель "наказывается" за слишком большие коэф.. Коэф. перед скоростью там получился просто коллосальный на много порядков больше, чем перед остальными, регуляризация же решила эту проблему
В отличие от L2-регуляризации, L1 обнуляет веса при некоторых признаках. Объяснение данному факту дается в одной из лекций курса.
Давайте пронаблюдаем, как меняются веса при увеличении коэффициента регуляризации $\alpha$ (в лекции коэффициент при регуляризаторе мог быть обозначен другой буквой).
In [14]:
# Код 3.1 (1 балл)
alphas = np.arange(1, 500, 50)
coefs_lasso = np.zeros((alphas.shape[0], X.shape[1])) # матрица весов размера (число регрессоров) x (число признаков)
coefs_ridge = np.zeros((alphas.shape[0], X.shape[1]))
# Для каждого значения коэффициента из alphas обучите регрессор Lasso
# и запишите веса в соответствующую строку матрицы coefs_lasso (вспомните встроенную в python функцию enumerate),
# а затем обучите Ridge и запишите веса в coefs_ridge.
print coefs_lasso.shape
num = 0
for alpha in alphas:
modelL1 = Lasso(alpha = alpha).fit(X, y)
modelL2 = Ridge(alpha = alpha).fit(X, y)
coefs_lasso[num] = modelL1.coef_
coefs_ridge[num] = modelL2.coef_
num = num + 1
Визуализируем динамику весов при увеличении параметра регуляризации:
In [15]:
plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso")
plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_ridge.T, df.columns):
plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Ridge")
plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
if feature == 'windspeed(mph)' or feature == 'windspeed(ms)' or feature == 'yr':
plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso - temp")
Out[15]:
Ответы на следующие вопросы можно давать, глядя на графики или выводя коэффициенты на печать.
Блок 3. Ответьте на вопросы (каждый 0.25 балла):
Далее будем работать с Lasso.
Итак, мы видим, что при изменении alpha модель по-разному подбирает коэффициенты признаков. Нам нужно выбрать наилучшее alpha.
Для этого, во-первых, нам нужна метрика качества. Будем использовать в качестве метрики сам оптимизируемый функционал метода наименьших квадратов, то есть Mean Square Error.
Во-вторых, нужно понять, на каких данных эту метрику считать. Нельзя выбирать alpha по значению MSE на обучающей выборке, потому что тогда мы не сможем оценить, как модель будет делать предсказания на новых для нее данных. Если мы выберем одно разбиение выборки на обучающую и тестовую (это называется holdout), то настроимся на конкретные "новые" данные, и вновь можем переобучиться. Поэтому будем делать несколько разбиений выборки, на каждом пробовать разные значения alpha, а затем усреднять MSE. Удобнее всего делать такие разбиения кросс-валидацией, то есть разделить выборку на K частей, или блоков, и каждый раз брать одну из них как тестовую, а из оставшихся блоков составлять обучающую выборку.
Делать кросс-валидацию для регрессии в sklearn совсем просто: для этого есть специальный регрессор, LassoCV, который берет на вход список из alpha и для каждого из них вычисляет MSE на кросс-валидации. После обучения (если оставить параметр cv=3 по умолчанию) регрессор будет содержать переменную mse_path_, матрицу размера len(alpha) x k, k = 3 (число блоков в кросс-валидации), содержащую значения MSE на тесте для соответствующих запусков. Кроме того, в переменной alpha_ будет храниться выбранное значение параметра регуляризации, а в coef_, традиционно, обученные веса, соответствующие этому alpha_.
Обратите внимание, что регрессор может менять порядок, в котором он проходит по alphas; для сопоставления с матрицей MSE лучше использовать переменную регрессора alphas_.
In [14]:
from sklearn.linear_model import LassoCV
In [18]:
# Код 3.2 (1 балл)
# Обучите регрессор LassoCV на всех параметрах регуляризации из alpha
# Постройте график _усредненного_ по строкам MSE в зависимости от alpha.
# Выведите выбранное alpha, а также пары "признак-коэффициент" для обученного вектора коэффициентов
alphas = np.arange(1, 100, 5)
clf = LassoCV(alphas = alphas)
model = clf.fit(X, y)
print model.mse_path_.shape
plt.plot(model.alphas_, model.mse_path_.mean(axis =1))
plt.xlabel('alpha')
plt.ylabel('mean MSE')
plt.show()
print 'выбранное alpha = ', model.alpha_,'\n', zip(df.columns[0:(df.shape[1]-1)], clf.coef_)
Итак, мы выбрали некоторый параметр регуляризации. Давайте посмотрим, какие бы мы выбирали alpha, если бы делили выборку только один раз на обучающую и тестовую, то есть рассмотрим траектории MSE, соответствующие отдельным блокам выборки.
In [19]:
# Код 3.3 (1 балл)
# Выведите значения alpha, соответствующие минимумам MSE на каждом разбиении (то есть по столбцам).
# На трех отдельных графиках визуализируйте столбцы .mse_path_
print model.alphas_[np.argmin(model.mse_path_[:,0])]
plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,0], 'r')
plt.xlabel('alpha')
plt.ylabel('1st MSE')
plt.legend()
plt.axhline(y= np.amin(model.mse_path_[:,0],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,0])]))
print model.alphas_[np.argmin(model.mse_path_[:,1])]
plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,1], color = 'y')
plt.xlabel('alpha')
plt.ylabel('2nd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,1],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,1])]))
print model.alphas_[np.argmin(model.mse_path_[:,2])]
plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,2], color = 'g')
plt.xlabel('alpha')
plt.ylabel('3rd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,2],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,2])]))
Out[19]:
На каждом разбиении оптимальное значение alpha свое, и ему соответствует большое MSE на других разбиениях. Получается, что мы настраиваемся на конкретные обучающие и контрольные выборки. При выборе alpha на кросс-валидации мы выбираем нечто "среднее", что будет давать приемлемое значение метрики на разных разбиениях выборки.
Наконец, как принято в анализе данных, давайте проинтерпретируем результат.
Блок 4. Ответьте на вопросы (каждый 0.5 балла):
In [19]:
zip(df.columns[0:(df.shape[1]-1)], clf.coef_)
Out[19]:
Итак, мы посмотрели, как можно следить за адекватностью линейной модели, как отбирать признаки и как грамотно, по возможности не настраиваясь на какую-то конкретную порцию данных, подбирать коэффициент регуляризации.
Стоит отметить, что с помощью кросс-валидации удобно подбирать лишь небольшое число параметров (1, 2, максимум 3), потому что для каждой допустимой их комбинации нам приходится несколько раз обучать модель, а это времязатратный процесс, особенно если нужно обучаться на больших объемах данных.