SecondBrain/90 Library/Machine Learning/Линейные модели.md

status, type, tags, created, updated, aliases

status	type	tags	created	updated	aliases
processing	concept	machine-learning linear-models regression classification	2026-02-26	2026-05-07	Линейные модели

Линейные модели

Общее представление о линейных моделях

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

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

Подробнее про линейную разделимость поговорим позже.

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

• Линейная регрессия используется для прогнозирования непрерывных значений, то есть для решения задачи регрессии.

• Логистическая регрессия используется для классификации

Этот урок мы посвятим рассмотрению модели линейной регрессии, а в следующем поговорим о логистической.

Линейная Регрессия

Постановка Задачи

Будем рассматривать задачу регрессии:

Пусть дана выборка пар объектов и ответов к ним:

Q = \{(x_i, y_i)\}_{i=1}^n =\{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}

где

• x_i = (x_{i1}, x_{i2}, ..., x_{im}), x_i ∈ X, X \subset \mathbb{R} ^{m} - множество объектов
• y_i ∈ Y, Y \subset \mathbb{R} - множество целевой переменной
• m - количество признаков
• n - размер выборки

Необходимо построить модель, восстанавливающую зависимость y от x:

y=f(x)

Общая Идея

Если исходить из предположения, что зависимость между y и x является линейной, то в качестве модели f(x) можно использовать линейную функцию вида:

f(x) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_m x_m

То есть прогнозы модели \hat{y} будут строиться на основе этой функции:

\hat{y} = f(x)= \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_m x_{m} = \beta_0 + \sum_{j=1}^m \beta_j \cdot x_{j}

где

• {x} = (x_{1}, x_{2}, ..., x_{m}) - признаковое описание объекта
• {y} - значение таргета для объекта
• \hat{y} - прогноз модели для объекта
• f - функция модели, принимающая на вход объект и возвращающая прогноз
• \beta_0, \beta_1, \beta_2, ..., \beta_m - коэффициенты/веса/параметры линейной регрессии, коэффициент \beta_0 часто еще называется свободным членом или смещением (bias) (не путать с bias из теории про смещение и разброс модели, это другое)

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

Рассмотрим эту формулировку на конкретном примере датасета.

import pandas as pd
import numpy as np
import plotly.express as px

data = pd.read_csv('https://raw.githubusercontent.com/merion-networks/data-science-course/refs/heads/main/unconv.csv')
data.head()

X = data.drop('Prod', axis=1)
y = data['Prod']

X.head()

y.head()

Признаки в данных:

В нашем случае имеется 1 таргет и 7 фичей. Обозначим их как:

• x_{1} - Well
• x_{2} - Por
• x_{3} - Perm
• x_{4} - AI
• x_{5} - Brittle
• x_{6} - TOC
• x_{7} - VR
• y - Prod

Тогда наша модель линейной регрессии будет иметь вид:

\hat{y} = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + \beta_3 x_{3} + \beta_4 x_{4} + \beta_5 x_{5} + \beta_6 x_{6} + \beta_7 x_{7}

Например, для i=1 (строка с индексом 0) получим:

\hat{y_1} = \beta_0 + \beta_1 \cdot 1 + \beta_2 \cdot 12.08 + \beta_3 \cdot 2.92 + \beta_4 \cdot 2.8 + \beta_5 \cdot 81.40 + \beta_6 \cdot 1.16 + \beta_7 \cdot 2.31

Для i=2 (строка с индексом 1) получим:

\hat{y_2} = \beta_0 + \beta_1 \cdot 1 + \beta_2 \cdot 12.38 + \beta_3 \cdot 3.53 + \beta_4 \cdot 3.22 + \beta_5 \cdot 46.17 + \beta_6 \cdot 0.89 + \beta_7 \cdot 1.88

и так далее...

Цель - найти такие коэффициенты \beta_i, чтобы ответы модели \hat{y_i} были как можно ближе к истинным ответам y_i.

О том как их искать и что значит "как можно ближе" поговорим далее.

Геометрическая Интерпретация

Рассмотрим частные случаи линейной регрессии и поговорим о ее геометрическом смысле.

2D вариант

Строим модель с одним признаком - по известной пористости скважин предсказывать неизвестную выработку газа.

Зависимость целевого признака от фактора представлена на диаграмме рассеяния.

Уравнение модели линейной регрессии будет иметь вид:

\hat{y} = \beta_0 + \beta_1 x

Представим, что мы нашли коэффициенты \beta_0 и \beta_1. Пусть коэффициенты составляют:

\beta_0 = -2.94, \beta_1=287.7

Тогда уравнение примет вид:

\hat{y} = -2.94 + 287.7 \cdot x

Геометрически данная запись задает уравнение прямой, которая пересекает ось ординат (y) в точке -2.94 (при x=0), а тангенс наклона прямой относительно оси x равен 287.7 (tg\alpha = 287.7).

Если подставлять значения конкретные значения пористости в модель, можно построить прямую, которая описывает исходную зависимость:

beta_0 = -2.94
beta_1 = 287

y_hat = beta_0 + beta_1 * X['Por']

y_hat

3D вариант

Теперь представим, что у нас не один фактор, а два. Например, помимо пористости скважины, мы дополнительно знаем ещё и о её хрупкости в процентах. То есть у нас теперь есть два фактора: x_1 — пористость и x_2 — хрупкость.

Можно отобразить зависимость добычи газа от этих факторов в трёхмерном пространстве в виде диаграммы рассеяния:

В таком случае в выражение для модели добавится ещё одна переменная и соответствующий ей коэффициент:

\hat{y} = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2}

Опять же представим, что мы нашли подходящие значения коэффициентов \beta_i, например:

\hat{y} = -2000 + 302.3 x_{1} + 31.8 x_{2}

Геометрически данное уравнение описывает плоскость в трёхмерном пространстве с осями x_1 и x_2, коэффициент \beta_0 — смещение плоскости по вертикальной оси, а коэффициенты \beta_1 и \beta_2 — коэффициенты (тангенсы угла) наклона этой плоскости к осям x_1 и x_2.

beta_0 = -2000
beta_1 = 302.3
beta_2 = 31.8

y_hat = beta_0 + beta_1 * X['Por'] + beta_2 * X['Brittle']
y_hat

MD вариант

А что если факторов не два, а больше: 3, 15, 100?

Тут-то мы и приходим к общему виду модели линейной регрессии, который вводили в самом начале:

\hat{y} = \beta_0 + \sum_{j=1}^m \beta_j \cdot x_{j}

В геометрическом смысле данное уравнение описывает плоскость в $(m+1)$-мерном пространстве (m факторов + 1 целевой признак отложены по осям координат). Такую плоскость называют гиперплоскостью.

Векторно-матричная Запись Модели

Запись модели линейной регрессии f(x):

\hat{y} = f(x) = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_m x_{m} = \beta_0 + \sum_{j=1}^m \beta_j \cdot x_{j}

можно значительно упростить, если обозначить коэффициенты и признаки как вектора:

• x = (x_{1}, x_{2}, ..., x_{m}) - вектор-описание $i$-ого объекта
• \beta = (\beta_1, \beta_2, ..., \beta_m)^T - вектор коэффициентов

То запись можно свести к виду скалярного произведения векторов:

\hat{y} = \beta_0 + \beta \cdot x

Можно пойти дальше и еще больше сократить запись, добавив \beta_0 в вектор коэффициентов, а в вектор объекта добавить 1:

• x = (1, x_{1}, x_{2}, ..., x_{m}) - вектор-описание объекта
• \beta = (\beta_0, \beta_1, \beta_2, ..., \beta_m)^T - вектор коэффициентов

Тогда:

\hat{y} = \beta \cdot x = \beta_0 \cdot 1 + \beta_1 \cdot x_{1} + \beta_2 \cdot x_{2} + ... + \beta_m \cdot x_{m}

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

Если вспомнить, что в датасете n объектов, то можно расписать уравнение регрессии для всего датасета,получим систему:

$$\begin{equation} \begin{cases} \hat{y_1} = \beta_0 + \beta_1 x_{11} + \beta_2 x_{12} + ... + \beta_m x_{1m} & \ \hat{y_2}= \beta_0 + \beta_1 x_{21} + \beta_2 x_{22} + ... + \beta_m x_{2m} & \ \cdots \ \hat{y_n} = \beta_0 + \beta_1 x_{n1} + \beta_2 x_{n2} + ... + \beta_m x_{nm} & \end{cases} \end{equation}$$

Если обозначить: $$X=\begin{pmatrix} 1 & x_{11} & x_{12} & ... & x_{1m} \ 1 & x_{21} & x_{22} & ... & x_{2m} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_{n1} & x_{n2} & ... & x_{nm} \end{pmatrix}$$

\hat{y}=\begin{pmatrix} \hat{y_1}\\ \hat{y_2}\\ ...\\ \hat{y_n} \end{pmatrix} \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_m\end{pmatrix}

То исходную запись можно перезаписать в матричном виде:

\hat{y} = X \cdot {\beta}

где:

• X - матрица признаков размера n\times (m+1), где n - количество объектов, а m - количество признаков
• \beta - вектор коэффициентов размера (m+1)\times 1
• \hat{y} - вектор предсказаний модели размера n\times 1

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

Пример:

# Создаём Вектор Из Единиц
ones = np.ones(X.shape[0])
# Добавляем Вектор К Таблице Первым Столбцом
X_ = np.column_stack([ones, X])

X_.round(1)

# Вектор-столбец Коэффициентов
beta = np.array([[-1232], [0], [230], [116], [-365], [25], [-78], [785]])

print(X_.shape)
print(beta.shape)

# Строим Прогноз
y_hat = X_ @ beta
# y_hat

Метод Наименьших Квадратов (OLS). MSE

Ключевой вопрос на который нам нужно ответить - как найти наилучшие коэффициенты \beta для нашей модели?

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

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

Говоря простым языком, мы должны научиться измерять качество модели и минимизировать её ошибку, как-то меняя обучаемые параметры. В случае линейной модели обучаемые параметры — это веса \beta.

Функция, показывающая как часто модель ошибается называется функцией потерь (функцией ошибки) или просто лоссом (loss function) и мы будем ее обозначать как L.

В качестве ошибки мы можем просто брать разницу между истинными ответами y и предсказаниями модели \hat{y} была минимальна. Чтобы не учитывать знак расхождения, можно взять модуль разницы между истинным значением и предсказанным

e_i = |y_i-\hat{y_i}|

Иллюстрация для понимания:

Однако, тут выясняется, что у модуля производная существует не везде (функция не дифференцируема в 0). Поэтому вместо модуля используют квадрат.

e_i = (y_i-\hat{y_i})^2

Ошибку на одном примере рассматривать смысла нет, минимизировать надо сумму ошибок по всей обучающей выборке. Тогда мы получим функцию потерь под названием суммарная квадратичная ошибка (Sum Squared Error, SSE):

L(y, \hat{y}) = SSE = \sum_{i=1}^n e_i^2 = \sum_{i=1}^n (y_i-\hat{y_i})^2

Часто на практике используют не просто сумму ошибок, а среднюю ошибку по всему датасету, хотя на результаты сходимости алгоритмов это никак не влияет. В таком случае функции потерь называются средняя квадратичная ошибка (Mean Squared Error, MSE):

L(y, \hat{y}) = MSE = \frac{1}{n} \sum_{i=1}^n e_i^2 = \frac{1}{n} \sum_{i=1}^n (y_i-\hat{y_i})^2

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

L(y, \hat{y}) = \frac{1}{n} ||y - \hat{y}||^2

К тому же \hat{y} - это функция, зависящая от X и \beta:

L(X, y, \beta) = ||y - X \cdot {\beta}||^2

Тогда задача сводится к оптимизации функционала:

L(X, y, \beta)= \frac{1}{n} ||y - X \cdot {\beta}||^2 \rightarrow \min_{\beta}

Запись L(X, y, \beta) означает, что функция потерь зависит от выборки (X, y) и параметров модели. Однако, выборка то у нас одна и та же. Поэтому в целях экономии места и упрощения формул, будем писать, что функция потерь зависит только от параметров модели L(\beta).

Оптимальными будут считаться такие параметры \beta_{opt}, которые дают минимум функции потерь:

\beta_{opt} = \arg \min_{\beta} L(\beta)

Рассмотрим примеры: возьмем 3 разных варианта коэффициентов \beta и рассчитаем MSE

def mse(y, y_hat):
    n = len(y)
    return ((y - y_hat) ** 2).sum() * 1 / n

# Варианты Коэффициентов
beta_1 = np.array([[-1232], [0], [230], [116], [-365], [25], [-78], [785]])
beta_2 = np.array([[-849], [50], [155], [67], [36], [0.17], [-89], [800]])
beta_3 = np.array([[0], [10], [245], [352], [-843], [53], [-36], [56]])

# Прогнозы
y_hat_1 = X_ @ beta_1
y_hat_2 = X_ @ beta_2
y_hat_3 = X_ @ beta_3

# Считаем Mse
print(mse(y.values, y_hat_1))
print(mse(y.values, y_hat_2))
print(mse(y.values, y_hat_3))

Полным перебором коэффициенты искать нам не придется. Для этого есть математические методы поиска минимума функции.

Метод для поиска коэффициентов регрессии, основанный на минимизации суммы квадратов ошибок (отклонений) был изобретен Гауссом ещё в 1795 году и позднее был назван методом наименьших квадратов (МНК). В английской литературе часто можно встретить аббревиатуру OLS (Ordinary Least Squares).

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

→ Существует теорема Гаусса-Маркова, которая говорит о том, что, если выполнены все условия теоремы, МНК всегда находит единственно оптимальные параметры.

Примеры для понимания

В простейшей ситуации, когда в модели используется всего 1 параметр \beta и в данных есть только 1 признак мы получим:

\hat{y} = x \cdot \beta L(\beta) = \frac{1}{n} \sum_{i=1}^n (y_i- x_i\cdot \beta)^2

Эта функция задает параболу:

Если добавить в модель еще один коэффициент \beta_0, то мы получим:

