Лабораторное занятие 12¶

В этой лабораторной работе мы познакомимся с инструментами scikit-learn для решения задач регрессии. На этом семинаре мы познакомимся с моделью для обучения с учителем: LinearRegression.

Задача одномерной линейной регрессии¶

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

В одномерной линейной регрессии прогнозируемое значение определяется как: $$y(x) = b + k \cdot x$$

где $k, b$ - подбираемые в ходе обучения модели коэффициенты.

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

"Строгое" математическое обсуждение¶

Пускай есть набор из $y_{i}, i = \overline{1, N}$, который характеризуется набором признаков $x_{i}, i = \overline{1, N}$. Набор признаков можно представить в виде таблицы с размерами $N \times 2$. Такую таблицу обозначим через $X$: $$ X = \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \dots \\ 1 & x_N \\ \end{pmatrix} $$ Каждому признаку можно сопоставить целевое значение $y_i, i = \overline{1, N}$. Тогда можно задаться целью решить следующую систему уравнений: $$X w = y,$$ где w — вектор параметров модели с $2$ значениями. Если же явно расписать матричное уравнение, то получим $$ \begin{cases} 1 w_0 + x_1 w_1 = y_1 \\ 1 w_0 + x_2 w_1 = y_1 \\ \dots \\ 1 w_0 + x_N w_1 = y_N \\ \end{cases} $$ Получается система линейных уравнений. Данная система уравнений может не иметь точного решения. Соотвественно, мы хотим подобрать параметры $w_0, w_1$ таким образом, чтобы оказаться как можно ближе к целевым значениям $y$. Задачу можно переформулировать: провести прямую $y = w_0 + w_1 x$, которая лучше всего проходит через набор точек $(x_i, y_i), i = \overline{1, N}$. Соотвествующую задачу мы можем решить с помощью метода наименьших квадратов (МНК). Математически МНК записывется как $$ Xw - y \rightarrow \min.$$ Сам по себе МНК минимизирует следующую функцию потерь (loss function): $$ \mathrm{MSE} = \min\limits_{w} \frac{1}{N}\sum\limits_{i}\left(a(x) - y_i\right)^2 = \min\limits_{w} \frac{1}{N}\sum\limits_{i}\left((w_1 x_i + w_0) - y_i\right)^2. $$ Решения для параметров $w_0$ и $w_1$ известны. Их можно прямо посчитать с использованием производных $$ w_1 = \frac{N \sum\limits_{i = 1}^{N} x_i y_i - \sum\limits_{i = 1}^{N} y_i \sum\limits_{i = 1}^{N} y_i}{N \sum\limits_{i = 1}^{N} x_i^2 - \left(\sum\limits_{i = 1}^{N} x_i\right)^2} $$ $$ w_0 = \frac{1}{N} \sum\limits_{i = 1}^{N} y_i - w_1 \frac{1}{N} \sum\limits_{i = 1}^{N} x_i $$

Данные про алмазы¶

Воспользуемся общедоступным набором данных про алмазы для построения линейных моделей.

Для этого загрузим данные с помощью библиотеки seaborn.

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

import pandas as pd

df = pd.read_csv("https://raw.githubusercontent.com/mwaskom/seaborn-data/refs/heads/master/diamonds.csv")

Примечание 2 (для пользователей MacOS): если у вас не работает и код выше, то воспользуйтесь кодом еще ниже:

import ssl
import pandas as pd

ssl._create_default_https_context = ssl._create_unverified_context

df = pd.read_csv("https://raw.githubusercontent.com/mwaskom/seaborn-data/refs/heads/master/diamonds.csv")

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

In [1]:

import seaborn as sns

print(f"Version of seaborn: {sns.__version__}")

df = sns.load_dataset("diamonds")
df

Version of seaborn: 0.13.2

Out[1]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	Ideal	E	SI2	61.5	55.0	326	3.95	3.98	2.43
1	0.21	Premium	E	SI1	59.8	61.0	326	3.89	3.84	2.31
2	0.23	Good	E	VS1	56.9	65.0	327	4.05	4.07	2.31
3	0.29	Premium	I	VS2	62.4	58.0	334	4.20	4.23	2.63
4	0.31	Good	J	SI2	63.3	58.0	335	4.34	4.35	2.75
...	...	...	...	...	...	...	...	...	...	...
53935	0.72	Ideal	D	SI1	60.8	57.0	2757	5.75	5.76	3.50
53936	0.72	Good	D	SI1	63.1	55.0	2757	5.69	5.75	3.61
53937	0.70	Very Good	D	SI1	62.8	60.0	2757	5.66	5.68	3.56
53938	0.86	Premium	H	SI2	61.0	58.0	2757	6.15	6.12	3.74
53939	0.75	Ideal	D	SI2	62.2	55.0	2757	5.83	5.87	3.64

