Множественная линейная регрессия. Часть 2. | Оптимизация

Линейная регрессия. Часть 2.

Все курсы > Оптимизация > Занятие 4 (часть 2)

Во второй части занятия перейдем к практике.

Содержание занятия

Продолжим работать в том же ноутбуке

Сквозной пример

Данные и постановка задачи

Обратимся к хорошо знакомому нам датасету недвижимости в Бостоне.

При этом нам нужно будет решить две основные задачи:

Задача 1. Научиться оценивать качество модели не только с точки зрения метрики, но и исходя из рассмотренных ранее допущений модели. Эту задачу мы решим в три этапа.

  • Этап 1. Построим базовую (baseline) модель линейной регрессии с помощью класса LinearRegression библиотеки sklearn и оценим, насколько выполняются рассмотренные выше допущения.
  • Этап 2. Попробуем изменить данные таким образом, чтобы модель в большей степени соответствовала этим критериям.
  • Этап 3. Обучим еще одну модель и посмотрим, как изменится результат.

Задача 2. С нуля построить модель множественной линейной регрессии и сравнить прогноз с результатом, полученным при решении первой задачи. При этом обучение модели мы реализуем двумя способами, а именно, через:

  • Метод наименьших квадратов
  • Метод градиентного спуска

Разделение выборки

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

Исследовательский анализ данных

Теперь давайте более внимательно посмотрим на имеющиеся у нас данные. Как вы вероятно заметили, все признаки в этом датасете количественные, за исключением переменной CHAS.

Посмотрим на распределение признаков с помощью boxplots.

boxplots количественных признаков

Посмотрим на распределение целевой переменной.

распределение целевой переменной

Посмотрим на корреляцию количественных признаков с целевой переменной.

корреляция количественных признаков с целевой переменной

Используем точечно-бисериальную корреляцию для оценки взамосвязи переменной CHAS и целевой переменной.

Обработка данных

Пропущенные значения

Посмотрим, есть ли пропущенные значения.

Выбросы

Удалим выбросы.

При удалении выбросов важно помнить, что полное отсутствие вариантивности в данных не позволит выявить взаимосвязи.

Масштабирование признаков

Приведем признаки к одному масштабу (целевую переменную трогать не будем).

Замечу, что метод наименьших квадратов не требует масштабирования признаков, градиентному спуску же напротив необходимо, чтобы все значения находились в одном диапазоне (подробнее в дополнительных материалах).

Кодирование категориальных переменных

После стандартизации переменная CHAS также сохранила только два значения.

Ее можно не трогать.

Построение модели

Создадим первую пробную (baseline) модель с помощью библиотеки sklearn.

baseline-модель

Оценка качества

Диагностика модели, метрики качества и функции потерь

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

Функция потерь и метрика качества могут совпадать, а могут и не совпадать. Важно понимать, что у них разное назначение.

  • Функция потерь — это часть алгоритма, нам важно, чтобы эта функция была дифференцируема (у нее была производная)
  • Производная метрики качества нас не интересует. Метрика качества должна быть адекватна решаемой задаче.

MSE, RMSE, MAE, MAPE

MSE и RMSE

Для оценки качества RMSE предпочтительнее, чем MSE, потому что показывает насколько ошибается модель в тех же единицах измерения, что и целевая переменная. Например, если диапазон целевой переменной от 80 до 100, а RMSE 20, то в среднем модель ошибается на 20-25 процентов.

В качестве практики напишем собственную функцию.

Сравним с sklearn.

MAE

Приведем формулу.

$$ MAE = \frac{\sum |y-\hat{y}|}{n} $$

Средняя абсолютная ошибка представляет собой среднее арифметическое абсолютной ошибки $\varepsilon = |y-\hat{y}| $ и использует те же единицы измерения, что и целевая переменная.

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

MAPE

Средняя абсолютная ошибка в процентах (mean absolute percentage error) по сути выражает MAE в процентах, а не в абсолютных величинах, представляя отклонение как долю от истинных ответов.

$$ MAPE = \frac{1}{n} \sum \vert \frac{y-\hat{y}}{y} \vert $$

Это позволяет сравнивать модели с разными единицами измерения между собой.

Коэффициент детерминации

В рамках вводного курса в ответах на вопросы к занятию по регрессии мы подробно рассмотрели коэффициент детерминации ($R^2$), его связь с RMSE, а также зачем нужен скорректированный $R^2$. Как мы знаем, если использовать, например, класс LinearRegression, то эта метрика содержится в методе .score().

Также можно использовать функцию r2_score() модуля metrics.

Для скорректированного $R^2$ напишем собственную функцию.

Диагностика модели

Теперь проведем диагностику модели в соответствии со сделанными ранее допущениями.

Анализ остатков и прогнозных значений

Напишем диагностическую функцию, которая сразу выведет несколько интересующих нас графиков и метрик, касающихся остатков и прогнозных значений.

диагностика модели