\hat{y} = \beta_0 + \beta x L(\beta) = \frac{1}{n} \sum_{i=1}^n (y_i- \beta_0 - x_i \cdot \beta )^2

Эта функция задает парабаллоид:

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

Чтобы решить задачу МНК и найти оптимальные параметры \beta_{opt}, обеспечивающие минимум функции есть два способа:

1. Аналитический
2. Численный

OLS: Аналитическое Решение

Аналитическое решение поиска минимума функции строится по следующему принципу:

1. Найти производную функции
2. Прировнять эту производную к 0
3. Решить полученное уравнение и получить ответ - точку минимума

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

В общем случае функция потерь L(\beta) является многомерной функцией (функцией от нескольких переменных), а производная функции заменяется на градиент

Градиент функции — это вектор, который состоит из частных производных по параметрам функции.

$$\nabla_{\beta} L(\beta) = \begin{pmatrix} \frac{\partial L}{\partial \beta_0} & \frac{\partial L}{\partial \beta_1} & ... & \frac{\partial L}{\partial \beta_m} \end{pmatrix}^T$$

где

• L — функция потерь, зависящая от параметров модели.
• \nabla_{\beta} — символ набла — символьное сокращение градиента, нижний индекс подчеркивает, что производные вычисляются по компонентам вектора \beta
• \frac{\partial L}{\partial \beta_i} - частная производная функции L(\beta) по переменной \beta_i

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

Например, для квадратичной функции потерь (squared loss) градиент в простейшем случае - 1-ом признаке и двух коэффициентах \beta_0 и \beta_1 будет иметь вид:

L(\beta) = \frac{1}{n} \sum_{i=1}^n(y_i - \hat{y_i})^2= \frac{1}{n} \sum_{i=1}^n(y_i - \beta_0 - x_i \cdot \beta_1)^2 \nabla_{\beta} L(\beta) = (\frac{\partial L}{\partial \beta_0}, \frac{\partial L}{\partial \beta_1}) \frac{\partial L}{\partial \beta_0} = - \frac{2}{n} \sum_{i=1}^n (y_i - β_0 - x_i \cdot β_1 ) \frac{\partial L}{\partial \beta_1} = - \frac{2}{n} \sum_{i=1}^n x_i (y_i - β_0 - x_i \cdot β_1)

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

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

L(\beta) = \frac{1}{n} ||y - X \cdot {\beta}||^2 \nabla_{\beta} L(\beta) = \frac{2}{n} \cdot X^T(X \cdot \beta - y)

где $$X=\begin{pmatrix} 1 & x_{11} & x_{12} & ... & x_{1m} \ 1 & x_{21} & x_{22} & ... & x_{2m} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_{n1} & x_{n2} & ... & x_{nm} \end{pmatrix}$$

\beta = (\beta_0, \beta_1, \beta_2, ..., \beta_m)^T

Вывод этой формулы для градиента можно найти здесь.

Как правило множитель \frac{2}{n} обычно опускают, так как константы не влияют на саму точку минимума.

Теперь для поиска точки минимума функции потерь остается только решить матричное уравнение:

\nabla_{\beta} L(\beta) = X^T(X \cdot \beta - y) = 0

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

Итоговая формула будет иметь вид:

\beta_{opt}=(X^T \cdot X)^{-1} \cdot X^T \cdot y

Данная формула матричная формула была получена Гауссом, и она позволяет очень просто и легко путем перемножения матриц найти оптимальные коэффициенты \beta_{opt}.

Аналитическое Решение: Пример Реализации

Реализуем матричную формулу МНК через numpy.

Для начала вспомним, что для вычисления свободного члена \beta_0 необходимо добавить в таблицу столбец, полностью состоящий из единиц. Такой столбец можно создать с помощью знакомой нам функции ones() из библиотеки numpy, а присоединить его к таблице X поможет функция column_stack().

def linear_regression_fit(X, y):
    """
    Функция для вычисления коэффициентов линейной регрессии
    """
# Создаём Вектор Из Единиц
    ones = np.ones(X.shape[0])
# Добавляем Вектор К Таблице Первым Столбцом
    X = np.column_stack([ones, X])

# Вычисляем Обратную Матрицу Q
    Q = np.linalg.inv(X.T @ X)
# Вычисляем Вектор Коэффициентов
    beta = Q @ X.T @ y
    return beta

beta_opt = linear_regression_fit(X, y)

print('Коэффициенты регрессии:')
for i, beta_i in enumerate(beta_opt):
    print(f'i={i}: {beta_i.round(2)}')

print(beta_opt)

Чтобы вычислить прогноз выработки газа для новых скважин нужно лишь скалярно перемножить матрицу с описанием объектов на вычисленные коэффициенты регрессии:

def linear_regression_predict(X, beta):
    """
    Функция для построения прогноза
    """
# Создаём Вектор Из Единиц
    ones = np.ones(X.shape[0])
# Добавляем Вектор К Таблице Первым Столбцом
    X = np.column_stack([ones, X])
    return X @ beta

# Новая Скважина С Неизвестной Выработкой
x = np.array([[195.  ,  17.13,   5.99,   2.88,  51.45,   1.77,   2.28]])

# Делаем Прогноз
y_hat = linear_regression_predict(x, beta_opt)

y_hat

Можно вычислить прогноз для всех скважин и рассчитать итоговую MSE

# Делаем Прогноз
y_hat = linear_regression_predict(X, beta_opt)

n = X.shape[0]

mse = ((y - y_hat)**2).sum() * 1 / n

print('MSE: {:.2f}'.format(mse))

MSE хороша для оптимизации, но плоха для интерпретации. Для оценки качества лучше использовать коэффициент детерминации (R^2).

Напоминание:

R^2 = 1 - \frac{{\sum_{i=1}^n (y_i - \hat{y_i})^2}}{\sum_{i=1}^n (y_i - y_{mean})^2}

где

• {y_i} - истинное значение таргета для $i$-ого объекта

• \hat{y_i} - прогноз для $i$-ого объекта

• y_{mean} - среднее значение по таргету

# Делаем Прогноз
y_hat = linear_regression_predict(X, beta_opt)

r2 = 1 - ((y - y_hat)**2).sum() / ((y - y.mean())**2).sum()

print('R^2-score: {:.2f}'.format(r2))

Конечно же, никто не строит линейную регрессию руками, используя формулу МНК. Все дата-сайентисты пользуются библиотеками, такими как sklearn.

Аналитическое Решение: Реализация В Sklearn

Все линейные модели, которые мы будем рассматривать, реализованы в модуле linear_model библиотеки sklearn. Давайте импортируем этот модуль:

from sklearn import linear_model #линейные модели
from sklearn import metrics #метрики

В модуле находится класс LinearRegression, который реализует поиск коэффициентов линейной регрессии по МНК.

Для обучения модели необходимо просто вызвать метод fit(), передав в него матрицу наблюдений X и вектор правильных ответов y.

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

# Создаём Объект Класса LinearRegression
lr = linear_model.LinearRegression()
# Обучаем Модель — Ищем Параметры По МНК
lr.fit(X, y)

Чтобы получить свободный член \beta_0 нужно обратиться по атрибуту intercept_.

Вектор оставшихся параметров \beta_1, \beta_2, ..., \beta_m хранится в атрибуте coef_:

lr.coef_

lr.intercept_

# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X.columns, 'Coefficients': lr.coef_})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df =pd.DataFrame({'Features': ['Intercept'], 'Coefficients': lr.intercept_})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

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

Модель обучена. Теперь можно сделать прогноз. Для этого есть метод predict(). В него необходимо передать матрицу наблюдений, для которых нужно сделать предсказание.

# Новая Скважина С Неизвестной Выработкой
x = np.array([
    [195.  ,  17.13,   5.99,   2.88,  51.45,   1.77,   2.28]
])

# Делаем Прогноз
y_hat = lr.predict(x)

y_hat

Аналогично можно сделать предсказание для всей выборки. В результате получим вектор предсказаний модели \hat{y}.

# Делаем Прогноз Для Всей Выборки
y_hat = lr.predict(X)

y_hat.shape

Рассчитаем качество полученного решения по метрике MSE и R^2:

print('MSE: {:.2f}'.format(metrics.mean_squared_error(y, y_hat)))

print('R^2-score: {:.2f}'.format(metrics.r2_score(y, y_hat)))

При необходимости можно обучить модель без свободного члена \beta_0. Для этого нужно установить параметр fit_intercept как False:

# Создаём Объект Класса LinearRegression
lr = linear_model.LinearRegression(fit_intercept=False)
# Обучаем Модель — Ищем Параметры По МНК
lr.fit(X, y)

# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X.columns, 'Coefficients': lr.coef_})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df =pd.DataFrame({'Features': ['Intercept'], 'Coefficients': lr.intercept_})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

Преимущества И Недостатки

Из формулы

\beta=(X^T \cdot X)^{-1} \cdot X^T \cdot y = Q \cdot X^T \cdot y

можно заметить недостатки аналитического решения:

1. Сложность обращения матриц

   Аналитическое решение предполагает вычисление обратной матрицы Q=(X^T \cdot X)^{-1}
   Обращение матриц очень дорогая операция. Вычислительно обращать большие матрицы дело сложное. Если мы посмотрим повнимательнее, то увидим, что если матрица X имеет размерность n\times m, то матрица X^T - m \times n, а значит произведение X^T \cdot X имеет размерность m \times m. То есть сложность обращения зависит не от количества данных, а от количества признаков m. Причем сложность вычисления обратной матрицы - кубическая O(m^3).
   Что это означает на практике?
   Если в датасете 1000 фич, то для вычисления обратной матрицы понадобится 1000^3=1 000 000 000 операций!!!
   Отсюда делаем вывод: классическое аналитическое решение работает на "широких" датасетах несколько дольше.
   Вычисления обратной матрицы можно немного ускорить используя специальные алгоритмы и путем распараллеливания вычисления элементов матрицы, что и реализовано в sklearn.

2. Несуществование обратной матрицы

   Из линейной алгебры известно, что, если матрица X^T \cdot X является вырожденной (det(X^T \cdot X) = 0), то у нее нет обратной, а значит и решения не существует.
   Напомним, что полностью вырожденные матрицы мы получаем в случае, если в данных присутствуют линейно-зависимые признаки (один из признаков является линейной комбинацией других). Пример - общая площадь здания, жилая и нежилая площадь

house_data = pd.DataFrame(
    {
        'Общая площадь': [105, 40, 83, 54, 64, 43, 109],
        'Жилая площадь': [89, 35, 67, 48, 50, 36, 79],
        'Нежилая площадь': [16,  5, 16,  6, 14,  7, 30],
        'Цена (млн. руб)': [100.3 ,  13.6 ,  51.85,  34.85,  52.7 ,  13.6 ,  89.25],
        'Цена (млн. $)': [1.18, 0.16, 0.61, 0.41, 0.62, 0.16, 1.05],
        'Цена аренды (тыс. руб)': [80, 35, 55, 50, 45, 40, 73]
    }
)
X = house_data.drop('Цена аренды (тыс. руб)', axis=1)
y = house_data['Цена аренды (тыс. руб)']

house_data.corr()

Определитель матрицы X^T \cdot X:

np.linalg.det(X.T @ X)

Попробуем построить линейную регрессию через матричную формулу:

linear_regression_fit(X, y)

Ожидаемо получаем ошибку сингулярной (вырожденной) матрицы. Обратной к матрице X^T \cdot X не существует.

А что нам на это скажет sklearn?

model = linear_model.LinearRegression()

model.fit(X, y)
pd.DataFrame({'feature': X.columns, 'coef': model.coef_})

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

Для этого используется метод под названием сингулярное разложение (SVD). О нём мы будем говорить отдельно в главе о методах понижения размерности.

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

Если вы хотите понять как это работает уже сейчас, ознакомьтесь с этой статьёй.

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

Однако, открытым остаётся вопрос: можно ли доверять коэффициентам, полученным таким способом, и интерпретировать их?

3. Проблема мультиколлинеарности
   Пусть вырожденность матрицы и решаемая проблема. Однако, это не отменяет того факта, что в данных может присутствовать мультиколлинеарность.
   Мультиколлинерность приводит к тому, что определитель матрицы det(X^T \cdot X) начинает стремиться к 0, коэффициенты модели стремится к +- бесконечности (\beta → |∞|), что приводит к неустойчивости коэффициентов линейных моделей и катострофическим последствиям.
   Коэффициенты, полученные с помощью аналитического решения на данных, в которых присутствует высокая мультиколинеарность, являются неустойчивыми: при малейшем изменении в целевой переменной, мы будем получать абсолютно другие решения.
   Стоит отметить коэффициенты также теряют всяческую интерпретацию.

Смотрим пример. Добавим в данные еще одну квартирку:

house_data = pd.DataFrame(
    {
        'Общая площадь': [105, 40, 83, 54, 64, 43, 109, 83],
        'Жилая площадь': [89, 35, 67, 48, 50, 36, 79, 78],
        'Нежилая площадь': [16,  5, 16,  6, 14,  7, 30, 5],
        'Цена (млн. руб)': [100.3 ,  13.6 ,  51.85,  34.85,  52.7 ,  13.6 ,  89.25, 65],
        'Цена (млн. $)': [1.18, 0.16, 0.61, 0.41, 0.62, 0.16, 1.05, 0.76],
        'Цена аренды (тыс. руб)': [80, 35, 55, 50, 45, 40, 73, 53]
    }
)
X = house_data.drop('Цена аренды (тыс. руб)', axis=1)
y = house_data['Цена аренды (тыс. руб)']

model = linear_model.LinearRegression()

model.fit(X, y)
pd.DataFrame({'feature': X.columns, 'coef': model.coef_})

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

Именно в этом и заключается нестабильность решения на данных с мультиколлинерностью - на разных обучающих выборках мы будем получать абсолютно разные коэффициенты при сильно-коррелированных факторах

OLS: Численное Решение

Часть недостатков аналитического решения можно решить, если для поиска минимума функции потерь L(\beta) использовать численные методы. Самым популярным и наиболее эффективным является метод градиентного спуска.

Градиентный Спуск

