Все курсы > Оптимизация > Занятие 4 (часть 2)
В первой части мы рассмотрели теорию линейной регрессии. Теперь перейдем к практике.
Откроем ноутбук к этому занятию⧉
Сквозной пример
Данные и постановка задачи
Обратимся к хорошо знакомому нам датасету недвижимости в Бостоне.
|
1 |
boston = pd.read_csv('/content/boston.csv') |
При этом нам нужно будет решить две основные задачи:
Задача 1. Научиться оценивать качество модели не только с точки зрения метрики, но и исходя из рассмотренных ранее допущений модели. Эту задачу мы решим в три этапа.
- Этап 1. Построим базовую (baseline) модель линейной регрессии с помощью класса LinearRegression библиотеки sklearn и оценим, насколько выполняются рассмотренные выше допущения.
- Этап 2. Попробуем изменить данные таким образом, чтобы модель в большей степени соответствовала этим критериям.
- Этап 3. Обучим еще одну модель и посмотрим, как изменится результат.
Задача 2. С нуля построить модель множественной линейной регрессии и сравнить прогноз с результатом, полученным при решении первой задачи. При этом обучение модели мы реализуем двумя способами, а именно, через:
- Метод наименьших квадратов
- Метод градиентного спуска
Разделение выборки
Мы уже не раз говорили про важность разделения выборки на обучаущую и тестовую части. Сегодня же, с учетом того, что нам предстоит изучить много нового материала, мы опустим этот этап и будем обучать и тестировать модель на одних и тех же данных.
Исследовательский анализ данных
Теперь давайте более внимательно посмотрим на имеющиеся у нас данные. Как вы вероятно заметили, все признаки в этом датасете количественные, за исключением переменной CHAS.
|
1 |
boston.info() |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
<class 'pandas.core.frame.DataFrame'> RangeIndex: 506 entries, 0 to 505 Data columns (total 14 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 CRIM 506 non-null float64 1 ZN 506 non-null float64 2 INDUS 506 non-null float64 3 CHAS 506 non-null float64 4 NOX 506 non-null float64 5 RM 506 non-null float64 6 AGE 506 non-null float64 7 DIS 506 non-null float64 8 RAD 506 non-null float64 9 TAX 506 non-null float64 10 PTRATIO 506 non-null float64 11 B 506 non-null float64 12 LSTAT 506 non-null float64 13 MEDV 506 non-null float64 dtypes: float64(14) memory usage: 55.5 KB |
|
1 2 |
# мы видим, что переменная CHAS категориальная boston.CHAS.value_counts() |
|
1 2 3 |
0.0 471 1.0 35 Name: CHAS, dtype: int64 |
Посмотрим на распределение признаков с помощью boxplots.
|
1 2 3 |
plt.figure(figsize = (10, 8)) sns.boxplot(data = boston.drop(columns = ['CHAS', 'MEDV'])) plt.show() |
Посмотрим на распределение целевой переменной.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
def box_density(x): # создадим два подграфика f, (ax_box, ax_kde) = plt.subplots(nrows = 2, # из двух строк ncols = 1, # и одного столбца sharex = True, # оставим только нижние подписи к оси x gridspec_kw = {'height_ratios': (.15, .85)}, # зададим разную высоту строк figsize = (10,8)) # зададим размер графика # в первом подграфике построим boxplot sns.boxplot(x = x, ax = ax_box) ax_box.set(xlabel = None) # во втором - график плотности распределения sns.kdeplot(x, fill = True) # зададим заголовок и подписи к осям ax_box.set_title('Распределение переменной', fontsize = 17) ax_kde.set_xlabel('Переменная', fontsize = 15) ax_kde.set_ylabel('Плотность распределения', fontsize = 15) plt.show() |
|
1 |
box_density(boston.iloc[:, -1]) |
Посмотрим на корреляцию количественных признаков с целевой переменной.
|
1 |
boston.drop(columns = 'CHAS').corr().MEDV.to_frame().style.background_gradient() |
Используем точечно-бисериальную корреляцию для оценки взамосвязи переменной CHAS и целевой переменной.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
def pbc(continuous, binary): # преобразуем количественную переменную в массив Numpy continuous_values = np.array(continuous) # классы качественной переменной превратим в нули и единицы binary_values = np.unique(binary, return_inverse = True)[1] # создадим две подгруппы количественных наблюдений # в зависимости от класса дихотомической переменной group0 = continuous_values[np.argwhere(binary_values == 0).flatten()] group1 = continuous_values[np.argwhere(binary_values == 1).flatten()] # найдем средние групп, mean0, mean1 = np.mean(group0), np.mean(group1) # а также длины групп и всего датасета n0, n1, n = len(group0), len(group1), len(continuous_values) # рассчитаем СКО количественной переменной std = continuous_values.std() # подставим значения в формулу return (mean1 - mean0) / std * np.sqrt( (n1 * n0) / (n * (n-1)) ) |
|
1 |
pbc(boston.MEDV, boston.CHAS) |
|
1 |
0.17543361629972762 |
Обработка данных
Пропущенные значения
Посмотрим, есть ли пропущенные значения.
|
1 |
boston.isnull().sum() |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
CRIM 0 ZN 0 INDUS 0 CHAS 0 NOX 0 RM 0 AGE 0 DIS 0 RAD 0 TAX 0 PTRATIO 0 B 0 LSTAT 0 MEDV 0 dtype: int64 |
Выбросы
Удалим выбросы.
|
1 2 3 4 5 6 7 8 9 |
from sklearn.ensemble import IsolationForest clf = IsolationForest(max_samples = 100, random_state = 42) clf.fit(boston) boston['anomaly'] = clf.predict(boston) boston = boston[boston.anomaly == 1] boston = boston.drop(columns = 'anomaly') boston.shape |
|
1 |
(402, 14) |
При удалении выбросов важно помнить, что полное отсутствие вариантивности в данных не позволит выявить взаимосвязи.
Масштабирование признаков
Приведем признаки к одному масштабу (целевую переменную трогать не будем).
|
1 |
boston.iloc[:, :-1] = (boston.iloc[:, :-1] - boston.iloc[:, :-1].mean()) / boston.iloc[:, :-1].std() |
Замечу, что метод наименьших квадратов не требует масштабирования признаков, градиентному спуску же напротив необходимо, чтобы все значения находились в одном диапазоне (подробнее в дополнительных материалах).
Кодирование категориальных переменных
После стандартизации переменная CHAS также сохранила только два значения.
|
1 |
boston.CHAS.value_counts() |
|
1 2 3 |
-0.182581 389 5.463391 13 Name: CHAS, dtype: int64 |
Ее можно не трогать.
Построение модели
Создадим первую пробную (baseline) модель с помощью библиотеки sklearn.
baseline-модель
|
1 2 3 4 5 6 7 |
X = boston.drop('MEDV', axis = 1) y = boston['MEDV'] from sklearn.linear_model import LinearRegression model = LinearRegression() y_pred = model.fit(X, y).predict(X) |
Оценка качества
Диагностика модели, метрики качества и функции потерь
Вероятно, вы заметили, что мы использовали MSE и для обучения модели, и для оценки ее качества. Возникает вопрос, есть ли отличие между функцией потерь и метрикой качества модели.
Функция потерь и метрика качества могут совпадать, а могут и не совпадать. Важно понимать, что у них разное назначение.
- Функция потерь — это часть алгоритма, нам важно, чтобы эта функция была дифференцируема (у нее была производная)
- Производная метрики качества нас не интересует. Метрика качества должна быть адекватна решаемой задаче.
MSE, RMSE, MAE, MAPE
MSE и RMSE
Для оценки качества RMSE предпочтительнее, чем MSE, потому что показывает насколько ошибается модель в тех же единицах измерения, что и целевая переменная. Например, если диапазон целевой переменной от 80 до 100, а RMSE 20, то в среднем модель ошибается на 20-25 процентов.
В качестве практики напишем собственную функцию.
|
1 2 3 4 5 6 7 8 9 |
# параметр squared = True возвращает MSE # параметр squared = False возвращает RMSE def mse(y, y_pred, squared = True): mse = ((y - y_pred) ** 2).sum() / len(y) if squared == True: return mse else: return np.sqrt(mse) |
|
1 |
mse(y, y_pred), mse(y, y_pred, squared = False) |
|
1 |
(9.980044349414223, 3.1591208190593507) |
Сравним с sklearn.
|
1 2 3 4 |
from sklearn.metrics import mean_squared_error # squared = False дает RMSE mean_squared_error(y, y_pred, squared = False) |
|
1 |
3.1591208190593507 |
MAE
Приведем формулу.
$$ MAE = \frac{\sum |y-\hat{y}|}{n} $$
Средняя абсолютная ошибка представляет собой среднее арифметическое абсолютной ошибки $\varepsilon = |y-\hat{y}| $ и использует те же единицы измерения, что и целевая переменная.
|
1 2 |
def mae(y, y_pred): return np.abs(y - y_pred).sum() / len(y) |
|
1 |
mae(y, y_pred) |
|
1 |
2.355651795233418 |
|
1 2 3 |
from sklearn.metrics import mean_absolute_error mean_absolute_error(y, y_pred) |
|
1 |
2.355651795233418 |
MAE часто используется при оценке качества моделей временных рядов.
MAPE
Средняя абсолютная ошибка в процентах (mean absolute percentage error) по сути выражает MAE в процентах, а не в абсолютных величинах, представляя отклонение как долю от истинных ответов.
$$ MAPE = \frac{1}{n} \sum \vert \frac{y-\hat{y}}{y} \vert $$
Это позволяет сравнивать модели с разными единицами измерения между собой.
|
1 2 |
def mape(y, y_pred): return 1/len(y) * np.abs((y - y_pred) / y).sum() |
|
1 |
mape(y, y_pred) |
|
1 |
0.11819866589612749 |
|
1 2 3 |
from sklearn.metrics import mean_absolute_percentage_error mean_absolute_percentage_error(y, y_pred) |
|
1 |
0.11819866589612749 |
Коэффициент детерминации
В рамках вводного курса в ответах на вопросы к занятию по регрессии мы подробно рассмотрели коэффициент детерминации ($R^2$), его связь с RMSE, а также зачем нужен скорректированный $R^2$. Как мы знаем, если использовать, например, класс LinearRegression, то эта метрика содержится в методе .score().
|
1 |
model.score(X, y) |
|
1 |
0.7965234359550825 |
Также можно использовать функцию r2_score() модуля metrics.
|
1 2 3 |
from sklearn.metrics import r2_score r2_score(y, y_pred) |
|
1 |
0.7965234359550825 |
Для скорректированного $R^2$ напишем собственную функцию.
|
1 2 3 4 5 6 7 |
def r_squared(x, y, y_pred): r2 = 1 - ((y - y_pred)** 2).sum()/((y - y.mean()) ** 2).sum() n, k = x.shape r2_adj = 1 - ((y - y_pred)** 2).sum()/((y - y.mean()) ** 2).sum() return r2, r2_adj |
|
1 |
r_squared(X, y, y_pred) |
|
1 |
(0.7965234359550825, 0.7965234359550825) |
Диагностика модели
Теперь проведем диагностику модели в соответствии со сделанными ранее допущениями.
Анализ остатков и прогнозных значений
Напишем диагностическую функцию, которая сразу выведет несколько интересующих нас графиков и метрик, касающихся остатков и прогнозных значений.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
from scipy.stats import probplot from statsmodels.graphics.tsaplots import plot_acf from statsmodels.stats.stattools import durbin_watson def diagnostics(y, y_pred): residuals = y - y_pred residuals_mean = np.round(np.mean(y - y_pred), 3) f, ((ax_rkde, ax_prob), (ax_ry, ax_auto), (ax_yy, ax_ykde)) = plt.subplots(nrows = 3, ncols = 2, figsize = (12, 18)) # в первом подграфике построим график плотности остатков sns.kdeplot(residuals, fill = True, ax = ax_rkde) ax_rkde.set_title('Residuals distribution', fontsize = 14) ax_rkde.set(xlabel = f'Residuals, mean: {residuals_mean}') ax_rkde.set(ylabel = 'Density') # во втором - график нормальной вероятности остатков probplot(residuals, dist = 'norm', plot = ax_prob) ax_prob.set_title('Residuals probability plot', fontsize = 14) # в третьем - график остатков относительно прогноза ax_ry.scatter(y_pred, residuals) ax_ry.set_title('Predicted vs. Residuals', fontsize = 14) ax_ry.set(xlabel = 'y_pred') ax_ry.set(ylabel = 'Residuals') # в четвертом - автокорреляцию остатков plot_acf(residuals, lags = 30, ax = ax_auto) ax_auto.set_title('Residuals Autocorrelation', fontsize = 14) ax_auto.set(xlabel = f'Lags \ndurbin_watson: {durbin_watson(residuals).round(2)}') ax_auto.set(ylabel = 'Autocorrelation') # в пятом - сравним прогнозные и фактические значения ax_yy.scatter(y, y_pred) ax_yy.plot([y.min(), y.max()], [y.min(), y.max()], "k--", lw = 1) ax_yy.set_title('Actual vs. Predicted', fontsize = 14) ax_yy.set(xlabel = 'y_true') ax_yy.set(ylabel = 'y_pred') # в шестом - сравним распределение истинной и прогнозной целевой переменных sns.kdeplot(y, fill = True, ax = ax_ykde, label = 'y_true') sns.kdeplot(y_pred, fill = True, ax = ax_ykde, label = 'y_pred') ax_ykde.set_title('Actual vs. Predicted Distribution', fontsize = 14) ax_ykde.set(xlabel = 'y_true and y_pred') ax_ykde.set(ylabel = 'Density') ax_ykde.legend(loc = 'upper right', prop = {'size': 12}) plt.tight_layout() plt.show() |
|
1 |
diagnostics(y, y_pred) |
Разберем полученную информацию.
- В целом остатки модели распределены нормально с нулевым средним значением
- Явной гетероскедастичности нет, хотя мы видим, что дисперсия не всегда равномерна
- Присутствует умеренная отрицательная корреляция
- График y_true vs. y_pred показывает насколько сильно прогнозные значения отклоняются от фактических. В идеальной модели (без шума, т.е. без случайных колебаний) точки должны были би стремиться находиться на диагонали, в более реалистичной модели нам бы хотелось видеть, что точки плотно сосредоточены вокруг диагонали.
- Распределение прогнозных значений в целом повторяет распределение фактических.
Мультиколлинеарность
Отдельно проведем анализ на мультиколлинеарность. Напишем соответствующую функцию.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
def vif(df, features): vif, tolerance = {}, {} # пройдемся по интересующим нас признакам for feature in features: # составим список остальных признаков, которые будем использовать # для построения регрессии X = [f for f in features if f != feature] # поместим текущие признаки и таргет в X и y X, y = df[X], df[feature] # найдем коэффициент детерминации r2 = LinearRegression().fit(X, y).score(X, y) # посчитаем tolerance tolerance[feature] = 1 - r2 # найдем VIF vif[feature] = 1 / (tolerance[feature]) # выведем результат в виде датафрейма return pd.DataFrame({'VIF': vif, 'Tolerance': tolerance}) |
|
1 |
vif(df = X.drop('CHAS', axis = 1), features = X.drop('CHAS', axis = 1).columns) |
Дополнительная обработка данных
Попробуем дополнительно улучшить некоторые из диагностических показателей.
VIF
Уберем признак с наибольшим VIF (RAD) и посмотрим, что получится.
|
1 2 3 4 5 6 7 8 9 10 11 12 |
vif(df = X, features = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'TAX', 'PTRATIO', 'B', 'LSTAT']) |
Показатели пришли в норму. Окончательно удалим RAD.
|
1 |
boston.drop(columns = 'RAD', inplace = True) |
Преобразование данных
Применим преобразование Йео-Джонсона.
|
1 2 3 4 5 |
from sklearn.preprocessing import PowerTransformer pt = PowerTransformer() boston = pd.DataFrame(pt.fit_transform(boston), columns = boston.columns) |
Отбор признаков
Посмотрим на линейную корреляцию Пирсона количественных признаков и целевой переменной.
|
1 |
boston_t.drop(columns = 'CHAS').corr().MEDV.to_frame().style.background_gradient() |
Также рассчитаем точечно-бисериальную корреляцию.
|
1 |
pbc(boston_t.MEDV, boston_t.CHAS) |
|
1 |
0.11362001957976946 |
Удалим признаки с наименьшей корреляцией, а именно ZN, CHAS, DIS и B.
|
1 |
boston.drop(columns = ['ZN', 'CHAS', 'DIS', 'B'], inplace = True) |
Повторное моделирование и диагностика
Повторное моделирование
Выполним повторное моделирование.
|
1 2 3 4 5 6 7 |
X = boston_t.drop(columns = ['ZN', 'CHAS', 'DIS', 'B', 'MEDV']) y = boston_t.MEDV from sklearn.linear_model import LinearRegression model = LinearRegression() y_pred = model.fit(X, y).predict(X) |
Оценка качества и диагностика
Оценим качество. Так как мы преобразовали целевую переменную, показатель RMSE не будет репрезентативен. Воспользуемся MAPE и $R^2$.
|
1 |
mape(y, y_pred) |
|
1 |
9.427283256967268 |
|
1 |
r_squared(X, y, y_pred) |
|
1 |
(0.7546883769637166, 0.7546883769637166) |
Отклонение прогнозного значения от истинного снизилось. $R^2$ немного уменьшился, чтобы бывает, когда мы пытаемся привести модель к соответствию допущениям. Проведем диагностику.
|
1 |
diagnostics(y, y_pred) |
Распределение остатков немного улучшилось, при этом незначительно усилилась их отрицательная автокорреляция. Распределение целевой переменной стало менее островершинным.
Данные можно было бы продолжить анализировать и улучшать, однако в рамках текущего занятия перейдем к механике обучения модели.
Коэффициенты
Выведем коэффициенты для того, чтобы сравнивать их с результатами построенных с нуля моделей.
|
1 |
model.intercept_, model.coef_ |
|
1 2 3 |
(9.574055157844797e-16, array([-0.09989392, 0.03965441, 0.1069877 , 0.23172172, -0.05561128, -0.16878987, -0.18057055, -0.49319274])) |
Обучение модели
Теперь когда мы поближе познакомились с понятием регрессии, разобрали функции потерь и изучили допущения, при которых модель может быть удачной аппроксимацией данных, пора перейти к непосредственному созданию алгоритмов.
Векторизация уравнения
Для удобства векторизуем приведенное выше уравнение множественной линейной регрессии
$$ y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} X = \begin{bmatrix} x_0 & x_1 & \ldots & x_j \\ x_0 & x_1 & \ldots & x_j \\ \vdots & \vdots & \vdots & \vdots \\ x_{0} & x_{1} & \ldots & x_{n,j} \end{bmatrix}, \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix}, \varepsilon = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} $$
где n — количество наблюдений, а j — количество признаков.
Обратите внимание, что мы создали еще один столбец данных $ x_0 $, который будем умножать на сдвиг $ \theta_0 $. Его мы заполним единицами.
В результате такого несложного преобразования значение сдвига не изменится, но мы сможем записать уравнение через умножение матрицы на вектор.
$$ y = X\theta + \varepsilon $$
Кроме того, как мы увидим ниже, так нам не придется искать отдельную производную для коэффициента $ \theta_0 $.
Схематично для модели с четырьмя наблюдениями (n = 4) и двумя признаками (j = 2) получаются вот такие матрицы.
Функция потерь
Как мы уже говорили, чтобы подобрать оптимальные коэффициенты $\theta$, нам нужен критерий или функция потерь. Логично измерять отклонение прогнозного значения от истинного.
$$ \varepsilon = X\theta-y $$
При этом опять же просто складывать отклонения или ошибки мы не можем. Положительные и отрицательные значения будут взаимоудаляться. Для решения этой проблемы можно, например, использовать модуль и это приводит нас к абсолютной ошибке или L1 loss.
Абсолютная ошибка, L1 loss
При усреднении на количество наблюдений мы получаем среднюю абсолютную ошибку (mean absolute error, MAE).
$$ MAE = \frac{\sum{|y-X\theta|}}{n} = \frac{\sum{|\varepsilon|}}{n} $$
Приведем пример такой функции на Питоне.
|
1 2 |
def L1(y_true, y_pred): return np.sum(np.abs(y_true - y_pred)) / y_true.size |
Помимо модуля ошибку можно возводить в квадрат.
Квадрат ошибки, L2 loss
В этом случай говорят про сумму квадратов ошибок (sum of squared errors, SSE) или сумму квадратов остатков (sum of squared residuals, SSR или residual sum of squares, RSS).
$$ SSE = \sum (y-X\theta)^2 $$
Как мы уже знаем, на практике вместо SSE часто используется MSE, или вернее half MSE для удобства нахождения производной.
$$ MSE = \frac{1}{2n} \sum (y-\theta X)^2 $$
Ниже код на Питоне.
|
1 2 |
def L2(y_true, y_pred): return np.sum((y_true - y_pred) ** 2) / y_true.size |
На практике у обеих функций есть сильные и слабые стороны. Рассмотрим L1 loss (MAE) и L2 loss (MSE) на графике.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
# для построения графиков мы используем x вместо y_true, y_pred # в качестве входящего значения def mse(x): return x ** 2 def mae(x): return np.abs(x) plt.figure(figsize = (10, 8)) x_vals = np.arange(-3, 3, 0.01) plt.plot(x_vals, mae(x_vals), label = 'MAE') plt.plot(x_vals, mse(x_vals), label = 'MSE') plt.legend(loc = 'upper center', prop = {'size': 14}) plt.grid() plt.show() |
Очевидно, при отклонении от точки минимума из-за возведения в квадрат L2 значительно быстрее увеличивает ошибку, поэтому если в данных есть выбросы при суммированнии они очень сильно влияют на этот показатель, хотя де-факто большая часть значений такого уровня потерь не дали бы.
Функция L1 не дает такой большой ошибки на выбросах, однако ее сложно дифференцировать, в точке минимума ее производная не определена.
Функция Хьюбера
Рассмотрим функцию Хьюбера (Huber loss), которая объединяет сильные стороны вышеупомянутых функций и при этом лишена их недостатков. Посмотрим на формулу.
$$ L_{\delta}= \left\{\begin{matrix} \frac{1}{2}(y-\hat{y})^{2} & if | y-\hat{y} | < \delta\\ \delta (|y-\hat{y}|-\frac1 2 \delta) & otherwise \end{matrix}\right. $$
Представим ее на графике.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
plt.figure(figsize = (10, 8)) def huber(x, delta = 1.): huber_mse = 0.5 * np.square(x) huber_mae = delta * (np.abs(x) - 0.5 * delta) return np.where(np.abs(x) <= delta, huber_mse, huber_mae) x_vals = np.arange(-3, 3, 0.01) plt.plot(x_vals, mae(x_vals), label = 'MAE') plt.plot(x_vals, mse(x_vals), label = 'MSE') plt.plot(x_vals, huber(x_vals, delta = 2), label = 'Huber') plt.legend(loc = 'upper center', prop = {'size': 14}) plt.grid() plt.show() |
Также приведем код этой функции.
|
1 2 3 4 5 6 |
def huber(y_pred, y_true, delta = 1.0): # пропишем обе части функции потерь huber_mse = 0.5 * (y_true - y_pred) ** 2 huber_mae = delta * (np.abs(y_true - y_pred) - 0.5 * delta) # выберем одну из них в зависимости от дельта return np.where(np.abs(y_true - y_pred) <= delta, huber_mse, huber_mae) |
На сегодняшнем занятии мы, как и раньше, в качестве функции потерь используем MSE.
Метод наименьших квадратов
Нормальные уравнения
Для множественной линейной регрессии коэффициенты находятся по следующей формуле
$$ \theta = (X^TX)^{-1}X^Ty $$
Разберем эту формулу подробнее. Сумму квадратов остатков (SSE) можно переписать как произведение вектора $ \hat{\varepsilon} $ на самого себя, то есть $ SSE = \varepsilon^{T} \varepsilon$. Помня, что $\varepsilon = y-X\theta $ получаем (не забывая транспонировать)
$$ (y-X\theta)^T(y-X\theta) $$
Раскрываем скобки
$$ y^Ty-y^T(X\theta)-(X\theta)^Ty+(X\theta)^T(X\theta) $$
Заметим, что $A^TB = B^TA$, тогда
$$ y^Ty-(X\theta)^Ty-(X\theta)^Ty+(X\theta)^T(X\theta)$$
$$ y^Ty-2(X\theta)^Ty+(X\theta)^T(X\theta) $$
Вспомним, что $(AB)^T = A^TB^T$, тогда
$$ y^Ty-2\theta^TX^Ty+\theta^TX^TX\theta $$
Теперь нужно найти градиент этой функции
$$ \nabla_{J(\theta)} \left( y^Ty-2\theta^TX^Ty+\theta^TX^TX\theta \right) $$
После дифференцирования относительно $\theta$ мы получаем следующий результат
$$ \nabla_{J(\theta)} = -2X^Ty+2X^TX\theta $$
Как мы помним, оптимум функции находится там, где производная равна нулю.
$$ -2X^Ty+2X^TX\theta = 0 $$
$$ -X^Ty+X^TX\theta = 0 $$
$$ X^TX\theta = X^Ty $$
Выражение выше называется нормальным уравнением (normal equation). Решив его для $\theta$ мы найдем аналитическое решение минимизации суммы квадратов отклонений.
$$ \theta = (X^TX)^{-1}X^Ty $$
По теореме Гаусса-Маркова, оценка через МНК является наиболее оптимальной (обладающей наименьшей дисперсией) среди всех методов построения модели.
Код на Питоне
Перейдем к созданию класса линейной регрессии наподобие LinearRegression библиотеки sklearn. Вначале напишем функцию гипотезы (т.е. функцию самой модели), снабдив ее еще одной функцией, которая добавляет столбец из единиц к признакам.
$$ h_{\theta}(x) = \theta X $$
|
1 2 3 |
def add_ones(x): # важно! метод .insert() изменяет исходный датафрейм return x.insert(0,'x0', np.ones(x.shape[0])) |
|
1 2 3 4 |
def h(x, theta): x = x.copy() add_ones(x) return np.dot(x, theta) |
Перейдем к функции, отвечающей за обучение модели.
$$ \theta = (X^TX)^{-1}X^Ty $$
|
1 2 3 4 5 6 7 8 9 10 11 |
# строчную `x` используем внутри функций и методов класса # заглавную `X` вне функций и методов def fit(x, y): x = x.copy() add_ones(x) xT = x.transpose() inversed = np.linalg.inv(np.dot(xT, x)) thetas = inversed.dot(xT).dot(y) return thetas |
Обучим модель и выведем коэффициенты.
|
1 2 |
thetas = fit(X, y) thetas[0], thetas[1:] |
|
1 2 3 |
(9.3718435789647e-16, array([-0.09989392, 0.03965441, 0.1069877 , 0.23172172, -0.05561128, -0.16878987, -0.18057055, -0.49319274])) |
Примечание. Замечу, что не все матрицы обратимы. В этом случае можно найти псевдообратную матрицу (pseudoinverse). Для этого в Numpy есть функция np.linalg.pinv().
Сделаем прогноз.
|
1 2 |
y_pred = h(X, thetas) y_pred[:5] |
|
1 |
array([1.24414666, 0.55999778, 1.48103299, 1.49481605, 1.21342788]) |
Создание класса
Объединим код в класс.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
class ols(): def __init__(self): self.thetas = None def add_ones(self, x): return x.insert(0,'x0', np.ones(x.shape[0])) def fit(self, x, y): x = x.copy() self.add_ones(x) xT = x.T inversed = np.linalg.inv(np.dot(xT, x)) self.thetas = inversed.dot(xT).dot(y) def predict(self, x): x = x.copy() self.add_ones(x) return np.dot(x, self.thetas) |
Создадим объект класса и обучим модель.
|
1 2 3 |
model = ols() model.fit(X, y) |
Выведем коэффициенты.
|
1 |
model.thetas[0], model.thetas[1:] |
|
1 2 3 |
(9.3718435789647e-16, array([-0.09989392, 0.03965441, 0.1069877 , 0.23172172, -0.05561128, -0.16878987, -0.18057055, -0.49319274])) |
Сделаем прогноз.
|
1 2 |
y_pred = model.predict(X) y_pred[:5] |
|
1 |
array([1.24414666, 0.55999778, 1.48103299, 1.49481605, 1.21342788]) |
Оценка качества
Оценим качество через MAPE и $R^2$.
|
1 |
mape(y, y_pred) |
|
1 |
9.427283256967264 |
|
1 |
r_squared(X, y, y_pred) |
|
1 |
(0.7546883769637167, 0.7546883769637167) |
Мы видим, что результаты аналогичны.
Метод градиентного спуска
В целом с этим методом мы уже хорошо знакомы. В качестве упражнения давайте реализуем этот алгоритм на Питоне для многомерных данных.
Нахождение градиента
Покажем расчет градиента на схеме.
В данном случае мы берем датасет из четырех наблюдений и двух признаков ($x_1$ и $x_2$) и соответственно используем три коэффициента ($\theta_0, \theta_1, \theta_2$).
Пошаговое построение модели
Начнем с функции гипотезы.
$$ h_{\theta}(x) = \theta X $$
|
1 2 |
def h(x, thetas): return np.dot(x, thetas) |
Объявим функцию потерь.
$$ J({\theta_j}) = \frac{1}{2n} \sum (y-\theta X)^2 $$
|
1 2 |
def objective(x, y, thetas, n): return np.sum((y - h(x, thetas)) ** 2) / (2 * n) |
Объявим функцию для расчета градиента.
$$ \frac{\partial}{\partial \theta_j} J(\theta) = -x_j(y — X\theta) \times \frac{1}{n} $$
где j — индекс признака.
|
1 2 |
def gradient(x, y, thetas, n): return np.dot(-x.T, (y - h(x, thetas))) / n |
Напишем функцию для обучения модели.
$$ \theta_j := \theta_j-\alpha \frac{\partial}{\partial \theta_j} J(\theta) $$
Символ := означает, что левая часть равенства определяется правой. По сути, с каждой итерацией мы обновляем веса, умножая коэффициент скорости обучения на градиент.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
def fit(x, y, iter = 20000, learning_rate = 0.05): x, y = x.copy(), y.copy() # функцию add_ones() мы объявили выше add_ones(x) thetas, n = np.zeros(x.shape[1]), x.shape[0] loss_history = [] for i in range(iter): loss_history.append(objective(x, y, thetas, n)) grad = gradient(x, y, thetas, n) thetas -= learning_rate * grad return thetas, loss_history |
Обучим модель, выведем коэффициенты и достигнутый (минимальный) уровень ошибки.
|
1 |
thetas, loss_history = fit(X, y, iter = 50000, learning_rate = 0.05) |
|
1 |
thetas[0], thetas[1:], loss_history[-1] |
|
1 2 3 4 |
(9.493787734953824e-16, array([-0.09989392, 0.03965441, 0.1069877 , 0.23172172, -0.05561128, -0.16878987, -0.18057055, -0.49319274]), 0.1226558115181417) |
Полученный результат очень близок к тому, что было найдено методом наименьших квадратов.
Прогноз
Сделаем прогноз.
|
1 2 3 4 |
def predict(x, thetas): x = x.copy() add_ones(x) return np.dot(x, thetas) |
|
1 2 |
y_pred = predict(X, thetas) y_pred[:5] |
|
1 |
array([1.24414666, 0.55999778, 1.48103299, 1.49481605, 1.21342788]) |
Создание класса
Объединим написанные функции в класс.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
class gd(): def __init__(self): self.thetas = None self.loss_history = [] def add_ones(self, x): return x.insert(0,'x0', np.ones(x.shape[0])) def objective(self, x, y, thetas, n): return np.sum((y - self.h(x, thetas)) ** 2) / (2 * n) def h(self, x, thetas): return np.dot(x, thetas) def gradient(self, x, y, thetas, n): return np.dot(-x.T, (y - self.h(x, thetas))) / n def fit(self, x, y, iter = 20000, learning_rate = 0.05): x, y = x.copy(), y.copy() self.add_ones(x) thetas, n = np.zeros(x.shape[1]), x.shape[0] # объявляем переменную loss_history (отличается от self.loss_history) loss_history = [] for i in range(iter): loss_history.append(self.objective(x, y, thetas, n)) grad = self.gradient(x, y, thetas, n) thetas -= learning_rate * grad # записываем обратно во внутренние атрибуты, чтобы передать методу .predict() self.thetas = thetas self.loss_history = loss_history def predict(self, x): x = x.copy() self.add_ones(x) return np.dot(x, self.thetas) |
Создадим объект класса, обучим модель, выведем коэффициенты и сделаем прогноз.
|
1 2 3 4 5 |
model = gd() model.fit(X, y, iter = 50000, learning_rate = 0.05) model.thetas[0], model.thetas[1:], model.loss_history[-1] |
|
1 2 3 4 |
(9.493787734953824e-16, array([-0.09989392, 0.03965441, 0.1069877 , 0.23172172, -0.05561128, -0.16878987, -0.18057055, -0.49319274]), 0.1226558115181417) |
|
1 2 |
y_pred = model.predict(X) y_pred[:5] |
|
1 |
array([1.24414666, 0.55999778, 1.48103299, 1.49481605, 1.21342788]) |
Оценка качества
|
1 |
mape(y, y_pred) |
|
1 |
9.427283256967264 |
|
1 |
r_squared(X, y, y_pred) |
|
1 |
(0.7546883769637167, 0.7546883769637167) |
Теперь рассмотрим несколько дополнительных соображений, касающихся построения модели линейной регрессии.
Диагностика алгоритма
Работу алгоритма можно проверить с помощью кривой обучения (learning curve).
- Ошибка постоянно снижается
- Алгоритм остановится, после истечения заданного количества итераций
- Можно задать пороговое значение, после которого он остановится (например, $10^{-1}$)
Построим кривую обучения.
|
1 2 |
plt.plot(loss_history) plt.show() |
|
1 2 |
plt.plot(loss_history[:100]) plt.show() |
Она также позволяет выбрать адекватный коэффициент скорости обучения.
Подведем итог
Сегодня мы подробно рассмотрели модель множественной линейной регрессиии. В частности, мы поговорили про построение гипотезы, основные функции потерь, допущения модели линейной регрессии, метрики качества и диагностику модели.
Кроме того, мы узнали как изнутри устроены метод наименьших квадратов и метод градиентного спуска и построили соответствующие модели на Питоне.
Отдельно замечу, что, изучив скорректированный коэффициент детерминации, мы начали постепенно обсуждать способы усовершенствования базовых алгоритмов и метрик. На последующих занятиях мы продолжим этот путь в двух направлениях: познакомимся со способами регуляризации функции потерь и начнем создавать более сложные алгоритмы оптимизации.
Но прежде предлагаю в деталях изучить уже знакомый нам алгоритм логистической регрессии.