Разберем полученную информацию.

  • В целом остатки модели распределены нормально с нулевым средним значением
  • Явной гетероскедастичности нет, хотя мы видим, что дисперсия не всегда равномерна
  • Присутствует умеренная отрицательная корреляция
  • График y_true vs. y_pred показывает насколько сильно прогнозные значения отклоняются от фактических. В идеальной модели (без шума, т.е. без случайных колебаний) точки должны были би стремиться находиться на диагонали, в более реалистичной модели нам бы хотелось видеть, что точки плотно сосредоточены вокруг диагонали.
  • Распределение прогнозных значений в целом повторяет распределение фактических.

Мультиколлинеарность

Отдельно проведем анализ на мультиколлинеарность. Напишем соответствующую функцию.

оценка мультиколлинеарности

Дополнительная обработка данных

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

VIF

Уберем признак с наибольшим VIF (RAD) и посмотрим, что получится.

удаление признаков с высоким VIF

Показатели пришли в норму. Окончательно удалим RAD.

Преобразование данных

Применим преобразование Йео-Джонсона.

Отбор признаков

Посмотрим на линейную корреляцию Пирсона количественных признаков и целевой переменной.

корреляция признаков и целевой переменной Пирсона после обработки данных

Также рассчитаем точечно-бисериальную корреляцию.

Удалим признаки с наименьшей корреляцией, а именно ZN, CHAS, DIS и B.

Повторное моделирование и диагностика

Повторное моделирование

Выполним повторное моделирование.

Оценка качества и диагностика

Оценим качество. Так как мы преобразовали целевую переменную, показатель RMSE не будет репрезентативен. Воспользуемся MAPE и $R^2$.

Отклонение прогнозного значения от истинного снизилось. $R^2$ немного уменьшился, чтобы бывает, когда мы пытаемся привести модель к соответствию допущениям. Проведем диагностику.

диагностика модели после обработки данных

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

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

Коэффициенты

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

Обучение модели

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

Векторизация уравнения

Для удобства векторизуем приведенное выше уравнение множественной линейной регрессии

$$ y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} X = \begin{bmatrix} x_0 & x_1 & \ldots & x_j \\ x_0 & x_1 & \ldots & x_j \\ \vdots & \vdots & \vdots & \vdots \\ x_{0} & x_{1} & \ldots & x_{n,j} \end{bmatrix}, \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix}, \varepsilon = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} $$

где n — количество наблюдений, а j — количество признаков.

Обратите внимание, что мы создали еще один столбец данных $ x_0 $, который будем умножать на сдвиг $ \theta_0 $. Его мы заполним единицами.

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

$$ y = X\theta + \varepsilon $$

Кроме того, как мы увидим ниже, так нам не придется искать отдельную производную для коэффициента $ \theta_0 $.

Схематично для модели с четырьмя наблюдениями (n = 4) и двумя признаками (j = 2) получаются вот такие матрицы.

умножение матрицы на вектор в модели линейной регрессии

Функция потерь

Как мы уже говорили, чтобы подобрать оптимальные коэффициенты $\theta$, нам нужен критерий или функция потерь. Логично измерять отклонение прогнозного значения от истинного.

$$ \varepsilon = X\theta-y $$

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

Абсолютная ошибка, L1 loss

При усреднении на количество наблюдений мы получаем среднюю абсолютную ошибку (mean absolute error, MAE).

$$ MAE = \frac{\sum{|y-X\theta|}}{n} = \frac{\sum{|\varepsilon|}}{n} $$

Приведем пример такой функции на Питоне.

Помимо модуля ошибку можно возводить в квадрат.

Квадрат ошибки, L2 loss

В этом случай говорят про сумму квадратов ошибок (sum of squared errors, SSE) или сумму квадратов остатков (sum of squared residuals, SSR или residual sum of squares, RSS).

$$ SSE = \sum (y-X\theta)^2 $$

Как мы уже знаем, на практике вместо SSE часто используется MSE, или вернее half MSE для удобства нахождения производной.

$$ MSE = \frac{1}{2n} \sum (y-\theta X)^2 $$

Ниже код на Питоне.

На практике у обеих функций есть сильные и слабые стороны. Рассмотрим L1 loss (MAE) и L2 loss (MSE) на графике.

сравнение MSE и MAE

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

Функция L1 не дает такой большой ошибки на выбросах, однако ее сложно дифференцировать, в точке минимума ее производная не определена.

Функция Хьюбера

Рассмотрим функцию Хьюбера (Huber loss), которая объединяет сильные стороны вышеупомянутых функций и при этом лишена их недостатков. Посмотрим на формулу.