Градиентный спуск (Gradient Descent) - это численный метод оптимизации, который заключается в итеративном обновлении оценок коэффициентов регрессии в направлении, которое уменьшает значение функции потерь L(\beta).

Метод основан на двух очень важных свойствах вектора градиента:

• Градиент — это вектор, который показывает направление наискорейшего роста функции.
  Вектор противоположный градиенту - \nabla называется вектором антиградиента, и он показывает в сторону наискорейшего убывания функции.

• Длина — это само значение скорости роста функции в точке. В точке минимума \beta^* длина вектора равна 0 (||\nabla L(\beta^*)||=0). Это свойство мы можем использовать в качестве критерия остановки нашего алгоритма.

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

Формально это записывается следующим образом:

\beta^{(k+1)} = \beta^{(k)} - \eta \nabla L(\beta^{(k)})

где

• \beta — это вектор параметров модели, координаты в пространстве, а индекс в круглых скобках сверху означает номер точки в пространстве.
• \eta - темп обучения

Запись \nabla L(\beta^{(k)}) означает, что градиент вычисляется в текущей точке под номером.

Темп обучения \eta — это основной параметр алгоритма. Управляя данным параметром (уменьшая и увеличивая его), мы управляем скоростью движения к точке минимума. Чем больше темп обучения, тем длиннее наши шаги и тем быстрее мы движемся, и наоборот.

Приведенную формулу можно рассматривать с физической точки зрения: новая координата \beta^{(k+1)} в пространстве параметров определяется как текущая координата \beta^{(k)} минус скорость роста в текущей точке \nabla L(\beta^{(k)}), помноженная на коэффициент «скольжения» \eta.

Формализация алгоритма:

1. Проинициализировать алгоритм.

   Задать начальные значения оценок коэффициентов регрессии \beta_0, \beta_1, ..., \beta_m. Например, можно инициализировать все параметры нулями или случайными значениями.
   Задать максимальное количество итераций K, чтобы ограничить алгоритм по времени выполненения.
   Задать границу сходимости ϵ.

2. Повторять K раз, k - номер итерации:

• Вычислить градиент функции потерь \nabla L(\beta) в текущей точке \beta^{(k)}:

$$\nabla L(\beta^{(k)}) = \begin{pmatrix} \frac{\partial L}{\partial \beta_0} & \frac{\partial L}{\partial \beta_1} & ... & \frac{\partial L}{\partial \beta_m} \end{pmatrix}^T_{\beta=\beta^{(k)}}$$

• Проверить условие остановки - равенство нулю длины вектора градиента:
  ||\nabla L(\beta^{(k)})|| = 0
  Если условие остановки выполняется - завершить вычисления и вернуть найденные значения \beta.
  На практике полного равенства градиента нулю достичь невозможно из-за численных вычислений, поэтому в качестве остановки задают минимальную границу \epsilon, ниже которой длина градиента считается достаточной, чтобы остановиться (например, 0.1, 0.01 или 0.001):

||\nabla L(\beta^{(k)})|| < \epsilon

• Обновить оценки коэффициентов регрессии по формуле:

\beta^{(k+1)} = \beta^{(k)} - \eta \nabla L(\beta^{(k)})

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

L(\beta) = \frac{1}{n} ||y - X \cdot {\beta}||^2 = \frac{1}{n} \sum_{i=1}^n(y_i - \beta \cdot x_i)^2 \nabla L(\beta) = \frac{2}{n} \cdot X^T(X \cdot \beta - y)

Иллюстрация работы градиентного спуска:

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

Градиентный Спуск: Пример Расчётов И Реализации *

В целях экономии ручных расчётов в реализации ниже опущен коэффициент \frac{2}{n}, так как он не влияет на направление сходимости: алгоритм всё равно движется к точке минимума, просто с другим масштабом шага. В реализации градиентного спуска в sklearn данный коэффициент учитывается при расчёте.

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

Пусть у нас есть датасет состоящий из 5 наблюдений:

№	x	y
1	1	2
2	2	4
3	3	5
4	4	7
5	5	8

Необходимо восстановить зависимость x от y с помощью модели линейной регрессии:

\hat{y_i} = f(x_i) = \beta_0 + \beta_1 x_i

Поиск коэффициентов \beta_0 и \beta_1 будем производить с помощью градиентного спуска.

Инициализируем параметры:

• максимальное количество итераций: K=100
• темп обучения: \eta = 0.01.
• граница сходимости: \epsilon = 0.1
• начальное значение вектора коэффициентов: \beta^{(0)} = \begin{pmatrix} β_0 \\ β_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}
• функция потерь: L(β) = \sum_{i=1}^n (y_i - \hat{y_i})^2 = \sum_{i=1}^n (y_i - β_0 - β_1 x_i)^2

Итерация 1 ($k=1)$

• Вычисляем градиент функции потерь при \beta^{(0)} = (0, 0)^T:

\nabla L(\beta^{(0)}) = \begin{pmatrix} \frac{\partial L}{\partial β_0} \\\ \frac{\partial L}{\partial β_1} \end{pmatrix} = \begin{pmatrix} - \sum_{i=1}^n (y_i - β_0 - β_1 x_i) \\\ - \sum_{i=1}^n x_i (y_i - β_0 - β_1 x_i) \end{pmatrix}= \begin{pmatrix} - \sum_{i=1}^5 (y_i - 0 - 0 \cdot x_i) \\\ - \sum_{i=1}^5 x_i (y_i - 0 - 0 \cdot x_i) \end{pmatrix} = \begin{pmatrix} - \sum_{i=1}^5 y_i \\\ - \sum_{i=1}^5 x_i y_i \end{pmatrix} = \begin{pmatrix} - (2 + 4 + 5 + 7 + 8) \\\ - (1 \cdot 2 + 2 \cdot 4 + 3 \cdot 5 + 4 \cdot 7 + 5 \cdot 8) \end{pmatrix}= \begin{pmatrix} - 26 \\\ - 93 \end{pmatrix}

• Проверяем условие остановки. Вычисляем длину вектора градиента:

||\nabla L(β^{(0)})|| = \sqrt{(\frac{\partial L}{\partial β_0})^2 + (\frac{\partial L}{\partial β_1})^2}=\sqrt{(-26)^2 + (-93)^2} \approx 96.6

Условие остановки не выполнено: ||\nabla L(β^{(0)})|| > \epsilon.

• Обновляем коэффициенты:

β^{(1)} = β^{(0)} - \eta \nabla L(β^{(0)}) = \begin{pmatrix} 0 \\\ 0 \end{pmatrix} - 0.01 \begin{pmatrix} - 26 \\\ - 93 \end{pmatrix} = \begin{pmatrix} 0.26 \\\ 0.93 \end{pmatrix}

Итерация 2 (k=2)

• Вычисляем градиент функции потерь при \beta^{(1)} = (0.26, 0.93)^T:

\nabla L(β^{(1)}) = \begin{pmatrix} \frac{\partial L}{\partial β_0} \\\ \frac{\partial L}{\partial β_1} \end{pmatrix} = \begin{pmatrix} - \sum_{i=1}^n (y_i - β_0 - β_1 x_i) \\\ - \sum_{i=1}^n x_i (y_i - β_0 - β_1 x_i) \end{pmatrix}= \begin{pmatrix} - \sum_{i=1}^5 (y_i - 0.26 - 0.93 \cdot x_i) \\\ - \sum_{i=1}^5 x_i (y_i - 0.26 - 0.93 \cdot x_i) \end{pmatrix}= \begin{pmatrix} - (2 - 0.26 - 0.93 \cdot 1 + 4 - 0.26 - 0.93 \cdot 2 + 5 - 0.26 - 0.93 \cdot 3 + 7 - 0.26 - 0.93 \cdot 4 + 8 - 0.26 - 0.93 \cdot 5) \\\ - (1 \cdot (2 - 0.26 - 0.93 \cdot 1) + 2 \cdot (4 - 0.26 - 0.93 \cdot 2) + 3 \cdot (5 - 0.26 - 0.93 \cdot 3) + 4 \cdot (7 - 0.26 - 0.93 \cdot 4) + 5 \cdot (8 - 0.26 - 0.93 \cdot 5)) \end{pmatrix}= \begin{pmatrix} -10.75 \\\ -37.95 \end{pmatrix}

• Проверяем условие остановки. Вычисляем длину вектора градиента:

||\nabla L(β^{(1)})|| = \sqrt{(\frac{\partial L}{\partial β_0})^2 + (\frac{\partial L}{\partial β_1})^2}=\sqrt{(-10.75)^2 + (-37.95)^2} \approx 39.44

Условие остановки не выполнено: ||\nabla L(β^{(0)})|| > \epsilon.

• Обновляем коэффициенты:

β^{(2)} = β^{(1)} - \eta \nabla L(β^{(1)})= \begin{pmatrix} 0.26 \\\ 0.93 \end{pmatrix} - 0.01 \begin{pmatrix} - 10.75 \\\ -37.95 \end{pmatrix}= \begin{pmatrix} 0.3675 \\\ 1.3095 \end{pmatrix}

Итерация 3 (k=3)

• Вычисляем градиент функции потерь при \beta^{(2)} = (0.3675, 1.3095)^T:

\nabla L(β^{(2)}) = \begin{pmatrix} \frac{\partial L}{\partial β_0} \\\ \frac{\partial L}{\partial β_1} \end{pmatrix} = \begin{pmatrix} - \sum_{i=1}^n (y_i - β_0 - β_1 x_i) \\\ - \sum_{i=1}^n x_i (y_i - β_0 - β_1 x_i) \end{pmatrix} = \begin{pmatrix} - \sum_{i=1}^5 (y_i - 0.3675 - 1.3095 \cdot x_i) \\\ - \sum_{i=1}^5 x_i (y_i - 0.3675 - 1.3095 \cdot x_i) \end{pmatrix} = \begin{pmatrix} - -4.52 \\\ -15.465 \end{pmatrix}

• Проверяем условие остановки. Вычисляем длину вектора градиента:

||\nabla L(β^{(1)})|| = \sqrt{(\frac{\partial L}{\partial β_0})^2 + (\frac{\partial L}{\partial β_1})^2}=\sqrt{(-4.52)^2 + (-15.465)^2} \approx 16.1

Условие остановки не выполнено: ||\nabla L(β^{(0)})|| > \epsilon.

• Обновляем коэффициенты:

β^{(3)} = β^{(2)} - \eta \nabla L(β^{(2)}) = \begin{pmatrix} 0.3675 \\\ 1.3095 \end{pmatrix} - 0.1 \begin{pmatrix} -4.52 \\\ -15.465 \end{pmatrix} = \begin{pmatrix} 0.4127 \\\ 1.46415 \end{pmatrix}

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

Например, уже на 10-ой итерации (k=10) мы с вами получим:

\nabla L(β^{(10)}) = \begin{pmatrix} -0.215 \\ 0.028 \end{pmatrix} ||\nabla L(β^{(10)})|| = 0.217 \beta = \begin{pmatrix} 0.457 \\ 1.567 \end{pmatrix}

А на 100-ой итерации наш алгоритм сойдется:

\nabla L(β^{(10)}) = \begin{pmatrix} -0.096 \\ 0.027 \end{pmatrix} ||\nabla L(β^{(10)})|| = 0.1 \beta = \begin{pmatrix} 0.587 \\ 1.531 \end{pmatrix}

В итоге мы получим значения для коэффициентов \beta:

\beta = \begin{pmatrix} 0.58692341 \\ 1.5313204 \end{pmatrix}

Таким образом, уравнение линейной регрессии будет иметь вид:

\hat{y_i} = f(x_i) = 0.58692341 + 1.5313204 \cdot x_i

Подставляя в уравнение точки x_i из обучающего набора данных мы получим следующие прогнозы:

№	x	y	\hat{y}
1	1	2	2.11824381
2	2	4	3.64956422
3	3	5	5.18088462
4	4	7	6.71220502
5	5	8	8.24352542

Пример ручной реализации алгоритма поиска коэффициентов с помощью градиентного спуска для примера выше:

import numpy as np


# Исходные Данные
X = np.array([[1], [2], [3], [4], [5]], dtype='float64')
y = np.array([2, 4, 5, 7, 8], dtype='float64')

# Здесь Должен Быть Шаг С Масштабированием
# Но Так Как У Нас Один Признак, То Масштабировать Его Не Нужно
# Scaler = preprocessing.StandardScaler()
# X = scaler.fit_transform(X)

# Добавляем Столбец Из Единиц В Матрицу Наблюдений Для Удобства Рассчета Градиентов
ones = np.ones(X.shape[0])
X = np.column_stack([ones, X])

# Инициализируем Параметры
max_iter = 1000  # Максимальное количество итераций
eta = 0.01  # Темп обучения
epsilon = 0.1  # Граница сходимости
beta = np.array([0, 0], dtype='float64')  # Начальное значение вектора коэффициентов

# Обучаем Модель
for k in range(max_iter):

# Вычисляем Градиент Функции Потерь
    gradient = X.T@ X @ beta - X.T @ y
# Проверяем Условие Остановки
    if np.linalg.norm(gradient) < epsilon:
        break
# Обновляем Коэффициенты
    beta = beta - eta * gradient
    print('Итерация:', k + 1)
    print('Градиент:',  gradient)
    print('Длина градиента:', np.linalg.norm(gradient))
    print('Коэффициенты: ', beta)

# Выводим Результаты
print("Коэффициенты линейной регрессии:")
print(beta)

# Делаем Прогнозы
y_pred = X @ beta

# Выводим Прогнозы
print("Прогнозы:")
print(y_pred)

Стохастический Градиентный Спуск

У классического градиентного спуска есть одна большая проблема — это сходимость алгоритма к точке истинного минимума. Если в случае аналитического решения все просто - подставил матрицы и вектора в формулу и гарантировано получил ответ, то в случае численных методов, алгоритм может попросту не сойтись к истинному минимуму.

Сходимость градиентного спуска зависит от многих факторов, главные из которых:

• сложности зависимости и сложности функции потерь;
• выбранный темп обучения;
• выбранная начальная точка (инициализация параметров);
• масштабирование признаков.

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