53940 rows × 10 columns

В наборе данных присутствует следующая информация:

carat: масса в каратах
cut: тип граненки
color: цвет
clarity: прозрачность
depth: некоторый линейный размер, %
table: некоторый линейный размер, %
price: цена, $
x: некоторый линейный размер, мм
y: некоторый линейный размер, мм
z: некоторый линейный размер, мм

Попробуем оценить цену price алмаза по его "объему" двумя разными способами:

Посчитав коэффициенты $w_0$ и $w_1$ по известным формулам
Реализовав градиентный спуск

Подготовка данных¶

Упражнение 1. Выделите отдельно данные о цене в переменную y и отдельно вычислите "объем" алмаза $x \times y \times z$ мм$^3$ в переменную x.

In [ ]:

# Заполните код
x = 
y =

Теперь графически отобразим зависимость цены алмаза от объема.

In [3]:

from matplotlib import pyplot as plt

plt.scatter(x, y, s=10)
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.tight_layout()

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

Упражнение 2. Теперь подсчитает коэффициенты $w_0$ и $w_1$ по ранее определенным формулам.

In [ ]:

N = len(x)

# впишите формулы вычисления для w_1 и w_0
w_1 = 
w_0 = 

precise_w = (w_0, w_1)
print(f"Best fit values: w_0 = {w_0:1.3f}, w_1 = {w_1:1.3f}")

Best fit values: w_0 = -2041.479, w_1 = 46.009

Добавим нашу прямую на предыдущий график.

In [5]:

import numpy as np

x_plot = np.array([20, 500])

plt.scatter(x, y, s=10, label="data")
plt.plot(x_plot, w_0 + w_1 * x_plot, label="fit", color="C2")
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.legend()
plt.xlim((0, 2e+3))
plt.ylim((0, 2e+4))
plt.tight_layout()

По виду данных создается впечатление, что наклон мог бы быть еще круче. Скорее всего наша модель сильно отклонилась из-за аномального выброса для объема > 3500 мм$^3$ (выброс хорошо виден на предыдущем графике).

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

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

Ранее была введена функция потерь MSE. Подбором параметров в нашей линейной модели мы пытаемся минимизировать значение MSE. Попробуем выписать производную MSE по параметрам $w$ $$ \mathrm{MSE}' = \nabla_w \mathrm{MSE} = \begin{pmatrix} \dfrac{d \mathrm{MSE}}{d w_0} \\ \dfrac{d \mathrm{MSE}}{d w_1} \end{pmatrix} = \begin{pmatrix} \dfrac{2}{N}\sum\limits_{i}\left((w_1 x_i + w_0) - y_i\right) \\ \dfrac{2 w_1}{N}\sum\limits_{i}\left((w_1 x_i + w_0) - y_i\right) \end{pmatrix} $$ Получили вектор, который состоит из производных по параметрам. Если подставить численные значения параметров $w$, то получим соответствующее значение производной в точке. Если же мы хотим минимизировать MSE, то нам нужно найти такие значения $w$, что каждая производная будет зануляться. Этого можно добиться при условии $$ w^{(i + 1)} = w^{(i)} - \gamma \nabla_w\mathrm{MSE} \rightarrow \begin{pmatrix} w_0^{(i + 1)} \\ w_1^{(i + 1)} \end{pmatrix} = \begin{pmatrix} w_0^{(i)} \\ w_1^{(i)} \end{pmatrix} - \gamma \begin{pmatrix} \dfrac{d \mathrm{MSE}}{d w_0} \\ \dfrac{d \mathrm{MSE}}{d w_1} \end{pmatrix}, $$ где $\gamma$ является параметром скорости обучения и подбирается эмперически.

На подумать: какое графическое представление у вектора выше? Каков его физический (?) смысл.

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

In [ ]:

w_0 = -2000
w_1 = 30
gamma = 1e-5
derivative_step = 1e-3
n_iterations = 100
N = len(x)

loss_function = lambda x, y, w_0, w_1: ((x * w_1 + w_0 - y)**2).sum() / len(x)

loss_values = []

for i in range(n_iterations):
    # Впишите формульные значения для градиентов grad0 и grad1
    grad0 = 
    grad1 = 
    w_0 -= 
    w_1 -= 
    loss_values.append(loss_function(x, y, w_0, w_1))

approx_w = (w_0, w_1)
print(f"Gradient fit values: w_0 = {w_0:1.3f}, w_1 = {w_1:1.3f}")

Gradient fit values: w_0 = -1999.575, w_1 = 45.686

Построим график зависимости значения функции потерь от номера итерации

In [7]:

plt.plot(loss_values)
plt.xlabel("iteration")
plt.ylabel("loss function")
plt.tight_layout()

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

