Практическое задание к уроку 1 (2 неделя).

Линейная регрессия: переобучение и регуляризация

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

Во всех ячейках, где написан комментарий с инструкциями, нужно написать код, выполняющий эти инструкции. Остальные ячейки с кодом (без комментариев) нужно просто выполнить. Кроме того, в задании требуется отвечать на вопросы; ответы нужно вписывать после выделенного слова "Ответ:".

Напоминаем, что посмотреть справку любого метода или функции (узнать, какие у нее аргументы и что она делает) можно с помощью комбинации Shift+Tab. Нажатие Tab после имени объекта и точки позволяет посмотреть, какие методы и переменные есть у этого объекта.


In [1]:
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 [2]:
df = pd.read_csv('bikes_rent.csv')
df.head()


Out[2]:
season yr mnth holiday weekday workingday weathersit temp atemp hum windspeed(mph) windspeed(ms) cnt
0 1 0 1 0 6 0 2 14.110847 18.18125 80.5833 10.749882 4.805490 985
1 1 0 1 0 0 0 2 14.902598 17.68695 69.6087 16.652113 7.443949 801
2 1 0 1 0 1 1 1 8.050924 9.47025 43.7273 16.636703 7.437060 1349
3 1 0 1 0 2 1 1 8.200000 10.60610 59.0435 10.739832 4.800998 1562
4 1 0 1 0 3 1 1 9.305237 11.46350 43.6957 12.522300 5.597810 1600

Для каждого дня проката известны следующие признаки (как они были указаны в источнике данных):

  • season: 1 - весна, 2 - лето, 3 - осень, 4 - зима
  • yr: 0 - 2011, 1 - 2012
  • mnth: от 1 до 12
  • holiday: 0 - нет праздника, 1 - есть праздник
  • weekday: от 0 до 6
  • workingday: 0 - нерабочий день, 1 - рабочий день
  • weathersit: оценка благоприятности погоды от 1 (чистый, ясный день) до 4 (ливень, туман)
  • temp: температура в Цельсиях
  • atemp: температура по ощущениям в Цельсиях
  • hum: влажность
  • windspeed(mph): скорость ветра в милях в час
  • windspeed(ms): скорость ветра в метрах в секунду
  • cnt: количество арендованных велосипедов (это целевой признак, его мы будем предсказывать)

Итак, у нас есть вещественные, бинарные и номинальные (порядковые) признаки, и со всеми из них можно работать как с вещественными. С номинальныеми признаками тоже можно работать как с вещественными, потому что на них задан порядок. Давайте посмотрим на графиках, как целевой признак зависит от остальных


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 балла):

  1. Каков характер зависимости числа прокатов от месяца?
    • ответ: не линейная зависимость
  2. Укажите один или два признака, от которых число прокатов скорее всего зависит линейно
    • ответ: yr, temp, atamp,

Давайте более строго оценим уровень линейной зависимости между признаками и целевой переменной. Хорошей мерой линейной зависимости между двумя векторами является корреляция Пирсона. В pandas ее можно посчитать с помощью двух методов датафрейма: corr и corrwith. Метод df.corr вычисляет матрицу корреляций всех признаков из датафрейма. Методу df.corrwith нужно подать еще один датафрейм в качестве аргумента, и тогда он посчитает попарные корреляции между признаками из df и этого датафрейма.


In [4]:
# Код 1.1 (0.5 балла)
# Посчитайте корреляции всех признаков, кроме последнего, с последним с помощью метода corrwith:
df[df.columns[:-1]].corrwith(df[df.columns[-1]])


Out[4]:
season            0.406100
yr                0.566710
mnth              0.279977
holiday          -0.068348
weekday           0.067443
workingday        0.061156
weathersit       -0.297391
temp              0.627494
atemp             0.631066
hum              -0.100659
windspeed(mph)   -0.234545
windspeed(ms)    -0.234545
dtype: float64

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

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


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 hum windspeed(mph) windspeed(ms) cnt
temp 1.000000 0.991702 0.126963 -0.157944 -0.157944 0.627494
atemp 0.991702 1.000000 0.139988 -0.183643 -0.183643 0.631066
hum 0.126963 0.139988 1.000000 -0.248489 -0.248489 -0.100659
windspeed(mph) -0.157944 -0.183643 -0.248489 1.000000 1.000000 -0.234545
windspeed(ms) -0.157944 -0.183643 -0.248489 1.000000 1.000000 -0.234545
cnt 0.627494 0.631066 -0.100659 -0.234545 -0.234545 1.000000