Примеры графиков функций, имеющих несколько минимумов (1 параметр - слева, 2 параметра - справа):

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

Проблема градиентного спуска заключается в том, что алгоритм может «застрять» в локальном минимуме и не выйти из него.

Другая проблема градиентного спуска - это количество вычислений, необходимых для вычисления градиента по всей выборке. То есть на каждом шаге мы прогоняем всю выборку через алгоритм. Если выборка будет очень большой (на миллионы записей), поиск оптимальных коэффициентов может занять немалое количество времени и оперативной памяти.

Вычислительная сложность всего алгоритма - O(N \cdot M \cdot K), где N - число объектов в датасете, M - число признаков, K - количество итераций. Да, это гораздо лучше, чем аналитическое решение через обратную матрицу, сложность которого O(N^2 \cdot M + M^3). Но можно и быстрее.

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

В контексте линейных моделей, как правило, используется стохастический градиентный спуск (Stochastic Gradient Descent, SGD).

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

То есть мы случайно выбираем несколько строк из таблицы и подставляем их в алгоритм, делаем шаг в сторону минимума и повторяем это множество раз, пока алгоритм не сойдётся к приемлемому значению или пока не закончатся итерации (в реализации всегда задаётся максимум итераций на случай, если алгоритм не сойдётся и будет «блуждать по холмам» вечно).

На рисунке ниже приведено сравнение графиков «блуждания» точки в пространстве функции потерь (вид сверху):

Благодаря таким случайным колебаниям у нас появляется больше возможностей «выкарабкаться» из локальных минимумов и дойти до глобального минимума.

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

Чтобы управлять шагами, как раз и существует параметр темпа обучения. Он позволяет управлять размером шага градиентного спуска.

• Слишком маленькие значения \eta уменьшают скорость сходимости, есть шанс не сойтись за отведенное количество итераций
• Слишком большие значения \eta дают эффект бесконечного блуждания - алгоритм коллеблется вокруг точки минимума, но все время "проезжает" ее из-за слишком высокой скорости

Наиболее эффективной оказывается следующая стратегия:

Будем брать большой шаг в начале обучения и уменьшать его постепенно, приближаясь к минимуму, чтобы не «выпрыгнуть» из точки минимума.

В реализации SGD в sklearn данная стратегия используется по умолчанию. Параметр регулируется в процессе обучения — он уменьшается с ростом числа итераций по формуле:

\eta_t = \frac{\eta_0}{t^p}

где

• \eta_0 — начальное значение темпа обучения,
• p — мощность уменьшения темпа (задаётся пользователем).

Ещё один важный момент, на который стоит обратить внимание при работе с градиентным спуском — это обязательное масштабирование факторов (приведение факторов к единому масштабу или к единым статистическим характеристикам), если их несколько.

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

Будем использовать реализацию стохастического градиентного спуска для линейной регрессии из библиотеки sklearn — SGDRegressor. Она находится в том же модуле linear_model.

Численное Решение: Реализация В Sklearn

data = pd.read_csv('https://raw.githubusercontent.com/merion-networks/data-science-course/refs/heads/main/unconv.csv')
data.head()
X = data.drop('Prod', axis=1)
y = data['Prod']

X

Будем использовать реализацию стохастического градиентного спуска для линейной регрессии из библиотеки sklearn — SGDRegressor. Она находится в том же модуле linear_model.

У класса SGDRegressor есть множество параметров. Наиболее важные:

• random_state- отвечает за число, на основе которого происходит генерация случайных чисел. Напомним, в SGD случайность присутствует в инициализации параметров и выборе части из набора данных.
• loss - функция потерь.
  По умолчанию используется squared_loss — уже привычная нам квадратичная ошибка. Но могут использоваться и несколько других. Например, значение "huber" определяет функцию потерь Хьюбера. Эта функция менее чувствительна к наличию выбросов, чем квадратичная ошибка. Полный список возможных функций потерь.
• max_iter - максимальное количество итераций, выделенное на сходимость. Значение по умолчанию — 1000.
• learning_rate - режим управления темпом обучения. Значение по умолчанию — 'invscaling'. Этот режим уменьшает темп обучения по формуле, которую мы рассматривали ранее: \eta_t = \frac{\eta_0}{t^p}.
  Если вы не хотите, чтобы темп обучения менялся на протяжении всего обучения, то можете выставить значение параметра на "constant".
• eta0 - начальное значение темпа обучения . Значение по умолчанию — 0.01. Если параметр learning_rate="constant", то значение этого параметра будет темпом обучения на протяжении всех итераций.
• power_t - значение мощности t уменьшения \eta в формуле \eta_t = \frac{\eta_0}{t^p}. То есть данный параметр отвечает за степень знаменателя (чем больше степень, тем быстрее уменьшается значение темпа обучения с каждой итерацией). Значение по умолчанию — 0.25.

Попробуем обучить модель на наших данных:

# Создаём Объект Класса Линейной Регрессии С SGD
sgd_lr = linear_model.SGDRegressor(
    random_state=42,
    eta0=0.1
)
# Обучаем Модель — Ищем Параметры По Методу SGD
sgd_lr.fit(X, y)

Оценим качество модели по R^2:

# Делаем Прогноз Для Всей Выборки
y_hat = sgd_lr.predict(X)

# Рассчитываем Коэффициент Детерминации
print('R^2-score: {:.2f}'.format(metrics.r2_score(y, y_hat)))

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

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

Зато можно посмотреть полученные в результате оптимизации коэффициенты:

# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X.columns, 'Coefficients': sgd_lr.coef_})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df =pd.DataFrame({'Features': ['Intercept'], 'Coefficients': sgd_lr.intercept_})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

Все коэффициенты имеют запредельные значения (9-11 степени числа 10). Это типичная картина расходящегося градиентного спуска: алгоритм не достиг точки минимума по каким-то причинам. Такие высокие значения коэффициентов означают, что модель является неустойчивой.

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

Воспользуемся классом StandardScaler из модуля preprocessing библиотеки sklearn и стандартизируем наши данные.

from sklearn import preprocessing

# Cоздаем Объект Для Стандартизации
scaler = preprocessing.StandardScaler()
# Выполняем Стандартизацию Данных
X_sc = scaler.fit_transform(X)
X_sc = pd.DataFrame(X_sc, columns=X.columns)
X_sc.describe()

Обучим модель еще раз:

# Создаём Объект Класса Линейной Регрессии С SGD
sgd_lr = linear_model.SGDRegressor(
    random_state=42,
    eta0=0.01
)
# Обучаем Модель — Ищем Параметры По Методу SGD
sgd_lr.fit(X_sc, y)

# Делаем Прогноз Для Всей Выборки
y_hat = sgd_lr.predict(X_sc)

# Рассчитываем Коэффициент Детерминации
print('R^2-score: {:.2f}'.format(metrics.r2_score(y, y_hat)))

Проверим коэффициенты:

# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X.columns, 'Coefficients': sgd_lr.coef_})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df = pd.DataFrame({'Features': ['Intercept'], 'Coefficients': sgd_lr.intercept_})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

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

При этом качество модели по метрике R^2 абсолютно идентично.

Проверим, что SGD находит те же коэффициенты, что аналитическое решение.

# Создаём Объект Класса Линейной Регрессии
lr = linear_model.LinearRegression()
# Обучаем Модель — Ищем Параметры
lr.fit(X_sc, y)

# Делаем Прогноз Для Всей Выборки
y_hat = sgd_lr.predict(X_sc)

# Рассчитываем Коэффициент Детерминации
print('R^2-score: {:.2f}'.format(metrics.r2_score(y, y_hat)))


# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X.columns, 'Coefficients': sgd_lr.coef_})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df = pd.DataFrame({'Features': ['Intercept'], 'Coefficients': sgd_lr.intercept_})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

Результаты идентичные.

Поэкспериментируем с параметрами SGD.

Попробуем взять большое значение темпа обучения.

# Создаём Объект Класса Линейной Регрессии С SGD
sgd_lr = linear_model.SGDRegressor(
    random_state=42,
    eta0=0.5,
    learning_rate='constant', #режим темпа обучения — константа
)
# Обучаем Модель — Ищем Параметры По Методу SGD
sgd_lr.fit(X_sc, y)

sgd_lr.score(X_sc, y)

R^2 < 0 - SGD разошёлся из-за слишком высокого темпа обучения.

А теперь слишком маленькое:

# Создаём Объект Класса Линейной Регрессии С SGD
sgd_lr = linear_model.SGDRegressor(
    random_state=42,
    eta0=0.00001,
    max_iter=100000,
    learning_rate='constant', #режим темпа обучения — константа
)
# Обучаем Модель — Ищем Параметры По Методу SGD
sgd_lr.fit(X_sc, y)

sgd_lr.score(X_sc, y)

Мы видим предупреждение (warning), которое говорит о том, что алгоритму не хватило количества итераций (max_iter), чтобы добраться до минимума. То есть SGD не дошёл до точки минимума из-за слишком низкого темпа обучения.

Преимущества И Недостатки

Преимущества

• Градиентный спуск — простой и мощный алгоритм оптимизации, который позволяет итеративно находить минимум функции потерь и тем самым находить оптимальные параметры модели.

• Благодаря своей простоте алгоритм обладает минимальной вычислительной сложностью и работает быстрее, чем метод наименьших квадратов, даже на огромных наборах данных с тысячами признаков.

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

  Примеры других функций потерь:
  • Log Loss (логистическая функция потерь):

L(y_i, \hat{y_i}) = log(1 + exp(-y_i \cdot \hat{y_i}))

* **Hinge Loss**:

L(y_i, \hat{y_i}) = max(0, 1 - y_i \cdot \hat{y_i})

* **Huber Loss**:

$$L(y, \hat{y}) = \begin{cases} \frac{1}{2}(y - \hat{y})^2, & \text{if } |y - \hat{y}| \leq \delta \ \delta(|y - \hat{y}| - \frac{1}{2}\delta), & \text{if } |y - \hat{y}| > \delta \end{cases}$$ \delta - некоторое число между 0 и 1

• Возможность инкрементального обучения. Есть возможность дообучить модель на новых данных в режиме реального времени. Повторный вызов fit() уточняет уже существующие параметры модели.

Недостатки

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

Полиномиальная Регрессия

Общее Представление

Полиномиальная регрессия (Polynomial Regression) — это более сложная модель, чем линейная регрессия. Вместо уравнения прямой используется уравнение полинома (многочлена). Степень полинома (максимальная степень) может быть сколь угодно большой: чем больше степень, тем сложнее модель.

Пример полиномиальных моделей

• Полином степени 2 с 1-им признаком x:

\hat{y} = \beta_0 + \beta_1 x_ + \beta_2 x^2

• Полином степени 2 с 2-мя признаками x_1 и x_2:

\hat{y} = \beta_0 + \beta_1 x_{1} + \beta_2 x_{1}^2 + \beta_3 x_{2} + \beta_4 x_{2}^2 + \beta_5 x_{1} x_{2}

• Полином степени 3 с 2-мя признаками x_1 и x_2:

\hat{y} = \beta_0 + \beta_1 x_{1} + \beta_2 x_{1}^2 + \beta_3 x_{2} + \beta_4 x_{2}^2 + \beta_5 x_{1} x_{2}+ \beta_6 x_{1}^2 x_{2} + \beta_7 x_{1} x_{2}^2 + \beta_8 x_{2}^3

Геометрически полиномиальная регрессия - это некоторая аппроксимирующая поверхность.

Примеры:

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

Однако, сложность модели является ее преимуществом и недостатком.

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

Это отсылает нас к необходимости валидации модели при использовании полиномиальной регрессии.

Полиномиальные Признаки

Заметим, что степени при $x$-ах можно тоже считать своего рода искусственными признаками в данных. Они называются полиномиальными признаками.

Например в уравнении модели:

\hat{y} = \beta_0 + \beta_1 x_{1} + \beta_2 x_{1}^2 + \beta_3 x_{2} + \beta_4 x_{2}^2 + \beta_5 x_{1} x_{2}

можно обозначить z_{1} = x_{1}, z_{i2} = x_{i1}^2, z_{3} = x_{2}, z_{4} = x_{2}^2, z_{5} = x_{1} x_{2}

Тогда получим обычное уравнение линейной регрессии f(z):

\hat{y} = \beta_0 + \beta_1 z_{1} + \beta_2 z_{2} + \beta_3 z_{3} + \beta_4 z_{4} + \beta_5 z_{5}

Поэтому полиномиальная регрессия — это та же линейная регрессия, построенная на исходном датасете с дополнительными признаками.

Создать полиномиальные признаки в sklearn можно с помощью с помощью объекта класса PolynomialFeatures из модуля preprocessing.

У класса есть два важных параметра:

• degree — степень полинома. По умолчанию используется степень 2.
• include_bias — включать ли в результирующую таблицу столбец из единиц (x в степени 0). По умолчанию стоит True, но лучше выставить его в значение False, так как столбец из единиц и так добавляется в методе наименьших квадратов.

Для того чтобы подогнать генератор и рассчитать количество комбинаций степеней, мы используем метод fit(), а чтобы сгенерировать новую таблицу признаков, в которую будут включены полиномиальные признаки, используется метод transform(), в который нужно передать выборки:

#Создаём генератор полиномиальных признаков 2 степени
poly = preprocessing.PolynomialFeatures(degree=2, include_bias=False)
poly.fit(X)

X_poly = poly.transform(X)

X_poly.shape

Чтобы узнать, какие столбцы каким комбинациям признаков соответствуют, можно использовать метод get_feature_names_out()

poly.get_feature_names_out()

Для удобства можно составить новый DataFrame из полученных признаков:

X_poly = pd.DataFrame(X_poly, columns=poly.get_feature_names_out())

X_poly.head()

Полиномиальная Регрессия: Реализация В Sklearn

data = pd.read_csv('https://raw.githubusercontent.com/merion-networks/data-science-course/refs/heads/main/unconv.csv')
data.head()
X = data.drop('Prod', axis=1)
y = data['Prod']

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

Предварительно давайте позаботимся о валидации и выделим 20% наших на тестовую выборку:

from sklearn.model_selection import train_test_split

#Разделяем выборку на тренировочную и тестовую в соотношении 80/20
#Устанавливаем random_state для воспроизводимости результатов
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=40)
#Выводим результирующие размеры таблиц
print('Train:', X_train.shape, y_train.shape)
print('Test:', X_test.shape, y_test.shape)

#Создаём генератор полиномиальных признаков
poly = preprocessing.PolynomialFeatures(degree=3, include_bias=False)
poly.fit(X_train)

#Генерируем полиномиальные признаки для тренировочной выборки
X_train_poly = poly.transform(X_train)
#Генерируем полиномиальные признаки для тестовой выборки
X_test_poly = poly.transform(X_test)

# Преобразуем Результаты В DataFrame Для Удобства

X_train_poly = pd.DataFrame(X_train_poly, columns=poly.get_feature_names_out())
X_test_poly = pd.DataFrame(X_test_poly, columns=poly.get_feature_names_out())

#Выводим результирующие размерности таблиц
print(X_train_poly.shape)
print(X_test_poly.shape)

Обучим две модели:

• линейную регрессию
• полиномиальную регрессию

Сравним качество полученных моделей по метрике MAPE

Напоминание:

MAPE = \frac{100 \%}{n} \sum_{i=1}^n \frac{|y_i - \hat{y_i}|}{|y_i|}

Простая линейная регрессия:

#Создаём объект класса LinearRegression
lr = linear_model.LinearRegression()
#Обучаем модель по МНК
lr.fit(X_train, y_train)

#Делаем предсказание для тренировочной выборки
y_train_predict = lr.predict(X_train)
#Делаем предсказание для тестовой выборки
y_test_predict = lr.predict(X_test)

print("Train MAPE: {:.3f}".format(metrics.mean_absolute_percentage_error(y_train, y_train_predict)  * 100))
print("Test MAPE: {:.3f}".format(metrics.mean_absolute_percentage_error(y_test, y_test_predict) * 100) )

Полиномиальная регрессия:

#Создаём объект класса LinearRegression
lr_poly = linear_model.LinearRegression()
#Обучаем модель по МНК
lr_poly.fit(X_train_poly, y_train)

#Делаем предсказание для тренировочной выборки
y_train_predict = lr_poly.predict(X_train_poly)
#Делаем предсказание для тестовой выборки
y_test_predict = lr_poly.predict(X_test_poly)

print("Train MAPE: {:.3f}".format(metrics.mean_absolute_percentage_error(y_train, y_train_predict)  * 100))
print("Test MAPE: {:.3f}".format(metrics.mean_absolute_percentage_error(y_test, y_test_predict)  * 100))

lr_poly.coef_

Видим, что качество полиномиальной модели превосходит качество линейной модели.

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

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

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

Регуляризация

Регуляризация — это способ уменьшения переобучения моделей машинного обучения путём намеренного увеличения смещения модели для уменьшения её разброса.

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

• предотвратить переобучение модели;
• включить в функцию потерь штраф за переобучение;
• обеспечить существование обратной матрицы (X^TX)^{-1};
• не допустить огромных коэффициентов модели.

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

L(\beta) \rightarrow 0

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

Идея регуляризации состоит в наложении ограничения на вектор весов (часто говорят — наложение штрафа за высокие веса).

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

\tilde{L}(\beta) = L(\beta) + \alpha ||\beta||_p \tilde{L}(\beta)\rightarrow \min_{\beta}

В случае квадратичной функции потерь получаем:

\frac{1}{n}||y - X \cdot {\beta}||^2 + \alpha ||\beta||_p \rightarrow \min_{\beta}

где ||\beta||_p - норма весов степени p вычисляемая как:

||\beta||_p = \sum_{i=0}^m|\beta_i|^p

На первый взгляд, ничего особо не изменилось по сравнению с изначальной задачей оптимизации. Добавилось только одно слагаемое — \alpha ||\beta||_p. Однако это слагаемое имеет важное значение:

• Во-первых, оно не дает найти истинный минимум функции потерь, так как минимум итоговой функции несколько смещен относительно истинного минимума
• Во-вторых, оно не допускает высоких значений параметров, так как, если ||\beta||_p будет слишком большой, минимума найдено не будет.

Таким образом, именно это маленькое слагаемое позволяет нам побороть переобучение.

Коэффициент \alpha принято называть коэффициентом регуляризации. Он отвечает за «силу» регуляризации. Чем он больше, тем меньшие значения может принимать слагаемое \alpha ||\beta||_p, то есть тем сильнее ограничения на норму весов (штраф).

Порядок нормы p в общем случае может быть любой — главное, чтобы она была больше 1. Однако на практике распространены только первая и вторая степени, так называемые $L_1$- и $L_2$-регуляризации.

Обратите внимание на сумму под знаком корня. У нас она начинается с 0. Однако иногда в литературе можно встретить i=1 (например здесь), то есть свободный член \beta_0 не регуляризируют. Мы будем придерживаться реализации в sklearn, в которой \beta_0 всё-таки включается в регуляризацию.

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

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

$L_1$-регуляризация (Lasso)

$L_1$-регуляризацией, Lasso (Least Absolute Shrinkage and Selection Operator), называется регуляризация, в которой порядок нормы p=1.

\tilde{L}(\beta) = L(\beta) + \alpha ||\beta||_1 \rightarrow \min_{\beta} или \tilde{L}(\beta) = \frac{1}{n} ||y - X \cdot {\beta}||^2 + \alpha \sum_{i=0}^m|\beta_i|\rightarrow \min_{\beta}

Таким образом, в случае -регуляризации мы ограничиваем сумму модулей весов модели. Напомним, такая величина называется нормой Манхэттена (расстоянием городских кварталов).

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

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

В sklearn $L_1$-регуляризация реализована в классе Lasso, а заданная выше оптимизационная задача решается алгоритмом координатного спуска (Coordinate Descent).

from sklearn.linear_model import Lasso

$L_2$-регуляризация (Ridge)

$L_2$-регуляризация (Ridge), или регуляризация по Тихонову — это регуляризация, в которой порядок нормы p=2.

Тогда, если подставить p=2 в наши формулы, то оптимизационная задача в случае будет иметь вид:

\tilde{L}(\beta) = L(\beta) + \alpha ||\beta||_2 \rightarrow \min_{\beta} или \tilde{L}(\beta) = \frac{1}{n} ||y - X \cdot {\beta}||^2 + \alpha \sum_{i=0}^m(\beta_i)^2\rightarrow \min_{\beta}

Видно, что норма порядка p=2 в случае $L_2$-регуляризации мы накладываем ограничение на длину вектора весов \beta.

У данной задачи даже есть аналитическое решение, полученное математиком Тихоновым, вот оно:

\beta=(X^T \cdot X + \alpha \cdot E)^{-1} \cdot X^T \cdot y

здесь E - единичная матрица размера (m+1) \times (m+1)

Преимущество этой формулы в том, что, если \alpha > 0, то матрица X^TX + \alpha E гарантированно является невырожденной, а значит у нее 100% есть обратная. Так получается за счёт того, что по диагонали матрицы мы добавляем поправки, которые создают линейную независимость между столбцами матрицы.

За реализацию линейной регрессии с $L_2$-регуляризацией в sklearn отвечает класс Ridge. Основной параметр модели, на который стоит обратить внимание — alpha, коэффициент регуляризации из формулы Тихонова.

from sklearn.linear_model import Ridge

Регуляризация: Сравнение Моделей

А теперь давайте применим регуляризацию на нашем датасете и сравним линейные модели.

Так как у всех моделей единый интерфейс (методы fit(), predict() и атрибут coef_), то для избежания дублирования кода можно написать небольшую функцию, которая будет выполнять обучение переданной в нее модели на обучающей выборке, оценивать качество модели на обучающей и тестовой выборках по переденной метрике, а также считать сумму модулей коэффициентов модели.

def validate_model(model, degree, X_train, y_train, X_test, y_test, metric):
    #Инициализируем стандартизатор StandardScaler
    scaler = preprocessing.StandardScaler()
    #Подгоняем параметры стандартизатора (вычисляем среднее и СКО)
    scaler.fit(X_train)
    #Производим стандартизацию тренировочной выборки
    X_train = scaler.transform(X_train)
    #Производим стандартизацию тестовой выборки
    X_test = scaler.transform(X_test)

    #Создаём генератор полиномиальных признаков
    poly = preprocessing.PolynomialFeatures(degree=degree, include_bias=False)
    poly.fit(X_train)
    #Генерируем полиномиальные признаки для тренировочной выборки
    X_train = poly.transform(X_train)
    #Генерируем полиномиальные признаки для тестовой выборки
    X_test = poly.transform(X_test)

# Обучаем Модель
    model.fit(X_train, y_train)
# Делаем Предсказание Для Тренировочной И Тестовой Выборок
    y_train_pred = model.predict(X_train)
    y_test_pred = model.predict(X_test)

# Вычисляем Метрики На Тренировочной И Тестовой Выборках
    train_score = metric(y_train, y_train_pred)
    test_score = metric(y_test, y_test_pred)

# Считаем Сумму Коэффициентов
    coef_sum = abs(model.coef_).sum()
    return model, degree, train_score, test_score, coef_sum

lr = linear_model.LinearRegression()
lasso_lr = linear_model.Lasso(alpha=0.1)
ridge_lr = linear_model.Ridge(alpha=10)

# Метрика Для Оценки Качества
metric = metrics.mean_absolute_percentage_error

scores_df = pd.DataFrame(
    data = [
# Линейная Регрессия
        validate_model(
            lr, 1,
            X_train, y_train,
            X_test, y_test,
            metric
        ),
# Линейная Регрессия С L1-регуляризацией
        validate_model(
            lasso_lr, 1,
            X_train, y_train,
            X_test, y_test,
            metric
        ),
# Линейная Регрессия С L2-регуляризацией
        validate_model(
            ridge_lr, 1,
            X_train, y_train,
            X_test, y_test,
            metric
        ),
# Полиномиальная Регрессия
        validate_model(
            lr, 3,
            X_train, y_train,
            X_test, y_test,
            metric
        ),
# Полиномиальная Регрессия С L1-регуляризацией
        validate_model(
            lasso_lr, 3,
            X_train, y_train,
            X_test, y_test,
            metric
        ),
# Полиномиальная Регрессия С L2-регуляризацией
        validate_model(
            ridge_lr, 3,
            X_train, y_train,
            X_test, y_test,
            metric
        )
    ],
    columns=['model', 'degree', 'train_score', 'test_score', 'coef_sum']
)

scores_df.sort_values(by=['test_score', 'train_score'])

Видим, что лучшей моделью оказалась модель Lasso со степенью полинома 3.

Давайте попробуем подобрать для нее коэффициент регуляризации \alpha.

alpha_list = np.arange(0.01, 10, 0.1)
# Создаём Пустые Списки, В Которые Будем Добавлять Результаты
train_scores = []
test_scores = []

metric = metrics.mean_absolute_percentage_error
for alpha in alpha_list:
# Создаём Объект Класса Линейной Регрессии С L1-регуляризацией
    lasso_lr = linear_model.Lasso(alpha=alpha, max_iter=10000)
# Обучаем Модель И Считаем Метрики
    model, degree, train_score, test_score, coef_sum = validate_model(
        lasso_lr, 3,
        X_train, y_train,
        X_test, y_test,
        metric
    )
    train_scores.append(train_score)
    test_scores.append(test_score)

scores_df = pd.DataFrame(
    {
        'alpha': alpha_list,
        'train_score': train_scores,
        'test_score': test_scores
    }
)

fig = px.line(
    data_frame=scores_df,
    x='alpha',
    y=['train_score', 'test_score']
)
fig.show()

lasso_lr = linear_model.Lasso(alpha=0.6, max_iter=10000)

Эластичная Регуляризация (Elastic-net) *

Последний вид регуляризации (хотя их на самом деле больше), который мы рассмотрим, называется Elastic-Net (эластичная сетка). Это комбинация L_1 - и L_2 -регуляризации.

\tilde{L}(\beta) = L(\beta) + \alpha \cdot {\lambda} \cdot ||\beta||_1 + \alpha \cdot \frac{(1-\lambda)}{2} \cdot ||\beta||_2 \rightarrow \min_{\beta} или \tilde{L}(\beta) = \frac{1}{n}||y - X \cdot {\beta}||^2 + \alpha \cdot {\lambda} \sum_{i=0}^m|\beta_i| + \alpha \cdot \frac{(1-\lambda)}{2} \sum_{i=0}^m(\beta_i)^2\rightarrow \min_{\beta}

Здесь коэффициенты и отвечают за вклад слагаемых регуляризации.

Если \alpha=0, получаем классическую МНК-задачу оптимизации. Если \alpha \neq 0, \lambda=1, получаем Lasso-регрессию. Если \alpha \neq 0, \lambda=0, получаем Ridge-регрессию с коэффициентом \frac{\alpha}{2}.

Аналитического решения у этой задачи нет, поэтому для её решения в sklearn, как и для модели Lasso, используется координатный спуск.

В sklearn эластичная сетка реализована в классе ElasticNet из пакета с линейными моделями — linear_model.

За коэффициент \alpha отвечает параметр alpha, за коэффициент \lambda — l1_ratio.

from sklearn.linear_model import ElasticNet

Преимущества И Недостатки Линейной Регрессии

Преимущества:

1. Простота интерпретации:
   Линейная регрессия предоставляет простую формулу, которая позволяет легко интерпретировать взаимосвязь между независимыми и зависимыми переменными.
2. Вычислительная эффективность:
   Линейная регрессия требует значительно меньше вычислительных ресурсов, чем более сложные модели, что делает ее хорошим выбором в качестве бейзлайновой модели.
3. Устойчивость к переобучению:
   Из-за своей простоты линейная регрессия является алгоритмом, устойчивым к переобучению, а благодаря регуляризации можно забыть о переобучении вовсе. Подбор гиперпараметров для такой линейной модели не так важен, как для более сложных моделей.
4. Масштабируемость прогнозов:
   Модель может неплохо предсказывать на данных, значения которых выходят из границ обучающего набора данных.