$$ L_{\delta}= \left\{\begin{matrix} \frac{1}{2}(y-\hat{y})^{2} & if | y-\hat{y} | < \delta\\ \delta (|y-\hat{y}|-\frac1 2 \delta) & otherwise \end{matrix}\right. $$

Представим ее на графике.

функция Хьюбера

Также приведем код этой функции.

На сегодняшнем занятии мы, как и раньше, в качестве функции потерь используем MSE.

Метод наименьших квадратов

Нормальные уравнения

Для множественной линейной регрессии коэффициенты находятся по следующей формуле

$$ \theta = (X^TX)^{-1}X^Ty $$

Разберем эту формулу подробнее. Сумму квадратов остатков (SSE) можно переписать как произведение вектора $ \hat{\varepsilon} $ на самого себя, то есть $ SSE = \varepsilon^{T} \varepsilon$. Помня, что $\varepsilon = y-X\theta $ получаем (не забывая транспонировать)

$$ (y-X\theta)^T(y-X\theta) $$

Раскрываем скобки

$$ y^Ty-y^T(X\theta)-(X\theta)^Ty+(X\theta)^T(X\theta) $$

Заметим, что $A^TB = B^TA$, тогда

$$ y^Ty-(X\theta)^Ty-(X\theta)^Ty+(X\theta)^T(X\theta)$$

$$ y^Ty-2(X\theta)^Ty+(X\theta)^T(X\theta) $$

Вспомним, что $(AB)^T = A^TB^T$, тогда

$$ y^Ty-2\theta^TX^Ty+\theta^TX^TX\theta $$

Теперь нужно найти частные производные этих функций

$$ \nabla_{\theta} J(\theta) = y^Ty-2\theta^TX^Ty+\theta^TX^TX\theta $$

После дифференцирования мы получаем следующую производную

$$ -2X^Ty+2X^TX\theta $$

Как мы помним, оптимум функции находится там, где производная равна нулю.

$$ -2X^Ty+2X^TX\theta = 0 $$

$$ -X^Ty+X^TX\theta = 0 $$

$$ X^TX\theta = X^Ty $$

Выражение выше называется нормальным уравнением (normal equation). Решив его для $\theta$ мы найдем аналитическое решение минимизации суммы квадратов отклонений.

$$ \theta = (X^TX)^{-1}X^Ty $$

По теореме Гаусса-Маркова, оценка через МНК является наиболее оптимальной (обладающей наименьшей дисперсией) среди всех методов построения модели.

Код на Питоне

Перейдем к созданию класса линейной регрессии наподобие LinearRegression библиотеки sklearn. Вначале напишем функцию гипотезы (т.е. функцию самой модели), снабдив ее еще одной функцией, которая добавляет столбец из единиц к признакам.

$$ h_{\theta}(x) = \theta X $$

Перейдем к функции, отвечающей за обучение модели.

$$ \theta = (X^TX)^{-1}X^Ty $$

Обучим модель и выведем коэффициенты.

Примечание. Замечу, что не все матрицы обратимы. В этом случае можно найти псевдообратную матрицу (pseudoinverse). Для этого в Numpy есть функция np.linalg.pinv().

Сделаем прогноз.

Создание класса

Объединим код в класс.

Создадим объект класса и обучим модель.

Выведем коэффициенты.

Сделаем прогноз.

Оценка качества

Оценим качество через MAPE и $R^2$.

Мы видим, что результаты аналогичны.

Метод градиентного спуска

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

Нахождение градиента

Покажем расчет градиента на схеме.

нахождение градиента в модели линейной регрессии

В данном случае мы берем датасет из четырех наблюдений и двух признаков ($x_1$ и $x_2$) и соответственно используем три коэффициента ($\theta_0, \theta_1, \theta_2$).

Пошаговое построение модели

Начнем с функции гипотезы.

$$ h_{\theta}(x) = \theta X $$

Объявим функцию потерь.

$$ J({\theta_j}) = \frac{1}{2n} \sum (y-\theta X)^2 $$

Объявим функцию для расчета градиента.

$$ \frac{\partial}{\partial \theta_j} J(\theta) = -x_j(y — X\theta) \times \frac{1}{n} $$

где j — индекс признака.

Напишем функцию для обучения модели.

$$ \theta_j := \theta_j-\alpha \frac{\partial}{\partial \theta_j} J(\theta) $$

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

Обучим модель, выведем коэффициенты и достигнутый (минимальный) уровень ошибки.

Полученный результат очень близок к тому, что было найдено методом наименьших квадратов.

Прогноз

Сделаем прогноз.

Создание класса

Объединим написанные функции в класс.

Создадим объект класса, обучим модель, выведем коэффициенты и сделаем прогноз.

Оценка качества

Теперь рассмотрим несколько дополнительных соображений, касающихся построения модели линейной регрессии.

Диагностика алгоритма

Работу алгоритма можно проверить с помощью кривой обучения (learning curve).

  • Ошибка постоянно снижается
  • Алгоритм остановится, после истечения заданного количества итераций
  • Можно задать пороговое значение, после которого он остановится (например, $10^{-1}$)

Построим кривую обучения.

зависимость ошибки от итераций
зависимость ошибки от итераций (первые 100 итераций)

Она также позволяет выбрать адекватный коэффициент скорости обучения.

Подведем итог

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

Кроме того, мы узнали как изнутри устроены метод наименьших квадратов и метод градиентного спуска и построили соответствующие модели на Питоне.

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

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

Дополнительные материалы к занятию.