На диагоналях, как и полагается, стоят единицы. Однако в матрице имеются еще две пары сильно коррелирующих столбцов: temp и atemp (коррелируют по своей природе) и два windspeed (потому что это просто перевод одних единиц в другие). Далее мы увидим, что этот факт негативно сказывается на обучении линейной модели.

Напоследок посмотрим средние признаков (метод mean), чтобы оценить масштаб признаков и доли 1 у бинарных признаков.


In [6]:
# Код 1.3 (0.5 балла)
# Выведите средние признаков
df.describe().iloc[[1]]


Out[6]:
season yr mnth holiday weekday workingday weathersit temp atemp hum windspeed(mph) windspeed(ms) cnt
mean 2.49658 0.500684 6.519836 0.028728 2.997264 0.683995 1.395349 20.310776 23.717699 62.789406 12.762576 5.70522 4504.348837

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

Проблема первая: коллинеарные признаки

Итак, в наших данных один признак дублирует другой, и есть еще два очень похожих. Конечно, мы могли бы сразу удалить дубликаты, но давайте посмотрим, как бы происходило обучение модели, если бы мы не заметили эту проблему.

Для начала проведем масштабирование, или стандартизацию признаков: из каждого признака вычтем его среднее и поделим на стандартное отклонение. Это можно сделать с помощью метода scale.

Кроме того, нужно перемешать выборку, это потребуется для кросс-валидации.


In [7]:
from sklearn.preprocessing import scale
from sklearn.utils import shuffle

In [8]:
df_shuffled = shuffle(df, random_state=123)
X = scale(df_shuffled[df_shuffled.columns[:-1]])
y = df_shuffled["cnt"]

Давайте обучим линейную регрессию на наших данных и посмотрим на веса признаков.


In [9]:
from sklearn.linear_model import LinearRegression

In [10]:
# Код 2.1 (1 балл)
# Создайте объект линейного регрессора, обучите его на всех данных и выведите веса модели 
# (веса хранятся в переменной coef_ класса регрессора).
# Можно выводить пары (название признака, вес), воспользовавшись функцией zip, встроенной в язык python
# Названия признаков хранятся в переменной df.columns
reg = LinearRegression()
reg.fit(X, y)
for weight, name in zip(reg.coef_, df.columns[:-1]):
    print('%s - %s' % (name, weight))


season - 570.870775171
yr - 1021.96892469
mnth - -141.305340472
holiday - -86.7575289679
weekday - 137.224644859
workingday - 56.3908513528
weathersit - -330.232193913
temp - 367.474410473
atemp - 585.553104521
hum - -145.608398517
windspeed(mph) - 1.24573948711e+13
windspeed(ms) - -1.24573948713e+13

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

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

$w = (X^TX)^{-1} X^T y$.

Если в X есть коллинеарные (линейно-зависимые) столбцы, матрица $X^TX$ становится вырожденной, и формула перестает быть корректной. Чем более зависимы признаки, тем меньше определитель этой матрицы и тем хуже аппроксимация $Xw \approx y$. Такая ситуацию называют проблемой мультиколлинеарности, вы обсуждали ее на лекции.

С парой temp-atemp чуть менее коррелирующих переменных такого не произошло, однако на практике всегда стоит внимательно следить за коэффициентами при похожих признаках.

Решение проблемы мультиколлинеарности состоит в регуляризации линейной модели. К оптимизируемому функционалу прибавляют L1 или L2 норму весов, умноженную на коэффициент регуляризации $\alpha$. В первом случае метод называется Lasso, а во втором --- Ridge. Подробнее об этом также рассказано в лекции.

Обучите регрессоры Ridge и Lasso с параметрами по умолчанию и убедитесь, что проблема с весами решилась.


In [11]:
from sklearn.linear_model import Lasso, Ridge