Недостатки модели:

1. Ограниченность допущений:
   Линейная регрессия основана на ряде предположений, включая линейность, нормальность распределения остатков, отсутствие мультиколлинеарности и другие. Как правило, в реальных задачах эти условия не выполняются и линейные модели показывают низкое качество, нежели на "игрушечных" задачах.
2. Неустойчивость к выбросам:
   Линейная регрессия как и все аналитические модели может постраиваться под выбросы в данных, особенно при отсутствии регуляризации.
3. Работа с категориальными признаками:
   Аналитические модели, к коим относятся все линейные модели предназначены для работы с числовыми характеристиками, а не с категориальными. Модель плохо понимает закодированные категориальные признаки и слабо учитывает их при построении уравнения.
4. Неустойчивость к высоким размерностям:
   Как правило, линейные модели плохо показывают себя в случае большого количества признаков в данных. Это связано с ухудшением качества сходимости методов оптимизации в пространствах больших размеростей (проклятие размерности)

Резюме

• Линейная регрессия - простой и эффективный алгоритм для решения задач регрессии, предполагающий, что зависимость между таргетом и фичами линейная
• Как правило для поиска коэффициентов линейной регрессии (параметров модели) используется метод наименьших квадратов (OLS), предполагающий минимизацию квадратичной функции потерь.
• Искать минимум функции потерь можно несколькими способами:
  • Аналитический метод, то есть "в лоб", используя матричную формулу OLS для оценки коэффициентов
  • Численными методами, такими как градиентный спуск (GD) и стохастический градиентный спуск (SGD)
• Полиномиальная регрессия - расширение линейной регрессии, которое позволяет моделировать нелинейные зависимости.
  Полиномиальная регрессия является мощной моделью, позволяющей решать сложные задачи, однако у данной модели есть недостаток в виде необходимости оптимального подбора степени полинома, так как высокие степени хоть и дают неплохие результаты, но и повышают риск переобучения модели.
• Регуляризация - метод борьбы с переобучением в линейных моделях. Он основан на идее внесения намеренного смещения в модель с целью уменьшить ее разброс.
  Регуляризация выполняется путем замены оригинальной функции потерь на функцию потерь с добавочным слагающим, зависящим от весов.
• Существует три основных типа регуляризации: L_1 (Lasso), L_2 (Ridge) и ElasticNet.
• $L_1$-регуляризация приводит к занулению коэффициентов при неинформативных признаках.
• $L_2$-регуляризация уменьшает все коэффициенты, но не обнуляет их.
• ElasticNet - комбинация L_1 и L_2 регуляризации.
• Выбор типа модели и коэффициента регуляризации - важный этап построения линейной модели.

Небольшое дополнение про категориальные признаки

Как мы видим, линейные модели рассчитаны на работу с числами.

А как быть, если в датасете присутствуют категориальные фичи, не являющихся числами?

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

Пример:

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

Самый простой способ – использовать one-hot кодирование (one-hot encoding). Пусть исходный признак имеет K уникальных значений. Давайте заменим категориальный столбец на K новых столбцов, которые принимают значения 1 и 0. Единица ставится, если $i$-ый объект принадлежит $k$-ой категории, в противном случае на месте ставится 0.

Конечно, one-hot кодирование – это самый наивный способ работы с категориальными признаками, и для более сложных фичей или фичей с большим количеством значений оно плохо подходит.

Дополнительные Материалы

• Лекции по машинному обучению от В.Л. Аббакумова на ютубчике
• Лекции по машинному обучению от Дмитрия Макарова (есть теория и базовая практика)
• Лекции по машинному обучению от Cornel Tech (английский, есть базовый код, МНОГО формул с выводом и доказательством)
• Учебник по машинному обучению от Яндекса (0 кода, но много математики и красивых картинок)

---

Линейная Разделимость Классов

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

Логистическая Регрессия

Постановка Задачи

Пусть дана выборка объектов и их классов:

\{(x_i, y_i)\}_{i=1}^n = \{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}

где

• x_i = (x_{i1}, x_{i2}, ..., x_{im}), x_i ∈ X, X \subset \mathbb{R} ^{m} - множество объектов
• y_i ∈ Y, Y =\{0, 1, ..., (k-1)\} - множество целевой переменной
• m - количество признаков
• n - размер выборки
• k - количество классов

Необходимо построить модель, восстанавливающую зависимость y от x: y=f(x)

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

Классы обозначим как 0 и 1, т.е Y=\{0, 1\}.

Будем рассматривать теорию на конкретном примере.

В качестве датасета возьмем датасет о раке молочной железы.

Датасет содержит 569 образцов тканей груди, каждый из которых описывается 30 признаками, такими как радиус опухоли, текстура, периметр, площадь, гладкость, компактность и др. Целевая переменная в этом наборе данных - это метка, указывающая на то, является ли опухоль доброкачественной (0) или злокачественной (1).

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

Полное описание датасета

Так как это широкоиспользуемый для экспериментов и обучения датасет, то его можно импортировать прямо из библиотеки sklearn из модуля datasets.

import pandas as pd
import numpy as np
import plotly.express as px

from sklearn import datasets

Импортируем данные:

cancer_data = datasets.load_breast_cancer(as_frame=True)

X = cancer_data.data
y = cancer_data.target

print(X.shape)
print(y.shape)

X.head()

y.head()

В нашем случае объектами x_i являются опухоли, характеризующиеся 30 признаками, т.е m=30, размер выборки n=569

Целевая переменная y_i представляет собой бинарную категориальную переменную Y=\{0, 1\}, где

• y=0 - доброкачественная опухоль
• y=1 - злокачественная опухоль

Необходимо обучить классификатор, который будет определять наличие рака по параметрам опухоли.

X.corrwith(y).sort_values()

Общая Идея

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

\hat{y} = f(x) = x \cdot \beta = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_m x_{m}

или в матричном виде:

\hat{y} = X \cdot \beta

где

$$X=\begin{pmatrix} 1 & x_{11} & x_{12} & ... & x_{1m} \ 1 & x_{21} & x_{22} & ... & x_{2m} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_{n1} & x_{n2} & ... & x_{nm} \end{pmatrix}$$

\hat{y}=(\hat{y_1}, \hat{y_2}, ..., \hat{y_n})^T \beta = (\beta_0, \beta_1, \beta_2, ..., \beta_m)^T

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

Если целевая переменная перекодирована в числа, например Y = \{0, 1\}, то ничто нам не мешает построить регрессию и восстановить зависимость между x и y.

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

L(y, \hat{y}) = \frac{1}{n} \sum_{i=1}^n (y_i-\hat{y_i})^2 = \frac{1}{n} \sum_{i=1}^n (y_i-x_i \cdot \beta)^2

# Строим Линейную Регрессию По 1-ому Признаку
px.scatter(
    data_frame = -X,
    x='worst perimeter',
    y=y,
    trendline="ols", # линия тренда определяется с помощью OLS
    trendline_scope="overall",
    height=400,
    width=800,
    title="Линейная регрессия для классификации"
)

И тут мы сталкиваемся с двумя принципиально важными вопросами:

• Кто сказал, что линейная регрессия гарантировано будет выдавать числа 0 и 1? Вполне возможна ситуация, когда результат линейной регрессии равен -10. К какому классу принадлежит в таком случае объект? Наверное, к классу 0, но кто нам это гарантировал?
• Что делать в том случае, если линейная регрессия выдала число между 0 и 1, например 0.57. К какому классу принадлежит такой объект? Наверное, 1? А почему?
• Все становится еще хуже, если мы рассматриваем не бинарную, а многоклассовую классификацию

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

Однако, решать задачу как-то надо.

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

Обозначим:

• p(x, y=0) - вероятность того, что объект x принадлежит классу 0
• p(x, y=1) - вероятность того, что объект x принадлежит классу 1

В нашем примере: p(x, y=0) - вероятность того, что опухоль x доброкачественная, p(x, y=1) - что она злокачественная.

Вероятности принадлежности к классам 0 и 1 должны обладать двумя важными свойствами:

• Это должны быть числа от 0 до 1: p \in [0, 1]
• В сумме они должны давать 1: p(x, y=0) + p(x, y=1)= 1

Для удобства математических выкладок, в случае бинарной классификации нам достаточно рассматривать только одну вероятность p(x, y=1), так как вероятность принадлежности к классу 0 можно посчитать как 1-p(x, y=1).

То есть мы должны предсказывать вероятность наличия эффекта (наличие рака в нашем случае).

В итоге мы добьёмся того, что будем предсказывать не дискретный категориальный, а непрерывный числовой признак, который лежит в диапазоне [0, 1]. А это уже знакомая нам задача регрессии.

Тут мы и приходим к модели логистической регрессии - регрессии вероятностей принадлежности к классу.

Логистическая Функция

Ключевой вопрос, который перед нами стоит - это как нам заставить модель линейной регрессии выдавать вероятности принадлежности к классам. По сути мы должны свести ответы модели из диапазона (-∞, +∞) к диапазону [0, 1].

Тут нам приходит на помощью логистическая функция (logistic function) \sigma(z), которая вычисляется как:

\sigma(z) = \frac{e^z}{1+e^z} = \frac{1}{1 + e^{-z}}

Здесь e — экспонента или число Эйлера. Это число является бесконечным, а его значение обычно принимают равным 2.718.

График зависимости этой функции от аргумента:

Данная функция чаще всего называется сигмоида.

Важные свойства сигмоиды:

• Значения сигмоиды \sigma(z) лежат в диапазоне от 0 до 1.
• Функция \sigma(z) \rightarrow 0 при z \rightarrow -∞ и \sigma(z) \rightarrow 1 при z \rightarrow +∞
• \sigma(z) < 0.5 при z < 0, \sigma(z) > 0.5 при z>0 и \sigma(z) = 0.5 при z=0.
• Сигмоида имеет очень простую и легко вычисляемую производную:

\frac{d\sigma(z)}{dz} = \sigma(z)(1-\sigma(z))

То есть, если значение сигмоиды в определенной точке уже вычислено, то мы можем использовать это же значение для расчёта значения производной (градиента) в этой же точке. Это сильно упрощает расчёты.

А теперь гвоздь программы:

Давайте в качестве аргумента z в \sigma(z) будем использовать выход модели линейной регрессии.

Получим:

f(x) = \sigma(x \cdot \beta) = \frac{1}{1 + e^{-x \cdot \beta}}

Полученное выражение как раз и называется моделью логистической регрессии (logistic regression).

Если мы поберем "хорошие" веса \beta, то результат логистической регрессии можно будет интерпретировать как вероятность принадлежности объекта x к классу 1: p(y=1) = f(x).

Вероятность принадлежности объекта к классу 0 тогда будет вычисляться как p(y=0) = 1 - p(y=1) = 1 - \sigma(x \cdot \beta).

Класс объекта определяется как тот, что имеет максимальную вероятностью:

\hat{y} = \arg \max_{y \in Y} (p)

В случае бинарной классификации это будет сответствовать тому, что можно если \sigma(x \cdot \beta) ≥ 0.5, то мы относим объект к классу 1, в противном случае - к классу 0: $$\hat{y_i} = \begin{cases} 1 , & \text{if } \sigma(x \cdot \beta) ≥ 0.5\ 0 , & \text{if } \sigma(x \cdot \beta) < 0.5\ \end{cases}$$

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

Пример

Представим, что мы построили модель логистической регрессии на задаче классификации опухолей на 2 признаках ('mean radius' и 'worst radius') и нашли оптимальные коэффициенты.

В результате обучения логистической регрессии мы получили следующие коэффициенты:

\beta= (\beta_0, \beta_1, \beta_2)^T = (16.59, -2.2, 1.37)^T

То есть уравнение линейной модели имеет вид:

\beta ⋅ x = -16.59 -2.2 \cdot x_{1} + 1.37 \cdot x_{2}

К нам приходит новая опухоль со следующими характеристиками: x^* = (17.5, 15). Хотим определить, является ли данная опухоль злокачественной опираясь только на 2 этих признака.

Смотрим на реализацию:

def linear(X, beta):
# Добавляем Столбец Из Единиц
    ones = np.ones((X.shape[0], 1))
    X = np.hstack((ones, X))
# Вычисляем Прогноз Линейной Модели
    z = X @ beta
    return z

def sigmoid(z):
# Вычисляем Сигмоиду
    sigma = 1 / (1 + np.exp(-z))
    return sigma

def predict_proba(X, beta):
# Делаем Предсказание Линейной Моделью
    z = linear(X, beta)
# Вычисляем Сигмоиду
    sigma = sigmoid(z)
# Вычисляем Вероятности
    p = np.array([1 - sigma, sigma])
    return p.reshape(-1)

def predict(X, beta):
    p = predict_proba(X, beta)
# Возвращаем Индекс Максимального Элемента
    y_hat = np.argmax(p)
    return y_hat

beta = np.array([[16.59], [-2.2], [1.37]])
x = np.array([[17.5, 15]])

p = predict_proba(x, beta)
y_hat = predict(x, beta)
print('Вероятности принадлежности к классам: {}'.format(p))
print('Предсказанный класс: {}'.format(y_hat))

А теперь возьмем другую опухоль: x = (13.5, 18):

x = np.array([[13.5, 18]])

p = predict_proba(x, beta)
y_hat = predict(x, beta)
print('Вероятности принадлежности к классам: {}'.format(p))
print('Предсказанный класс: {}'.format(y_hat))

Геометрическая Интерпретация

Разберёмся с геометрией.

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

z = \beta ⋅ x = \beta_0 + \beta_1 \cdot x_{1} + \beta_2 \cdot x_{2}

Для конкретики пусть

\beta= (\beta_0, \beta_1, \beta_2)^T = (16.59, -2.2, 1.37)^T

Тогда:

z = -16.59 -2.2 \cdot x_{1} + 1.37 \cdot x_{2}

Если рассматривать уравнение линейной регрессии отдельно от сигмоиды, то это уравнение задает плоскость с коэффициентами \beta_0, \beta_1 и \beta_2 . Эта плоскость проходит таким образом, чтобы отделить объекты разных классов друг от друга.