In [8]:

plt.scatter(x, y, s=10, label="data")
plt.plot(x_plot, approx_w[0] + approx_w[1] * x_plot, label="gradient descent", color="C1")
plt.plot(x_plot, precise_w[0] + precise_w[1] * x_plot, label="best fit", color="C2", linestyle="--")
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.legend()
plt.xlim((0, 2e+3))
plt.ylim((0, 2e+4))
plt.tight_layout()

Прямые накладываются друг на друга. Приближенный метод градиентного спуска дает схожий ответ с честным решением МНК.

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

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

In [ ]:

Модель линейной регрессии в sklearn¶

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

Попробуем решить задачу о цене алмаза по его объему с помощью модели из пакета scikit-learn.

In [9]:

from sklearn.linear_model import LinearRegression

model_l = LinearRegression()

У класса LinearRegression есть параметр fit_intercept: bool. Он позволяет включить или исключить из подгонки свободный параметр (в предыдущих терминах это будет параметр $w_0$).

В предыдущей ячейке мы создали модель. Теперь же попробуем ее обучить на полном наборе данных. Для этого воспользуемся методом .fit() нашей модели.

In [10]:

X_l = x.copy()
X_l = X_l.to_numpy().reshape((-1, 1))  # Метод fit(...) принимает только 2d массивы признаков

model_l.fit(X_l, y)
model_w = (model_l.intercept_, *model_l.coef_)

print(f"Model fit values: w_0 = {model_w[0]:1.3f}, w_1 = {model_w[1]:1.3f}")

Model fit values: w_0 = -2041.479, w_1 = 46.009

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

In [11]:

plt.scatter(x, y, s=10, label="data")
plt.plot(x_plot, approx_w[0] + approx_w[1] * x_plot, label="gradient descent", color="C1")
plt.plot(x_plot, precise_w[0] + precise_w[1] * x_plot, label="best fit", color="C2", linestyle="--")
plt.plot(x_plot, model_w[0] + model_w[1] * x_plot, label="linear model", color="C3", linestyle=":")
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.legend()
plt.xlim((0, 2e+3))
plt.ylim((0, 2e+4))
plt.tight_layout()

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

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

In [12]:

X_q = np.stack((x, x**2), axis=1)

model_q = LinearRegression()

model_q.fit(X_q, y)

model_q_w = (model_q.intercept_, *model_q.coef_)

print(f"Model fit values: w_0 = {model_q_w[0]:1.3f}, w_1 = {model_q_w[1]:1.3f}, w_2 = {model_q_w[2]:1.3f}")

Model fit values: w_0 = -2541.626, w_1 = 51.988, w_2 = -0.012

In [13]:

plt.scatter(x, y, s=10, label="data")
plt.plot(x_plot, precise_w[0] + precise_w[1] * x_plot, label="best fit", color="C1", linestyle="-")
plt.plot(x_plot, model_w[0] + model_w[1] * x_plot, label="linear model", color="C2", linestyle="--")
plt.plot(x_plot, model_q_w[0] + model_q_w[1] * x_plot + model_q_w[2] * x_plot**2, label="quadratic model", color="C3", linestyle=":")
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.legend()
plt.xlim((0, 2e+3))
plt.ylim((0, 2e+4))
plt.tight_layout()

Квадратичная модель незначительно отклоняется от линейной. Как же сравнить степень незначительности?

Для ответа на этот вопрос мы можем посчитать значение MSE для линейной и квадратичной моделей.

In [14]:

from sklearn.metrics import mean_squared_error

mse_l = mean_squared_error(y, model_l.predict(X_l))
mse_q = mean_squared_error(y, model_q.predict(X_q))

print(
    f"Linear approximation MSE: {mse_l:1.3f}\n"
    f"Quadratic approximation MSE: {mse_q:1.3f}"
)

Linear approximation MSE: 2955511.787
Quadratic approximation MSE: 2481050.863

Квадратичная модель дает меньшее значение. Если продолжать увеличивать степень многочлена, то значение MSE будет уменьшаться, но это может быть опасно переобучением.

Упражнение 5. Попробуйте в ячейке ниже написать модель для многочлена 4ой степени и посчитать для нее MSE.

In [ ]:

model_f = LinearRegression()

# Составьте данные для многочлена 4ой степени и вычислите MSE для этой модели
X_f = 
model_f.fit(X_f, y)
mse_f = mean_squared_error(y, model_f.predict(X_f))

model_f_w = (model_f.intercept_, *model_f.coef_)

print(f"Forth polynomial approximation MSE: {mse_f:1.3f}")

Forth polynomial approximation MSE: 2123287.961

Изобразим новую модель.

In [16]:

x_plot = np.arange(0, 1000, 10)
plt.scatter(x, y, s=10, label="data")
plt.plot(x_plot, precise_w[0] + precise_w[1] * x_plot, label="best fit", color="C1", linestyle="-")
plt.plot(x_plot, model_w[0] + model_w[1] * x_plot, label="linear model", color="C2", linestyle="--")
plt.plot(x_plot, model_f_w[0] + model_f_w[1] * x_plot + model_f_w[2] * x_plot**2 + model_f_w[3] * x_plot**3 + model_f_w[4] * x_plot**4, label="forth model", color="C3", linestyle=":")
plt.xlabel(r"Volume, mm$^3$")
plt.ylabel("Price")
plt.legend()
plt.xlim((0, 2e+3))
plt.ylim((0, 2e+4))
plt.tight_layout()

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

Задача многомерной линейной регрессии¶

Пускай есть набор из $N$ данных, который характеризуется некоторой "строкой" из $M$ признаков $x^{(j)}, j = \overline{1, M}$. Набор можно представить в виде таблицы с размерами $N \times M$. Такую таблицу обозначим через $X$: $$ X = \begin{pmatrix} x^{(1)}_1 & x^{(2)}_1 & \dots & x^{(M)}_1 \\ x^{(1)}_2 & x^{(2)}_2 & \dots & x^{(M)}_2 \\ \dots & \dots & \dots & \dots \\ x^{(1)}_N & x^{(2)}_N & \dots & x^{(M)}_N \\ \end{pmatrix} $$ Каждой строке сопоставим некоторое целевое значение $y_i, i = \overline{1, N}$. Тогда можно задаться целью решить следующую систему уравнений: $$X w = y,$$ где w — вектор параметров модели. Мы хотим подобрать параметры $w$ таким образом, чтобы оказаться как можно ближе к целевому вектору $y$. Данное "желание" можно исполнить, если найти минимум значения функции потерь (loss function). Для задачи линейной регрессии типично используется среднеквадратичная ошибка (MSE) $$ \mathrm{MSE} = \min\limits_{k_0, k_1, \dots, k_M} \frac{1}{N}\left(\left(k_{1} x^{(1)}_i + k_{2} x^{(2)}_i + \dots k_{M} x^{(M)}_i + k_0\right) - y_i\right)^2. $$

Эту формулу можно представить и в матричном виде: $$ Q(w) = N \mathrm{MSE}^2 = \min\limits_{k_0, k_1, \dots, k_M} \left((k_{1} x^{(1)}_i + k_{2} x^{(2)}_i + \dots k_{M} x^{(M)}_i + k_0) - y_i\right)^2 = (Xw - y)^T (Xw - y).$$ Для минимизации множества параметров требуется найти такое значение для вектора $w$, что производная от $Q(w)$ по вектору параметров $w$ будет равна 0. Можно обобщить на $M$-мерный случай понятие градиента: $$ \nabla_x f(x) = \left(\dfrac{d f}{d x^{(1)}}, \dots, \dfrac{d f}{d x^{(M)}}\right)^T. $$ Хитрыми махинациями и магией можно получить, что искомый вектор параметров $w$ выражается как $$ w = (X^T X)^{-1} X^T y $$

Но многомерный градиент в данной лабораторной работе реализовывать не будем... а могли.

Работа с данными большей размерности¶

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

Выберем целевой параметр (цена алмаза) и параметры модели (все оставшиеся)

In [17]:

feature_names = ["carat", "cut", "color", "clarity", "depth", "table", "x", "y", "z"]
target_name = "price"

df.head()

Out[17]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	Ideal	E	SI2	61.5	55.0	326	3.95	3.98	2.43
1	0.21	Premium	E	SI1	59.8	61.0	326	3.89	3.84	2.31
2	0.23	Good	E	VS1	56.9	65.0	327	4.05	4.07	2.31
3	0.29	Premium	I	VS2	62.4	58.0	334	4.20	4.23	2.63
4	0.31	Good	J	SI2	63.3	58.0	335	4.34	4.35	2.75

Среди данных есть строчные/категориальные. Их нужно будет преобразовать

In [18]:

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 53940 entries, 0 to 53939
Data columns (total 10 columns):
 #   Column   Non-Null Count  Dtype   
---  ------   --------------  -----   
 0   carat    53940 non-null  float64 
 1   cut      53940 non-null  category
 2   color    53940 non-null  category
 3   clarity  53940 non-null  category
 4   depth    53940 non-null  float64 
 5   table    53940 non-null  float64 
 6   price    53940 non-null  int64   
 7   x        53940 non-null  float64 
 8   y        53940 non-null  float64 
 9   z        53940 non-null  float64 
dtypes: category(3), float64(6), int64(1)
memory usage: 3.0 MB