In [12]:
# Код 2.2 (0.5 балла)
# Обучите линейную модель с L1-регуляризацией и выведите веса
las = Lasso()
las.fit(X, y)
for weight, name in zip(las.coef_, df.columns[:-1]):
    print('%s - %s' % (name, weight))


season - 560.241616031
yr - 1019.46349407
mnth - -128.730627037
holiday - -86.1527813337
weekday - 137.347893905
workingday - 55.2123706414
weathersit - -332.369856962
temp - 376.36323621
atemp - 576.530793505
hum - -144.129155003
windspeed(mph) - -197.139689402
windspeed(ms) - -2.8050001994e-08

In [13]:
# Код 2.3 (0.5 балла)
# Обучите линейную модель с L2-регуляризацией и выведите веса
rid = Ridge()
rid.fit(X, y)
for weight, name in zip(rid.coef_, df.columns[:-1]):
    print('%s - %s' % (name, weight))


season - 563.064572252
yr - 1018.94837879
mnth - -131.873320282
holiday - -86.7460979971
weekday - 138.005111179
workingday - 55.9031103751
weathersit - -332.349788499
temp - 386.457889192
atemp - 566.34704706
hum - -145.071327299
windspeed(mph) - -99.2594410818
windspeed(ms) - -99.2594411543

Проблема вторая: неинформативные признаки

В отличие от L2-регуляризации, L1 обнуляет веса при некоторых признаках. Объяснение данному факту дается в одной из лекций курса.

Давайте пронаблюдаем, как меняются веса при увеличении коэффициента регуляризации $\alpha$ (в лекции коэффициент при регуляризаторе мог быть обозначен другой буквой).


In [17]:
# Код 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.
for idx, alpha in enumerate(alphas):
    las = Lasso(alpha=alpha)
    las.fit(X, y)
    coefs_lasso[idx] = las.coef_
    
    rid = Ridge(alpha=alpha)
    rid.fit(X, y)
    coefs_ridge[idx] = rid.coef_

Визуализируем динамику весов при увеличении параметра регуляризации:


In [18]:
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")


Out[18]:
<matplotlib.text.Text at 0x115a6b7b8>

Ответы на следующие вопросы можно давать, глядя на графики или выводя коэффициенты на печать.

Блок 2. Ответьте на вопросы (каждый 0.25 балла):

  1. Какой регуляризатор (Ridge или Lasso) агрессивнее уменьшает веса при одном и том же alpha?
    • Ответ: Lasso
  2. Что произойдет с весами Lasso, если alpha сделать очень большим? Поясните, почему так происходит.
    • Ответ: Веса обнуляются. Так происходит потому что функция L1 не гладкая и плохо поддается оптимизации, поетому значения скативаются в точку локального мінімума где значение все весов минимальное.
  3. Можно ли утверждать, что Lasso исключает один из признаков windspeed при любом значении alpha > 0? А Ridge? Ситается, что регуляризатор исключает признак, если коэффициент при нем < 1e-3.
    • Ответ: Да, можна утверждать, что Lasso исключает один из признаков windspeed при любом значении alpha > 0. Нет, для Ridge такое утверждение не верно.
  4. Какой из регуляризаторов подойдет для отбора неинформативных признаков?
    • Ответ: Lasso

Далее будем работать с 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 [19]:
from sklearn.linear_model import LassoCV

In [27]:
# Код 3.2 (1 балл)
# Обучите регрессор LassoCV на всех параметрах регуляризации из alpha
# Постройте график _усредненного_ по строкам MSE в зависимости от alpha. 
# Выведите выбранное alpha, а также пары "признак-коэффициент" для обученного вектора коэффициентов
alphas = np.arange(1, 100, 5)
las = LassoCV(alphas=alphas, cv=3)
las.fit(X, y)
avg = np.average(las.mse_path_, axis = 1)
indexes = np.argsort(las.alphas_)
plt.plot(alphas, avg[indexes])
plt.title('Зависимость MSE от alpha')
plt.xlabel('Alpha')
plt.ylabel('MSE')

print('Minimum MSE: %s\nAlpha: %s' % (np.min(avg), alphas[indexes[np.argmin(avg)]]))