Если визуализировать данные в 2D, а классы обозначить цветом, то уравнение будет задавать разделяющую прямую.

Причем, смещение этой прямой относительно оси x_2 будет равно \frac{\beta_0}{\beta_2}.

А коэффициент наклона: -\frac{\beta_1}{\beta_2}

Математически подстановка в уравнение плоскости точки, которая не принадлежит ей (находится ниже или выше), означает вычисление расстояния от этой точки до плоскости.

• Если точка находится ниже плоскости, расстояние будет отрицательным (z < 0).
• Если точка находится выше плоскости, расстояние будет положительным (z > 0).
• Если точка находится на самой плоскости, z = 0.

Подстановка отрицательных чисел в сигмоиду приведёт к вероятности p(y=1) > 0.5, а постановка положительных — к вероятности p(y=1) < 0.5.

Чем больше расстояние от точки, находящейся выше разделяющей плоскости, до самой плоскости, тем больше оценка p(y=1). И наоборот - чем больше расстояние от точки, находящейся ниже плоскости, тем ниже p(y=1) (как следствие выше p(y=0))

Можно построить тепловую карту, которая показывает, чему равны вероятности в каждой точке пространства:

Для случая зависимости целевого признака от трёх факторов мы получаем уравнение:

z = \beta ⋅ x = \beta_0 + \beta_1 \cdot x_{1} + \beta_2 \cdot x_{2} + \beta_3 \cdot x_{3}

Уравнение задаёт плоскость в четырёхмерном пространстве. Но если вспомнить, что y — категориальный признак и классы можно обозначить цветом, то получится перейти в трёхмерное пространство:

В общем случае линейная модель задает разделяющую плоскость в $m$-мерном пространстве (гиперплоскость).

Метод Максимального Правдоподобия (MLE). Log Loss

Для поиска параметров нам нужна правильная функция потерь L(\beta), оптимизировав которую мы с вами найдем оптимальные параметры модели \beta.

Ранее мы с вами сказали о том, что квадратичная функция (MSE), используемая в МНК (OLS) нам не подходит по ряду причин.

А что тогда подходит?

Если для регрессионных моделей основной метод построения - МНК, то для моделей, имеющих вероятностную интерпретацию в большинстве случаев используется метод максимального правдоподобия (Maximum Likelihood Estimation — MLE).

Правдоподобие — это оценка того, насколько вероятно получить истинное значение целевой переменной y при данных x и параметрах \beta.

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

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

L(\beta) = \frac{1}{n} \sum_{i=1}^n (y_i \cdot log(p_{i}) + (1-y_i) \cdot log(1 - p_{i}))

где

• y_i - истинный класс $i$-ого объекта (0 или 1)
• p_i = p_i(y=1) = \sigma(\beta \cdot x_i) - вероятность принадлежности $i$-ого объекта к классу 1, вычисляемые через сигмоиду
• log - логарифм, как правило используется натуральный логарифм по основанию e

Чем выше значение функции L(\beta), тем больше вероятности p_i сооветствуют истинным классам в данных, то есть, тем лучше модель, а значит тем более подходящими являются параметры \beta.

Мы хотим найти такие параметры \beta, при которых правдоподобие было бы максимальным:

L(\beta) \rightarrow \max_{\beta}

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

L(\beta) = - \frac{1}{n} \sum_{i=1}^n (y_i \cdot log(p_{i}) + (1-y_i) \cdot log(1 - p_{i}))

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

Разберемся как работает выражение:

- (y_i \cdot log(p_{i}) + (1-y_i) \cdot log(1 - p_{i}))

• Сценарий №1:
  Представим, что нам на вход пришел объект x_i, причем y_i = 1, т.е злокачественная опухоль. Тогда мы получаем, что второе слагаемое обнуляется:

- (1 \cdot log(p_{i}) + (1-1) \cdot log(1 - p_{i})) = - log(p_{i})

То есть, чем выше вероятность p_{i}=\sigma(x_i \cdot \beta) мы предскажем, тем меньше будет ошибка. То есть для объекта класса 1 сигмоида должна выдавать высокие значения.

• Сценарий №2:
  Представим, что к нам пришел объект x_i с y_i = 0, т.е. доброкачественная опухоль. Тогда обнуляется первое слагаемое:

- (0 \cdot log(p_{i}) + (1-0) \cdot log(1 - p_{i})) = - log(1 - p_{i})

В этом случае, чем больше будет вероятность p_{i}=\sigma(x_i \cdot \beta), тем больше будет ошибка. То есть для объектов класса 0 сигмоида должна выдавать низкие значения p_i

Таким образом, наша цель - подобрать такие параметры \beta при которых среднее значение log loss по всей выборке было бы минимальным:

L(\beta) \rightarrow \min_{\beta}

Маленькое дополнение №1.

Так как, гипотетически, при накоплении компьютерных вычислений может получиться так, что \sigma(z) = 0 или \sigma(z) = 1, а по правилам математики log(0) не существует, то в программной реализации вычисления log loss вероятности предварительно приводят к диапазону [\epsilon, 1 - \epsilon], где \epsilon - очень маленькое число, близкое к 0.

Маленькое дополнение №2.

В некоторых обучающих материалах (например, здесь и здесь) по бинарной классификации вы можете найти обозначение классов 0 и 1 как -1 и +1. В таком случае логистическая функция потерь будет иметь очень простую формулу:

L(\beta) = - \frac{1}{n} \sum_{i=1}^n log(1 + e^{-y_i \cdot x_i \cdot \beta})

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

Пример

Рассмотрим пример расчёта функции log loss.

Пусть у нас есть выборка из четырёх объектов (опухолей) с двумя признаками: x_{1} и x_{2}.

В результате поиска коэффициентов мы получили, что нам подходят следующие коэффициенты:

\beta= (\beta_0, \beta_1, \beta_2)^T = (16.59, -2.2, 1.37)^T

Тогда уравнение плоскости:

z = -16.59 -2.2 \cdot x_{2} + 1.37 \cdot x_{1}

По формуле сигмоиды мы рассчитали вероятности принадлежности к классам для всех 4-ех объектов и получили:

| i | $\large y$| \large p |$\large (1-p)$| | | | | 1 | 0 | 0.2 | 0.8 | | 2 |0 | 0.8 |0.2 | | 3 |1 | 1 |0 | | 4 |1 | 0.6 |0.4 |

Рассчитаем logloss:

def log_loss(y, p):
# Вычисляем Размер Выборки
    n = y.shape[0]
# Зажимаем Значения Вероятностей, Чтобы Не Получить Логарифм От 0
    p = np.clip(p, 1e-10, 1 - 1e-10)
# Вычисляем Значение Функции Потерь
    loss = -np.sum(y * np.log(p) + (1 -y) * np.log(1-p))
    return loss / n

y = np.array([0, 0, 1, 1])

p = np.array([0.2, 0.8, 1, 0.6])

log_loss(y, p)

Кстати, в sklearn log loss можно вычислить с помощью функции log_loss из модуля metrics. Сравним библиотечный результат с нашей реализацией:

from sklearn import metrics

metrics.log_loss(y, p)

Поиск Параметров

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

Формально она записывается как:

z = x \cdot \beta = \sum_{j=0}^n \beta_j \cdot x_{j} p = \sigma(z) = \frac{1}{1 + e^{-z}} L(\beta) = - \frac{1}{n} \sum_{i=1}^n y_i \cdot log(p_{i}) + (1 - y_i) \cdot log(1 - p_{i}) L(\beta) \rightarrow \min_{\beta}

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

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

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

\beta^{(k+1)} = \beta^{(k)} - \eta \nabla L(\beta^{(k)})

Вывод градиента для логистической функции потерь здесь или здесь (раздел "Логистическая регрессия/Вывод формулы градиента")

Приведем только конечную формулу в векторно-матричной форме:

\nabla L(\beta) = - \frac{1}{n} X^T(\sigma(X \cdot \beta) - y)

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

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

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

Регуляризация

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

Реализация регуляризации для логистической регрессии в sklearn немного отличается от регуляризации для линейной регрессии. Общий вид функции потерь с регуляризацией выглядит следующим образом:

\tilde{L}(\beta) = C \cdot L(\beta) + ||\beta||_p

где

• C — коэффициент, обратный коэффициенту регуляризации. Чем больше C, тем меньше «сила» регуляризации.
• ||\beta||_p = \sum_{j=0}^m|\beta_j|^p - норма вектора весов порядка p

При $L_1$-регуляризации получим:

\tilde{L}(\beta) = C \cdot L(\beta) + \sum_{i=0}^m|\beta_i|

При $L_2$-регуляризации получим:

\tilde{L}(\beta) = C \cdot L(\beta) + \sum_{i=0}^m(\beta_i)^2 \tilde{L}(\beta) \rightarrow \min_{\beta}

Доказательство Идентичности Обычной регуляризации *

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

\tilde{L}(\beta) = L(\beta) + \alpha ||\beta||_p

где \alpha - коэффициент регуляризации.

В нашем же случае:

\tilde{L}(\beta) = C \cdot L(\beta) + \sum_{i=0}^m(\beta_i)^2

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

Можно разделить всё соотношение на некоторый коэффициент C > 0:

\frac{1}{C} ( C \cdot L(\beta) + ||\beta||_p) = L(\beta) + \frac{1}{C} ||\beta||_p

В математическом анализе доказывается, что точка минимума функций f(x) и \frac{f(x)}{C} - это одна и та же точка.

То есть нет разницы между поиска минимума функции L(\beta) и \frac{1}{C} L(\beta).

Тогда, если обозначить \alpha = \frac{1}{C}, то мы получим:

\tilde{L}(\beta) = L(\beta) + \frac{1}{C} ||\beta||_p =L(\beta) + \alpha ||\beta||_p

Вывод: коэффициент C — это коэффициент, обратный коэффициенту регуляризации \alpha.

Реализация В Sklearn

Посмотрим на реализацию логистической регрессии в sklearn.

Примеры кастомной реализации логистической регрессии вы можете найти здесь (в ООП) и здесь (в виде функций).

Заранее позаботимся о валидации модели и разделим выборку на тренировочную и тестовую:

cancer_data = datasets.load_breast_cancer(as_frame=True)

X = cancer_data.data
y = cancer_data.target

print(X.shape)
print(y.shape)

X.head()

from sklearn import model_selection

#Разделяем выборку на тренировочную и тестовую в соотношении 80/20
#Устанавливаем random_state для воспроизводимости результатов

X_train, X_test, y_train, y_test = model_selection.train_test_split(
    X, y,
    test_size=0.2,
    random_state=1,
    stratify=y #стратификация
)
#Выводим результирующие размеры таблиц
print('Train:', X_train.shape, y_train.shape)
print('Test:', X_test.shape, y_test.shape)

Проверим стратификацию данных:

y_train.value_counts(normalize=True)

y_test.value_counts(normalize=True)

X_train.head()

Логистическая регрессия — линейная модель, поэтому она находится в уже знакомом нам модуле linear_model из библиотеки sklearn.

from sklearn import linear_model #линейные модели

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

Основные параметры Logistic Regression:

• random_state — число, на основе которого происходит генерация случайных чисел.
• penalty — метод регуляризации. Возможные значения:
  • 'l1' — $L_1$-регуляризация;
  • 'l2' — $L_2$-регуляризация (используется по умолчанию);
  • 'elasticnet' — эластичная сетка ($L_1$+L_2);
  • None — отсутствие регуляризации.
• C — коэффициент обратный коэффициенту регуляризации. Чем меньше C, тем сильнее регуляризация.
• solver — численный метод оптимизации функции потерь logloss, может быть:
  • 'sag' — стохастический градиентный спуск (требует масштабирования);
  • 'saga' — модификация предыдущего, которая поддерживает работу с негладкими функциями (требует масштабирования);
  • 'newton-cg' — метод Ньютона с модификацией сопряжённых градиентов (не требует масштабирования);
  • 'lbfgs' — метод Бройдена — Флетчера — Гольдфарба — Шанно (не требует масштабирования; используется по умолчанию, так как из всех методов теоретически обеспечивает наилучшую сходимость);
  • 'liblinear' — метод покоординатного спуска (не требует масштабирования).
• max_iter — максимальное количество итераций, выделенных на сходимость.

Для обучения модели необходимо просто вызвать метод fit(), передав в него матрицу наблюдений X и вектор правильных ответов y.

# Создаём Объект Класса LogisticRegression
log_reg = linear_model.LogisticRegression(
    random_state=42,
    max_iter=10000,
    #solver='liblinear'
)
# Обучаем Модель, Минимизируя Logloss С Помощью Численных Методов
log_reg.fit(X_train, y_train)
# Выводим Результирующие Коэффициенты

Коэффициенты уравнения гиперплоскости:

• .intercept_ - свободный член (коэффициент \beta_0)
• .coef_ - коэффициенты \beta_1, \beta_2, ..., \beta_m

log_reg.intercept_

log_reg.coef_

# Составляем Таблицу Из Признаков И Их Коэффициентов
coef_df = pd.DataFrame({'Features': X_train.columns, 'Coefficients': log_reg.coef_.reshape(-1)})
# Составляем Строчку Таблицы Со Свободным Членом
intercept_df = pd.DataFrame({'Features': ['Intercept'], 'Coefficients': log_reg.intercept_.reshape(-1)})

coef_df = pd.concat([intercept_df, coef_df], ignore_index=True)
coef_df

Предсказать вероятности принадлежности к классам можно с помощью метода predict_proba:

log_reg.predict_proba(X_test).round(2)[-10: ]

Сами классы объектов с помощью метода predict:

# Определяем Классы На Тренировочной Выборке
y_train_hat = log_reg.predict(X_train)
# Определяем Классы На Тестовой Выборке
y_test_hat = log_reg.predict(X_test)

Оценим качество полученного решения по метрике accuracy:

from sklearn.metrics import accuracy_score

print('Train accuracy score: {:.2f}'.format(accuracy_score(y_train, y_train_hat)))
print('Test accuracy score: {:.2f}'.format(accuracy_score(y_test, y_test_hat)))

Преимущества И Недостатки Логистической Регрессии