Конвертация категориальных данных в числовые¶

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

In [19]:

from sklearn.preprocessing import LabelEncoder

encoder = LabelEncoder()

for column in df.select_dtypes(["category"]):
    df[column] = encoder.fit_transform(df[column])

df.head()

Out[19]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	2	1	3	61.5	55.0	326	3.95	3.98	2.43
1	0.21	3	1	2	59.8	61.0	326	3.89	3.84	2.31
2	0.23	1	1	4	56.9	65.0	327	4.05	4.07	2.31
3	0.29	3	5	5	62.4	58.0	334	4.20	4.23	2.63
4	0.31	1	6	3	63.3	58.0	335	4.34	4.35	2.75

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

In [20]:

sns.pairplot(df);

In [21]:

correlation_df = df[feature_names].corr()
sns.heatmap(correlation_df, vmin=-1, vmax=1, annot=True)
plt.figure()
sns.heatmap(correlation_df[np.abs(correlation_df) > 0.75], vmin=-1, vmax=1, annot=True);

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

In [22]:

X = df[feature_names].to_numpy()
y = df[target_name].to_numpy()

model_md = LinearRegression()
model_md.fit(X, y)

mse_md = mean_squared_error(model_md.predict(X), y)
print(f"Multi-dimensional approximation MSE: {mse_md:1.3f}")

Multi-dimensional approximation MSE: 1829128.158

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

In [23]:

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

model_md = LinearRegression()
model_md.fit(X_train, y_train)

mse_md = mean_squared_error(model_md.predict(X_test), y_test)
print(f"Multi-dimensional approximation MSE: {mse_md:1.3f}")

Multi-dimensional approximation MSE: 1750176.393

Среди данных могли попасться дубликаты строк.

Упражнение 6. Попробуйте проверить набор данных на уникальность значений. Если обнаружатся одинаковые строки, то удалите их из набора данных df с помощью метода .drop_duplicates(). Обучите модель линейной регрессии на отфильтрованных данных.

In [24]:

df.shape, df.drop_duplicates().shape

Out[24]:

((53940, 10), (53794, 10))

In [ ]:

df_unique =   # Вставьте сюда свой код
X = df_unique[feature_names].to_numpy()
y = df_unique[target_name].to_numpy()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

model_md = LinearRegression()
model_md.fit(X_train, y_train)

mse_md = mean_squared_error(model_md.predict(X_test), y_test)
print(f"Multi-dimensional approximation MSE: {mse_md:1.3f}")

Multi-dimensional approximation MSE: 1765296.815

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

Для более достоверной оценки качества модели используется кросс-валидация, при которой сначала выполняем разделение на обучающую и тестовую выборки. Далее обучающую часть выборки делим на k примерно (в случае если на число объектов не делится нецело на k) равных частей (это будет k-fold кросс-валидация).

На первой итерации обучаем на 1-4 частях данных, и на 5-й оцениваем качество модели.
На следующей итерации обучаем модель на 1-ой и 2-5-й частях данных и на 2-й части тестируем модель.
и т. д.

Смотрете картинку ниже.

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

Читать продолжение в источнике...

Воспользуемся методикой кросс-валидации:

In [26]:

from sklearn.model_selection import cross_validate
from sklearn.model_selection import KFold

kf = KFold(n_splits=5, random_state=None, shuffle=True)

model_cv = LinearRegression()
X = df[feature_names].to_numpy()
y = df[target_name].to_numpy()

cv_res = cross_validate(
    model_cv,
    X,
    y,
    scoring='neg_mean_squared_error',
    cv=kf,
)

print(f"Test mse errors are {cv_res['test_score']}")
print(f"Mean test mse = {cv_res['test_score'].mean()}")

Test mse errors are [-1825572.45352771 -1776426.34873347 -1851196.2680185  -2349542.61371931
 -1890815.65735265]
Mean test mse = -1938710.668270329

Метод K-ближайших соседей в задаче регрессии¶

Рассмотрим модель K-ближайших соседей (K-nearest neighbors KNN). Интуитивное понимание модели — для выбранного набора данных можно вычислить расстояние между точками. Принадлежность точки к определенному классу определяется большинством соседей этой точки.

Метод kNN можно использовать для решения задач регрессии. Алгоритм определяет значение целевого параметра по среднему значению ближайших соседей.

In [27]:

from sklearn.neighbors import KNeighborsRegressor

model_knn = KNeighborsRegressor(n_neighbors=5, metric="minkowski")

Попробуем обучить модель на имеющихся данных.

In [28]:

X = df[feature_names].to_numpy()
y = df[target_name].to_numpy()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

model_knn.fit(X_train, y_train)

mse_knn = mean_squared_error(model_knn.predict(X_test), y_test)

print(f"KNN model MSE: {mse_knn:1.3f}")