Minimum MSE: 790397.72811
Alpha: 6

Итак, мы выбрали некоторый параметр регуляризации. Давайте посмотрим, какие бы мы выбирали alpha, если бы делили выборку только один раз на обучающую и тестовую, то есть рассмотрим траектории MSE, соответствующие отдельным блокам выборки.


In [32]:
# Код 3.3 (1 балл)
# Выведите значения alpha, соответствующие минимумам MSE на каждом разбиении (то есть по столбцам).
# На трех отдельных графиках визуализируйте столбцы .mse_path_
def build_by_column_index(las, index):
    plt.figure(figsize=(8, 5))
    indexes = np.argsort(las.alphas_)
    avg = las.mse_path_[:, index]
    plt.plot(alphas, avg[indexes])
    plt.title('Зависимость MSE от alpha')
    plt.xlabel('Alpha')
    plt.ylabel('MSE')

    print('Minimum MSE: %s\nAlpha: %s' % (np.min(avg), alphas[indexes[np.argmin(avg)]]))
    
build_by_column_index(las, 0)
build_by_column_index(las, 1)
build_by_column_index(las, 2)


Minimum MSE: 843336.181499
Alpha: 41
Minimum MSE: 772598.495628
Alpha: 6
Minimum MSE: 745668.605961
Alpha: 1

In [35]:
las = Lasso(alpha=6)
las.fit(X, y)
for weight, name in sorted(zip(las.coef_, df.columns[:-1])):
    print('%s - %s' % (name, weight))


weathersit - -330.55985674
windspeed(mph) - -191.771408471
hum - -140.007405501
mnth - -100.039526144
holiday - -83.2939598753
windspeed(ms) - -2.68542954683e-08
workingday - 51.5570856141
weekday - 132.504465491
temp - 370.67985503
season - 532.018982841
atemp - 581.396931065
yr - 1015.06022264

На каждом разбиении оптимальное значение alpha свое, и ему соответствует большое MSE на других разбиениях. Получается, что мы настраиваемся на конкретные обучающие и контрольные выборки. При выборе alpha на кросс-валидации мы выбираем нечто "среднее", что будет давать приемлемое значение метрики на разных разбиениях выборки.

Наконец, как принято в анализе данных, давайте проинтерпретируем результат.

Блок 3. Ответьте на вопросы (каждый 0.5 балла):

  1. В последней обученной модели выберите 4 признака с наибольшими (положительными) коэфициентами (и выпишите их), посмотрите на визуализации зависимостей cnt от этих признаков, которые мы рисовали в блоке "Знакомство с данными". Видна ли возрастающая линейная зависимость cnt от этих признаков по графикам? Логично ли утверждать (из здравого смысла), что чем больше значение этих признаков, тем больше людей захотят взять велосипеды?
    • Ответ: yr, atemp, season, temp. Да, для этих признаков видна линейная зависимость для cnt от них. Да, логично утвреждать что чем больше значение этих признаков, тем больше людей захотят взять велосипеды.
  2. Выберите 3 признака с наибольшими по модулю отрицательными коэффициентами (и выпишите их), посмотрите на соответствующие визуализации. Видна ли убывающая линейная зависимость? Логично ли утверждать, что чем больше величина этих признаков, тем меньше людей захотят взять велосипеды?
    • Ответ: weathersit, windspeed(mph), hum, mnth. Да, логично утвреждать что чем менше значение этих признаков, тем больше людей захотят взять велосипеды.
  3. Выпишите признаки с коэффициентами, близкими к нулю (< 1e-3). Как вы думаете, почему модель исключила их из модели (вновь посмотрите на графики)? Верно ли, что они никак не влияют на спрос на велосипеды?
    • Ответ: windspeed(ms). Нет, данный признак влияет на спрос на велосипеды, просто он уже был учтен в параметру windspeed(mph) и тепер является избыточным для модели.

Заключение

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

Стоит отметить, что с помощью кросс-валидации удобно подбирать лишь небольшое число параметров (1, 2, максимум 3), потому что для каждой допустимой их комбинации нам приходится несколько раз обучать модель, а это времязатратный процесс, особенно если нужно обучаться на больших объемах данных.