Преимущества:

✔️ Простота интерпретации:

Модель имеет простую и понятную геометрическую интерпретацию, которую можно донести заказчику.

✔️ Вычислительная эффективность:

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

✔️ Устойчивость к переобучению:

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

Недостатки:

⛔️ Требование линейной-разделимости:

Главное предположение логистической регрессии заключается в том, что классы являются линейно-разделимыми. Как правило, в реальных задачах это не так.

⛔️ Неустойчивость к выбросам:

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

⛔️ Работа с категориальными признаками:

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

⛔️ Неустойчивость к высоким размерностям:

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

Полиномиальная Логистическая Регрессия

Общая Идея

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

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

Идея бинарной полиномиальной логистической регрессии (binary polynomial logistic regression) заключается в том, чтобы использовать полином внутри сигмоиды и соответственно создать нелинейную границу между двумя классами.

Например, уравнение полинома 2-ой степени для случая 2-ух признаков будет иметь вид:

z = \beta_0 + \beta_1 x_{1} + \beta_2 x_{1}^2 + \beta_3 x_{2} + \beta_4 x_{2}^2 + \beta_5 x_{1} x_{2}

При этом сам принцип построения модели логистической регрессии никак не меняется.

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

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

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

Реализация В Sklearn

Смотрим пример:

from sklearn import preprocessing

#Создаём генератор полиномиальных признаков 2 степени
poly = preprocessing.PolynomialFeatures(degree=2, include_bias=False)
poly.fit(X_train)


#Генерируем полиномиальные признаки для тренировочной выборки
X_train_poly = poly.transform(X_train)
#Генерируем полиномиальные признаки для тестовой выборки
X_test_poly = poly.transform(X_test)

X_train_poly.shape, X_test_poly.shape

# Создаём Объект Класса LogisticRegression
poly_log_reg = linear_model.LogisticRegression(
    C=0.5, # коэффициент регуляризации
    penalty='l1', # регуляризация L1
    solver='liblinear', # численный метод оптимизации
    random_state=42,
    max_iter=10000
)
# Обучаем Модель, Минимизируя Logloss С Помощью Численных Методов
poly_log_reg.fit(X_train_poly, y_train)


# Определяем Классы На Тренировочной Выборке
y_train_hat = poly_log_reg.predict(X_train_poly)
# Определяем Классы На Тестовой Выборке
y_test_hat = poly_log_reg.predict(X_test_poly)


print('Train accuracy score: {:.2f}'.format(accuracy_score(y_train, y_train_hat)))
print('Test accuracy score: {:.2f}'.format(accuracy_score(y_test, y_test_hat)))

Оценка Качества Классификации

Ранее мы с вами использовали метрику accuracy для оценки качества классификации. Но в задачах классификации часто используется и другие метрики. Познакомимся с ними.

Ошибка I И II Рода С Точки Зрения Классификация

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

Будем считать класс 1 положительным исходом или наличием эффекта (positive), а класс 0 — отрицательным (negative).

В нашем примере positive исходами считаются злокачественные опухоли, а negative - доброкачественные.

Тогда у нас есть два случая, в которых мы можем допустить ошибку:

• Ошибка I рода ($\alpha$-ошибка): ложноположительный результат. То есть мы предсказали наличие эффекта, а эффекта на самом деле нет.
  Для нашего примера - предсказали, что пациент болен раком, хотя это не так.
• Ошибка II рода ($\beta$-ошибка): ложноотрицательный результат. То есть мы предсказали отсутствие эффекта, а эффект на самом деле есть.
  Для нашего примера - предсказали, что опухоль доброкачественная, хотя на самом деле она злокачественная.

Подробнее про статистические ошибки здесь.

Ошибки I и II рода — это поможет нам понять суть метрик классификации. Давайте перейдём к ним.

Confusion Matrix

Матрица ошибок (confusion matrix) показывает все возможные исходы совпадения и несовпадения предсказания модели с действительностью. Используется для расчёта других метрик.

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

• Истинно положительные (True Positive, TP) — это объекты, обозначенные моделью как класс 1 и действительно принадлежащие к классу 1.
• Ложноположительные (False Positive, FP) — это объекты, обозначенные моделью как класс 1, но в действительности принадлежащие к классу 0. То есть это объекты, для которых модель совершила ошибку I рода.
• Истинно отрицательные (True Negative, TN) — это объекты, обозначенные моделью как класс 0 и действительно принадлежащие к классу 0.
• Ложноотрицательные (False Negative, FN) — это объекты, обозначенные моделью как класс 0, но в действительности принадлежащие к классу 1. То есть это объекты, для которых модель совершила ошибку II рода.

Для случая бинарной классификации матрица ошибок будет следующей:

| | $\large \widehat{y}=0$| \large \widehat{y}=1 | | | | | | \large y=0 | 90 | 10 | | \large y=1 |5 | 5 |

Тогда accuracy:

\large accuracy = \frac{5 + 90}{5 + 90 + 10 + 5} = 0.864

Однако представим, что мы построили классификатор, который просто предсказывает все письма как «не спам», то есть True Negative = 100, False Negative = 10, True Positive = 0, False Positive = 0.

Матрица ошибок будет иметь вид:

| | $\large \widehat{y}=0$| \large \widehat{y}=1 | | ||||| | | | |--- | | \large i=1 | 0.85 | 0.14 | 0.01 | | \large i=2 |0.1 | 0.8 |0.1 | | \large i=3 |0.5 | 0.25 |0.25 | | \large i=4 |1 | 0 |0 |

Тогда кросс-энтропию можно переписать в матричном виде:

L(\beta) = \frac{1}{N} \sum (Y \circ P)

где \circ - поэлементное умножение матриц с последующим сложением всех элементов.

Смотрим на реализацию:

def cross_entropy(y_true, y_pred):
# Выполняем One-Hot Кодировку Для Y
  y_true = np.eye(y_pred.shape[1])[y_true]
# Зажимаем Значения Вероятностей, Чтобы Не Получить Логарифм От 0
  y_pred = np.clip(y_pred, 1e-10, 1 - 1e-10)
# Рассчитываем Значение Функции Потерь
  loss = - np.sum(y_true * np.log(y_pred)) / len(y_true)
  return loss

# Истинные Метки Классов
y_true = np.array([0, 1, 2])

# Предсказанные Вероятности
y_pred = np.array([[0.8, 0.1, 0.1],
                   [0.1, 0.8, 0.1],
                   [0.1, 0.1, 0.8]])

# Считаем Кросс-энтропию
loss = cross_entropy(y_true, y_pred)
print(loss)

В sklearn кросс-энтропия считается с помощью той же самой функции log_loss из модуля metrics:

from sklearn.metrics import log_loss

# Истинные Метки Классов
y_true = np.array([0, 1, 2])

# Предсказанные Вероятности
y_pred = np.array([[0.8, 0.1, 0.1],
                   [0.1, 0.8, 0.1],
                   [0.1, 0.1, 0.8]])

log_loss(y_true, y_pred)

Поиск Параметров

Теперь, у нас есть почти все компоненты, чтобы обобщить логистическую регрессию на мультиклассовый случай.

Последний штрих - нужно лишь определить то, как мы будем искать заветные параметры \beta.

Опять же для поиска минимума можно использовать алгоритм градиентного спуска:

\beta^{(k+1)} = \beta^{(k)} - \eta \nabla L(\beta^{(k)})

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

Ниже представлена конечная формула в матричном виде:

\nabla L(\beta) = -\frac{1}{N}X^T(Y - P)

где

• X - матрица наблюдений размера N \times (M+1) (с добавочным столбцом из 1)
• Y - матрица целевой переменной в One-Hot кодировке размера N \times K
• P - матрица вероятностей размера N \times K, каждый элемент которой вычисляется как

P_{ik} = σ(X \beta_k) =\frac{e^{X \beta_k}}{\sum_{j=0}^{K-1} e^{X \beta_j}}

Реализация В Sklearn

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

import pandas as pd
from sklearn import linear_model
from sklearn import metrics
from sklearn import preprocessing
import seaborn as sns

# Загружаем Данные О Пингвинах Из Библиотеки Seaborn
penguins_data = sns.load_dataset('penguins')

# Выводим Первые Пять Строк Датасета
penguins_data.head()

penguins_data.info()

Описание данных:

• species — класс пингвина ('Adelie', 'Chinstrap', 'Gentoo'), целевой признак;
• island — остров, на котором живёт пингвин ('Torgersen', 'Biscoe', 'Dream');
• bill_length_mm — длина клюва в миллиметрах;
• bill_depth_mm — толщина клюва в миллиметрах;
• flipper_length_mm — длина крыльев;
• body_mass_g — масса;
• sex — пол ('Male', 'Female').

Наша цель — предсказать класс пингвина.

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

# Удаляем Пропущенные Значения Из Датасета
penguins_data = penguins_data.dropna()

# Разделяем Данные На Признаки И Целевую Переменную
X = penguins_data.drop('species', axis=1)
y = penguins_data['species']

# Выводим Первые Пять Строк Признаков
X.head()

Посмотрим на распределение классов в данных:

px.histogram(penguins_data, x='species', height=300)

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

px.scatter_matrix(
    penguins_data,
    dimensions=['bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'body_mass_g'],
    color='species'
)

Данные содержат строковые категориальные столбцы — island и sex. Логистическая регрессия не умеет работать со строковыми значениями. Необходимо произвести кодирование категориальных признаков.

В случае, если вы работаете с датасетом единожды и не собираетесь применять кодирование к иным данным, то можно воспользоваться простейшим вариантом кодирования категорий - функцией get_dummies() из библиотеки pandas.

Данная функция позволяет выполнять One Hot-кодирование в одно действие. Однако обладает неприятным недостатком - кодировку нельзя запомнить и применить на другие наборы данных.

X_dummies = pd.get_dummies(X)
X_dummies.head()

Создаём модель логистической регрессии. Так как классов в наших данных более 2-ух, то значение параметра multi_class выставляем на 'multinomial' (мультиклассовая классификация).

#Создаём объект класса LogisticRegression
log_reg = linear_model.LogisticRegression(
    multi_class='multinomial', #мультиклассовая классификация
    max_iter=1000, #количество итераций, выделенных на сходимость
    random_state=42 #генерация случайных чисел
)

#Обучаем модель
log_reg.fit(X_dummies, y)
#Делаем предсказание вероятностей
y_pred_proba = np.round(log_reg.predict_proba(X_dummies), 2)
#Делаем предсказание класса
y_pred = log_reg.predict(X_dummies)

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

#Создаём DataFrame из вероятностей
y_pred_proba_df = pd.DataFrame(
    y_pred_proba,
    columns=['Adelie', 'Chinstrap', 'Gentoo']
)
#Создаём DataFrame из предсказанных классов
y_pred_df = pd.DataFrame(
    y_pred,
    columns=['Predicted Class']
)
#Объединяем таблицы по вертикальной оси
y_df = pd.concat([y_pred_proba_df, y_pred_df], axis=1)
#Выбираем пять случайных строк
y_df.sample(5, random_state=2)

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

Например, для строки под номером 7 логистическая регрессия предсказала три вероятности: 0.76 — вероятность принадлежности к классу Adelie, 0.24 — к классу Chinstrap и 0 — к классу Gentoo. На основе этих вероятностей было сделано предсказание и модель отнесла пингвина в строке 7 к классу Adelie.

Давайте посмотрим, как в таком случае будет выглядеть отчёт о метриках:

print(metrics.classification_report(y, y_pred))

Для мультиклассовой классификации к отчёту просто добавляется новая строка, соответствующая третьему классу.

Из отчёта видно, что наша модель идеально решила задачу классификации (все метрики равны 1), то есть классы оказались линейно разделимыми.

Резюме

• Логистическая регрессия - один из основных алгоритмов машинного обучения, используемых для классификации
• Логистическая регрессия - это регрессия на вероятности принадлежности к классам объектов. Ключевая идея модели заключается в использовании сигмоидальной функции поверх классической линейной модели, позволяющей сводить выход к вероятностям.
• Геометрическая интерпретация логистической регрессии - разделяющая плоскость (гиперплоскость) в пространстве признаков
• Для поиска параметров модели используется логистическая функция потерь (log loss), получение которой основывается на статистическом методе максимального правдоподобия
• Поиск параметров осуществляется численными методами, в том числе методом градиентного спуска
• Логистическая регрессия как и все линейные модели допускает использование регуляризации
• Расширением модели является ее полиномиальная версия, построенная путем генерации полиномиальных признаков в данных. Такая модель способна к классификации линейно неразделимых объектов, однако имеет высокую склонность к переобучению при увеличении степени полинома
• Модель логистической регрессии легко адаптируется на случай многоклассовой классификации путем замены сигмоиды на функцию softmax и расширением логистической функции до кросс-энтропии. Принцип обучения при этом не меняется.
• Качество классификации может быть измерено с помощью разных метрик, наиболее популярные: accuracy, precision, recall и $F_1$-мера

Дополнительные Материалы

• Лекции по машинному обучению от Дмитрия Макарова (есть теория и базовая практика)
• Лекции по машинному обучению от Cornel Tech (английский, есть базовый код, МНОГО формул с выводом и доказательством)
• Учебник по машинному обучению от Яндекса (0 кода, но много математики и красивых картинок)

Связанные заметки

• Нейронные сети — персептрон и логистическая регрессия используют ту же идею линейной комбинации признаков и весов
• Обратное распространение ошибки — градиентный спуск и цепное правило обобщают идею оптимизации параметров на многослойные сети
• Принципы дискриминантного анализа — статистическая формализация той же задачи классификации; обучающая выборка + матрица "объект-свойство" = постановка задачи ML
• Функция потерь — C → min в ДА и MSE → min в МНК — один и тот же принцип минимизации ошибки классификации
• Постановка задачи снижения размерности признакового пространства — PCA как предобработка перед линейными моделями
• ВР - основы — AR-модель — это линейная регрессия, где регрессоры — лаговые значения той же переменной; оценивается МНК
• Логические операции над нечеткими множествами — нечёткие границы классов vs линейная разделимость в пространстве признаков