KNN model MSE: 827841.171

Оценим влияние количества соседей на MSE.

In [29]:

mse_knn = []
for n_neighbors in range(1, 20):
    model_knn = KNeighborsRegressor(n_neighbors=n_neighbors, metric="euclidean")
    model_knn.fit(X_train, y_train)

    mse_knn.append(mean_squared_error(model_knn.predict(X_test), y_test))

mse_knn = np.array(mse_knn)
plt.scatter(range(1, 20), mse_knn)
plt.xlim((0, 21))
plt.xlabel("N neighbors")
plt.ylabel("MSE")
plt.tight_layout()

print(f"Optimal number of neighbors: {mse_knn.argmin() + 1}")

Optimal number of neighbors: 5

Предварительный вывод: оптимальное количество соседей — 5.

Нормирование данных¶

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

Уменьшить влияние масштаба можно с помощью нормирования.

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

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

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

In [ ]:

import pandas as pd
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

# Вставьте код
X_train_transformed = scaler.fit_transform
X_test_transformed = scaler.transform

sns.pairplot(pd.DataFrame(X_train_transformed, columns=feature_names))

model_l_transformed = KNeighborsRegressor()

# Вставьте код
model_l_transformed.fit
y_pred_transformed =

mse_transformed = mean_squared_error(y_test, y_pred_transformed)


model_l = KNeighborsRegressor()

# Вставьте код
model_l.fit
y_pred = 

mse = mean_squared_error(y_test, y_pred)

print(f"MSE (data was transformed): {mse_transformed:<1.7f}")
print(f"MSE (no transformation):    {mse:>14.7f}")

MSE (data was transformed): 649250.5048276
MSE (no transformation):    827841.1705784

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

Для оценки качества модели линейной регрессии $f$ используют метрику $R^2$: $$ R^2 = 1 - \dfrac{ \sum\limits_{i=1}^{N}(y_i - f(x_i))^2 }{ \sum\limits_{i=1}^{N}(y_i - \overline{y})^2 }$$ Ее можно интерпретировать следующим образом:

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

Метрика $R^2$ может быть применима для оценки моделей с одинаковым количеством признаков. Если признаков разное количество, то требуется применяять друугие определения метрик. Смотрите тут.

Упражнение 8: попробуйте вычислить значение $R^2$ метрики для моделей model_l_transformed и model_l из предыдущего упражнения.

In [ ]:

from sklearn.metrics import r2_score

# Ваш код здесь

Out[ ]:

(0.9582166873855468, 0.9467232660302995)

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

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

In [ ]:

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import Normalizer
from sklearn.metrics import r2_score, mean_squared_error


dataset = load_diabetes(as_frame=True)
df = dataset.frame

print(dataset.DESCR)

.. _diabetes_dataset:

Diabetes dataset
----------------

Ten baseline variables, age, sex, body mass index, average blood
pressure, and six blood serum measurements were obtained for each of n =
442 diabetes patients, as well as the response of interest, a
quantitative measure of disease progression one year after baseline.

**Data Set Characteristics:**

:Number of Instances: 442

:Number of Attributes: First 10 columns are numeric predictive values

:Target: Column 11 is a quantitative measure of disease progression one year after baseline

:Attribute Information:
    - age     age in years
    - sex
    - bmi     body mass index
    - bp      average blood pressure
    - s1      tc, total serum cholesterol
    - s2      ldl, low-density lipoproteins
    - s3      hdl, high-density lipoproteins
    - s4      tch, total cholesterol / HDL
    - s5      ltg, possibly log of serum triglycerides level
    - s6      glu, blood sugar level

Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times the square root of `n_samples` (i.e. the sum of squares of each column totals 1).

Source URL:
https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html

For more information see:
Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) "Least Angle Regression," Annals of Statistics (with discussion), 407-499.
(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)

In [34]:

df.head()

Out[34]:

	age	sex	bmi	bp	s1	s2	s3	s4	s5	s6	target
0	0.038076	0.050680	0.061696	0.021872	-0.044223	-0.034821	-0.043401	-0.002592	0.019907	-0.017646	151.0
1	-0.001882	-0.044642	-0.051474	-0.026328	-0.008449	-0.019163	0.074412	-0.039493	-0.068332	-0.092204	75.0
2	0.085299	0.050680	0.044451	-0.005670	-0.045599	-0.034194	-0.032356	-0.002592	0.002861	-0.025930	141.0
3	-0.089063	-0.044642	-0.011595	-0.036656	0.012191	0.024991	-0.036038	0.034309	0.022688	-0.009362	206.0
4	0.005383	-0.044642	-0.036385	0.021872	0.003935	0.015596	0.008142	-0.002592	-0.031988	-0.046641	135.0

In [ ]:

Задание 2. В качестве задачи на линейную регрессию возьмем классический набор данных калифорнийской недвижимости (California housing). Вам предоставляется информация о характеристиках домов в одном из калифорнийском уезде. Попробуйте определить влияние параметров на целевую метрику "MedHouseVal".

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

In [ ]:

from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import Normalizer
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.neighbors import KNeighborsRegressor


dataset = fetch_california_housing(as_frame=True)

df = dataset.frame

print(dataset.DESCR)

.. _california_housing_dataset:

California Housing dataset
--------------------------

**Data Set Characteristics:**

:Number of Instances: 20640

:Number of Attributes: 8 numeric, predictive attributes and the target

:Attribute Information:
    - MedInc        median income in block group
    - HouseAge      median house age in block group
    - AveRooms      average number of rooms per household
    - AveBedrms     average number of bedrooms per household
    - Population    block group population
    - AveOccup      average number of household members
    - Latitude      block group latitude
    - Longitude     block group longitude

:Missing Attribute Values: None

This dataset was obtained from the StatLib repository.
https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html

The target variable is the median house value for California districts,
expressed in hundreds of thousands of dollars ($100,000).

This dataset was derived from the 1990 U.S. census, using one row per census
block group. A block group is the smallest geographical unit for which the U.S.
Census Bureau publishes sample data (a block group typically has a population
of 600 to 3,000 people).

A household is a group of people residing within a home. Since the average
number of rooms and bedrooms in this dataset are provided per household, these
columns may take surprisingly large values for block groups with few households
and many empty houses, such as vacation resorts.

It can be downloaded/loaded using the
:func:`sklearn.datasets.fetch_california_housing` function.

.. rubric:: References

- Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions,
  Statistics and Probability Letters, 33:291-297, 1997.

In [ ]:

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

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

Дерево решений¶

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

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

Данные про алмазы¶

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

Для этого загрузим данные с помощью библиотеки seaborn.

In [36]:

import seaborn as sns

print(f"Version of seaborn: {sns.__version__}")

df = sns.load_dataset("diamonds")
df

Version of seaborn: 0.13.2

Out[36]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	Ideal	E	SI2	61.5	55.0	326	3.95	3.98	2.43
1	0.21	Premium	E	SI1	59.8	61.0	326	3.89	3.84	2.31
2	0.23	Good	E	VS1	56.9	65.0	327	4.05	4.07	2.31
3	0.29	Premium	I	VS2	62.4	58.0	334	4.20	4.23	2.63
4	0.31	Good	J	SI2	63.3	58.0	335	4.34	4.35	2.75
...	...	...	...	...	...	...	...	...	...	...
53935	0.72	Ideal	D	SI1	60.8	57.0	2757	5.75	5.76	3.50
53936	0.72	Good	D	SI1	63.1	55.0	2757	5.69	5.75	3.61
53937	0.70	Very Good	D	SI1	62.8	60.0	2757	5.66	5.68	3.56
53938	0.86	Premium	H	SI2	61.0	58.0	2757	6.15	6.12	3.74
53939	0.75	Ideal	D	SI2	62.2	55.0	2757	5.83	5.87	3.64

53940 rows × 10 columns

В наборе данных присутствует следующая информация:

carat: масса в каратах
cut: тип граненки
color: цвет
clarity: прозрачность
depth: некоторый линейный размер, %
table: некоторый линейный размер, %
price: цена, $
x: некоторый линейный размер, мм
y: некоторый линейный размер, мм
z: некоторый линейный размер, мм

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

In [37]:

feature_names = ["carat", "cut", "color", "clarity", "depth", "table", "x", "y", "z"]
target_name = "price"

df.head()

Out[37]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	Ideal	E	SI2	61.5	55.0	326	3.95	3.98	2.43
1	0.21	Premium	E	SI1	59.8	61.0	326	3.89	3.84	2.31
2	0.23	Good	E	VS1	56.9	65.0	327	4.05	4.07	2.31
3	0.29	Premium	I	VS2	62.4	58.0	334	4.20	4.23	2.63
4	0.31	Good	J	SI2	63.3	58.0	335	4.34	4.35	2.75

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

In [38]:

from sklearn.preprocessing import LabelEncoder

encoder = LabelEncoder()

for column in df.select_dtypes(["category"]):
    df[column] = encoder.fit_transform(df[column])

df.head()

Out[38]:

	carat	cut	color	clarity	depth	table	price	x	y	z
0	0.23	2	1	3	61.5	55.0	326	3.95	3.98	2.43
1	0.21	3	1	2	59.8	61.0	326	3.89	3.84	2.31
2	0.23	1	1	4	56.9	65.0	327	4.05	4.07	2.31
3	0.29	3	5	5	62.4	58.0	334	4.20	4.23	2.63
4	0.31	1	6	3	63.3	58.0	335	4.34	4.35	2.75

Обучим дерево на всех возможных признаках.

In [39]:

from sklearn import tree

model_tree = tree.DecisionTreeRegressor(random_state=1, max_depth=3, criterion="squared_error")

In [40]:

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X = df[feature_names].to_numpy()
y = df[target_name].to_numpy()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

model_tree.fit(X_train, y_train)

y_pred = model_tree.predict(X_test)

mean_squared_error(y_test, y_pred)

Out[40]:

1817470.888303215

Дерево решений показывает значительно меньшее значение MSE (для линейной регрессии значение данной метрики/loss-функции было равно $2\cdot 10^6$). Подумайте как это объяснить.

Попробуем визуализировать дерево в текстовом формате.

In [41]:

print(tree.export_text(model_tree))

|--- feature_0 <= 1.00
|   |--- feature_7 <= 5.53
|   |   |--- feature_7 <= 4.99
|   |   |   |--- value: [788.69]
|   |   |--- feature_7 >  4.99
|   |   |   |--- value: [1695.89]
|   |--- feature_7 >  5.53
|   |   |--- feature_0 <= 0.89
|   |   |   |--- value: [2726.66]
|   |   |--- feature_0 >  0.89
|   |   |   |--- value: [3961.79]
|--- feature_0 >  1.00
|   |--- feature_7 <= 7.19
|   |   |--- feature_3 <= 3.50
|   |   |   |--- value: [5149.85]
|   |   |--- feature_3 >  3.50
|   |   |   |--- value: [7467.35]
|   |--- feature_7 >  7.19
|   |   |--- feature_7 <= 7.86
|   |   |   |--- value: [10937.97]
|   |   |--- feature_7 >  7.86
|   |   |   |--- value: [14914.64]

Из текстового формата видно, что первый отбор происходит по самой первой характеристике алмаза - массе алмаза в каратах. После этого выбор осуществляется по линейной характеристике x и так далее. Видно, что в дереве присутствует всего 8 конечных значений, на которых у нас и получается значение MSE выше. Если бы дерево решений имело бы больше конечных решений (листьев), то MSE было бы еще меньше! Но подход с большим числом листьев может привести к переобучению.

Попробуем визуализировать дерево в графическом формате.

In [42]:

from matplotlib import pyplot as plt

plt.figure(figsize=(10, 6))
tree.plot_tree(model_tree);

Если не удается рассмотреть картинку выше, то

Сделайте ее больше или
Сохрание в pdf-формате и откройте сохраненный файл.

Рассмотрим влияние глубины дерева на значение MSE.

In [43]:

mses = []
for depth in range(1, 20):
    model_tree = tree.DecisionTreeRegressor(random_state=1, max_depth=depth, criterion="squared_error")
    X = df[feature_names].to_numpy()
    y = df[target_name].to_numpy()
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)
    model_tree.fit(X_train, y_train)
    y_pred = model_tree.predict(X_test)
    mses.append(mean_squared_error(y_test, y_pred))

plt.scatter(range(1, 20), mses)
plt.xlabel("Tree depth")
plt.ylabel("MSE")
plt.tight_layout()

Видно, что чем больше глубина дерева, тем меньше значение MSE. Но видно еще и "насыщение" нашего алгоритма: уже после глубины 10 MSE не так быстро убывает, а для глубин >12 будто бы растет. Это и является косвенным признаком "переобучения" нашего алгоритма.

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

In [44]:

import numpy as np

n_samples = 100
np.random.seed(1)

X = np.linspace(0, 20, n_samples, endpoint=False).reshape((-1, 1))
y = np.sin(X) + 0.01 * X + 0.1 * np.random.normal(X)
y = y.reshape(-1)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

plt.scatter(X, y, color="C3", facecolor="none", s=5)
plt.ylabel("target")
plt.xlabel("feature")

for depth in (1, 3, 9):
    model_tree = tree.DecisionTreeRegressor(max_depth=depth, random_state=1)
    model_tree.fit(X_train, y_train)
    y_pred = model_tree.predict(sorted(X_test))
    plt.plot(sorted(X_test), y_pred, label=f"depth = {depth}, MSE = {mean_squared_error(y_test, y_pred):1.3f}", linestyle="--")

plt.legend()
plt.tight_layout()

Для глубины дерева 1 мы видим ступенчатую зависимость: если значение признака x меньше 12, то ответ модели примерно 0.6. Если больше 12, то что-то близкое к 2.0.

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

Для глубины 9 модель начинает подтсраиваться под особенности в данных.

Еще один примечательный факт: метрика MSE получается минимальной для случая глубины дерева 1. Но и это можно объяснить, попробуйте подумать как.