Глубокое погружение в множественную линейную регрессию с примерами в Python
В первой части этой статьи мы формально определили задачу линейной регрессии и показали, как решать простые задачи линейной регрессии, когда набор данных содержит только один признак. Во второй части статьи мы обсудим задачи множественной линейной регрессии, где набор данных может содержать любое количество признаков.
Сначала мы обобщим решение в закрытой форме, которое мы нашли для простой линейной регрессии, на любое количество признаков. Затем мы предложим альтернативный подход к решению задач линейной регрессии, основанный на градиентном спуске, и обсудим плюсы и минусы этого подхода по сравнению с использованием решения в закрытой форме. Кроме того, мы изучим классы Scikit-Learn, которые реализуют оба подхода, и продемонстрируем, как их использовать на реальном наборе данных.
Определение множественной линейной регрессии
Напомним, что в задачах регрессии нам дается набор из n помеченных примеров: D = {(x₁, y₁), (x₂, y₂), … , (xₙ, yₙ)}, где xᵢ — m-мерный вектор, содержащий признакииз примера i,и yᵢ – реальное значение, представляющее метку этого примера.
В задачах линейной регрессии мы предполагаем, что существует линейная связь между вектором признаков x и меткой y, поэтому наша модельная гипотеза принимает следующую форму:
Наша цель — найти параметры w этой модели, которые минимизируют сумму квадратов невязок:
В предыдущей части статьи мы показали, как найти оптимальное w для случая m = 1, используя нормальные уравнения. Теперь мы расширим эти уравнения для любого количества функций m.
Чтобы упростить вывод нормальных уравнений для общего случая, мы сначала определяем матрицу X, которая содержит значения всех признаков в наборе данных, включая члены пересечения:
Эта матрица называется матрицей дизайна. Каждая строка в матрице плана представляет отдельный образец, а столбцы представляют собой независимые переменные. Размеры матрицы равны n × (m + 1), где n — количество выборок, а m количество признаков.
Кроме того, мы определяем вектор y как n-мерный вектор, содержащий все целевые значения:
Эти определения позволяют нам записать функцию стоимости обычного метода наименьших квадратов (OLS) в следующей матричной форме:
Доказательство:
Прежде всего отметим, что:
Скалярное произведение вектора с самим собой uᵗu — это просто сумма квадратов всех его компонентов, поэтому мы имеем:
Закрытое решение
Как и в случае простой линейной регрессии, функция J(w) является выпуклой, поэтому она имеет один локальный минимум, который также является глобальным минимумом. Чтобы найти этот глобальный минимум, мы вычисляем градиент J(w) относительно w и устанавливаем его равным нулю.
Градиент J(w) относительно w равен:
Во время этого доказательства мы использовали некоторые основные правила матричного исчисления, которые объясняются в этой статье:
Теперь мы устанавливаем этот градиент равным нулю, чтобы получить нормальные уравнения:
Таким образом, оптимальный w*, который минимизирует функцию стоимости метода наименьших квадратов, равен:
Обратите внимание, что здесь мы предполагали, что столбцы X линейно независимы (т. е. X имеет полный ранг столбца), в противном случае XᵗX необратима, и нет единственного решения для w*.
Когда столбцы X линейно зависимы, мы называем это явление мультиколлинеарностью. Математически набор переменных является совершенно мультиколлинеарным, если для всех выборок i:
где λₖ — константы, а xᵢₖ — значение признака k в образце i.
На практике полная мультиколлинеарность встречается редко (например, она может быть вызвана случайным дублированием переменной в данных). Однако даже меньшая степень мультиколлинеарности, когда два или более признаков сильно коррелируют друг с другом (но не полностью коррелированы), может вызвать проблемы как при подборе модели (коэффициенты становятся очень чувствительными к небольшим изменениям в данных), так и при интерпретации. результаты (становится трудно определить, какие функции оказывают наибольшее влияние на прогнозы модели).
Заинтересованный читатель может найти больше информации о проблеме мультиколлинеарности и способах ее решения в этой статье в Википедии.
Пример множественной линейной регрессии
Чтобы продемонстрировать использование закрытого решения, давайте построим модель линейной регрессии для набора данных о жилье в Калифорнии, доступного в модуле sklearn.datasets. Целью этого набора данных является прогнозирование средней стоимости дома в данном районе (домовом квартале) в Калифорнии на основе 8 различных характеристик этого района, таких как средний доход или среднее количество комнат на домохозяйство.
Сначала мы импортируем необходимые библиотеки и инициализируем случайное начальное число, чтобы получить воспроизводимые результаты:
import numpy as np import pandas as pd import matplotlib.pyplot as plt np.random.seed(0)
Затем мы получаем набор данных:
from sklearn.datasets import fetch_california_housing data = fetch_california_housing() X, y = data.data, data.target feature_names = data.feature_names
Чтобы изучить набор данных, мы объединяем функции (X) и метки (y) в Pandas DataFrame и отображаем первые строки из таблицы:
mat = np.column_stack((X, y)) df = pd.DataFrame(mat, columns=np.append(feature_names, 'MedianValue')) df.head()
Мы можем дополнительно исследовать набор данных, вызвав метод info() DataFrame, который предоставляет информацию о типе столбцов и о том, содержат ли они какие-либо пропущенные значения:
df.info()
К счастью, этот набор данных содержит только числовые признаки и не содержит пропущенных значений. Поэтому здесь не требуется предварительная обработка данных (решение в закрытой форме не требует нормализации данных).
Теперь мы разделим данные на 80% обучающий набор и 20% тестовый набор:
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)
Давайте теперь напишем общую функцию для поиска оптимального w* для любого заданного набора данных, используя решение в закрытой форме, которое мы нашли ранее:
def closed_form_solution(X, y):
w = np.linalg.inv(X.T @ X) @ X.T @ y
return w
Мы можем реализовать закрытое решение в одной строке кода!
Прежде чем мы сможем использовать эту функцию для поиска оптимального w, нам нужно добавить столбец из 1 к матрице X_train, чтобы представить термин перехвата. Это легко сделать с помощью функции np.column_stack():
X_train_b = np.column_stack((np.ones(len(X_train)), X_train))
Теперь мы готовы подогнать нашу модель к тренировочному набору:
w = closed_form_solution(X_train_b, y_train) w
Оптимальный w:
array([-3.68585691e+01, 4.33333407e-01, 9.29324337e-03, -9.86433739e-02,
5.93215487e-01, -7.56192502e-06, -4.74516383e-03, -4.21449336e-01,
-4.34166041e-01])
Этот вектор состоит из 9 компонентов: по одному для каждой из 8 функций в наборе данных и дополнительного веса для точки пересечения (w₀).
Давайте теперь оценим модель на обучающем и тестовом наборах. Важно оценить вашу модель на обоих из них, потому что большое расхождение между результатами обучения и тестами может указывать на то, что ваша модель переоснащена.
Начнем с нахождения показателя R² на тренировочном наборе. Для этого мы сначала получаем предсказания модели на обучающих примерах, умножая матрицу X_train_b на вектор w:
y_train_pred = X_train_b @ w
Теперь мы используем функцию r2_score() из sklearn.metrics, чтобы найти оценку R² в обучающем наборе:
from sklearn.metrics import r2_score
train_score = r2_score(y_train, y_train_pred)
print(f'R2 score (train): {train_score:.5f}')
Результат, который мы получаем:
R2 score (train): 0.60890
Сделаем то же самое на тестовом наборе. Нам нужно добавить столбец из 1 к X_test, прежде чем мы получим прогнозы модели:
X_test_b = np.column_stack((np.ones(len(X_test)), X_test))
y_test_pred = X_test_b @ w
test_score = r2_score(y_test, y_test_pred)
print(f'R2 score (test): {test_score:.5f}')
Результат, который мы получаем:
R2 score (test): 0.59432
Оценка невысокая, что указывает на то, что взаимосвязь между признаками и меткой может быть нелинейной. В таких случаях лучшие результаты могут дать модели нелинейной регрессии, такие как деревья регрессии или k-ближайших соседей.
Упражнение
Допустим, мы случайно продублировали каждую точку в наборе данных, а затем снова запустили линейную регрессию. Как это повлияет на веса модели?
Подсказка. Подумайте, что произойдет с матрицей проекта X и вектором меток y, и как эти изменения повлияют на нормальные уравнения.
Решение можно найти в конце этой статьи.
Линейная регрессия в Scikit-Learn
Scikit-Learn предоставляет класс LinearRegression, который также реализует закрытое решение обычной задачи наименьших квадратов.
По умолчанию этот класс автоматически добавляет столбец с единицами в матрицу проектирования, поэтому вам не нужно добавлять его вручную, как мы делали ранее (если только в конструкторе вы не установили для параметра fit_intercept значение False).
Давайте используем этот класс, чтобы подогнать модель регрессии к тому же набору данных:
from sklearn.linear_model import LinearRegression reg = LinearRegression() reg.fit(X_train, y_train)
Подходящий вектор параметров w хранится в двух атрибутах этого класса:
- coef_ — это массив, содержащий все веса, кроме точки пересечения.
- intercept_ – термин перехвата (w₀).
Распечатаем их:
print(reg.intercept_) print(reg.coef_)
Результат:
-36.858569106801234 [ 4.33333407e-01 9.29324337e-03 -9.86433739e-02 5.93215487e-01 -7.56192502e-06 -4.74516383e-03 -4.21449336e-01 -4.34166041e-01]
Мы получаем точно такие же коэффициенты, как и при вычислении в NumPy.
Метод score() этого класса возвращает показатель R² модели. Для этого требуется только матрица X и вектор y набора данных, по которому вы хотите получить оценку (поэтому нет необходимости вычислять прогнозы модели). Например, мы можем получить оценку R² для обучающего и тестового наборов следующим образом:
train_score = reg.score(X_train, y_train)
print(f'R2 score (train): {train_score:.5f}')
test_score = reg.score(X_test, y_test)
print(f'R2 score (test): {test_score:.5f}')
R2 score (train): 0.60890
R2 score (test): 0.59432
Как и ожидалось, мы получаем те же оценки R², что и раньше.
Анализ ошибок регрессии
В дополнение к оценке общей производительности модели мы часто хотим исследовать поведение наших ошибок регрессии. Например, нормально ли распределяются ошибки вокруг 0 или они искажены? Существуют ли какие-либо входные данные, для которых наша модель имеет особенно высокие ошибки прогнозирования?
Ответы на эти вопросы помогут нам найти источник этих ошибок. Например, если ошибки обычно не распределяются вокруг 0, это может указывать на то, что модель линейной регрессии не подходит для нашего набора данных, и нам нужно попробовать другие модели регрессии (например, полиномиальную регрессию). Или, если наша модель имеет особенно высокие ошибки предсказания для некоторых выборок, они могут быть выбросами, и нам нужно выяснить, откуда они берутся.
График, который поможет вам ответить на эти вопросы, называется графиком остатков. На этом графике показаны невязки по оси y и предсказанные значения модели по оси x.
Давайте напишем функцию для создания этого графика:
def plot_residuals(y_train_pred, y_train, y_test_pred, y_test):
plt.scatter(y_train_pred, y_train_pred - y_train, s=2, marker='o', c='b', label='Training')
plt.scatter(y_test_pred, y_test_pred - y_test, s=2, marker='s', c='m', label='Test')
xmin = min(y_train_pred.min(), y_test_pred.min())
xmax = max(y_train_pred.max(), y_test_pred.max())
plt.hlines(y=0, xmin=xmin, xmax=xmax, color='black')
plt.xlim(xmin, xmax)
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend()
Теперь мы можем вызвать эту функцию, чтобы показать остатки как на обучающем, так и на тестовом наборах:
plot_residuals(y_train_pred, y_train, y_test_pred, y_test)
Мы видим, что большинство ошибок симметрично распределены вокруг 0, но есть некоторые выбросы на дальних концах входного диапазона, которые могут потребовать дальнейшего изучения.
Упражнение
Загрузите набор данных Students Marks dataset с Kaggle. Постройте модель линейной регрессии, чтобы предсказать оценку учащегося на основе времени обучения и количества пройденных курсов. Вычислите среднеквадратичную ошибку и показатель R² модели как в обучающем, так и в тестовом наборах. Постройте остатки по сравнению с предсказанными значениями. Что вы можете узнать из этого графика в наборе данных?
Градиентный спуск
Хотя решение в закрытой форме дает нам прямой способ найти оптимальные параметры регрессионной модели, оно имеет несколько недостатков:
- Решение в закрытой форме неэффективно с точки зрения вычислений, когда у нас есть большое количество функций, поскольку оно требует вычисления обратной величины XᵗX, которая представляет собой m × m матрица (m — количество признаков). Вычисление обратной матрицы имеет сложность времени выполнения O(m³) в большинстве реализаций.
- Для этого необходимо иметь в памяти всю матрицу проекта X, что не всегда возможно, если у нас очень большой набор данных.
- Он не поддерживает интерактивное (пошаговое) обучение, так как любое изменение в матрице проекта X требует повторного вычисления обратной величины XᵗX.
К счастью, есть альтернативный подход к нахождению оптимального w — градиентный спуск. Градиентный спуск — это итеративный подход к поиску минимума функции, при котором мы делаем небольшие шаги в направлении, противоположном градиенту, чтобы приблизиться к минимуму:
Чтобы использовать градиентный спуск для нахождения минимума стоимости метода наименьших квадратов, нам нужно вычислить частные производные J(w) по каждому из весов .
Частная производная J(w) по любому из весов wⱼ равна:
Следовательно, правило обновления градиентного спуска:
где α — скорость обучения, которая контролирует размер шага (0 ‹ α‹ 1).
Вместо того, чтобы обновлять каждый компонент w по отдельности, мы можем обновить весь вектор w за один шаг:
Градиентный спуск может применяться в одном из следующих режимов:
- Пакетный градиентный спуск — веса обновляются после того, как мы вычисляем ошибку на всем обучающем наборе.
- Стохастический градиентный спуск (SGD) — шаг градиентного спуска выполняется после каждого обучающего примера. В этом случае правило обновления градиентного спуска принимает следующий вид:
SGD обычно сходится быстрее, чем пакетный градиентный спуск, поскольку он продвигается после каждого примера, а также поддерживает онлайн-обучение, поскольку может обрабатывать новые точки данных по одной за раз. С другой стороны, SGD менее стабилен, чем пакетный градиентный спуск, и его сходимость к глобальному оптимуму не гарантируется (хотя на практике он очень близок к оптимуму).
Обратите внимание: всякий раз, когда вы используете градиентный спуск, вы должны убедиться, что ваш набор данных нормализован (иначе градиентный спуск может выполнять шаги разного размера в разных направлениях, что сделает его нестабильным).
Класс SGDRegressor
Класс SGDRegressor в Scikit-Learn реализует подход SGD для подбора модели линейной регрессии. Конкретный тип регрессионной модели определяется параметром потери. По умолчанию его значение равно ‘squared_error’, что означает соответствие обычной модели наименьших квадратов. Другие параметры включают huber и epsilon_insensitive (функция потерь, используемая в регрессии опорных векторов), которые будут обсуждаться в следующих статьях.
Помимо функции потерь конструктор этого класса имеет еще несколько важных параметров:
- штраф — используемый тип регуляризации (по умолчанию «l2»).
- alpha — коэффициент регуляризации (по умолчанию 0,0001).
- max_iter — максимальное количество эпох в обучающих данных (по умолчанию 1000).
- learning_rate — график скорости обучения для обновлений веса (по умолчанию «invscaling»).
- eta0 —используемая начальная скорость обучения (по умолчанию 0,01).
- early_stopping — следует ли останавливать обучение, когда оценка проверки не улучшается (по умолчанию False).
- validation_fraction — доля обучающего набора, отведенная для проверки (по умолчанию 0,1).
Поскольку нам нужно нормализовать наши данные перед использованием SGD, мы создадим конвейер, состоящий из двух шагов:
- StandardScaler, который нормализует функции, удаляя среднее значение и масштабируя их до единичной дисперсии.
- SGDRegressor с настройками по умолчанию.
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import SGDRegressor
from sklearn.pipeline import Pipeline
pipeline = Pipeline([
('scaler', StandardScaler()),
('reg', SGDRegressor())
])
Давайте подгоним конвейер к нашему обучающему набору:
pipeline.fit(X_train, y_train)
А теперь оценим модель на обучающей и тестовой выборках:
train_score = pipeline.score(X_train, y_train)
print(f'R2 score (train): {train_score:.5f}')
test_score = pipeline.score(X_test, y_test)
print(f'R2 score (test): {test_score:.5f}')
Оценки, которые мы получаем:
R2 score (train): -485.79460 R2 score (test): -4485.30424
Это очень плохие показатели! Что только что произошло?
Когда вы получаете такие плохие оценки с градиентным спуском, это обычно означает, что ваша скорость обучения слишком высока, что заставляет алгоритм колебаться между двумя сторонами минимума:
Давайте уменьшим скорость обучения с 0,01 до 0,001, изменив параметр eta0 SGDRegressor:
pipeline.set_params(reg__eta0=0.001)
Давайте подгоним конвейер к обучающему набору и переоценим его:
pipeline.fit(X_train, y_train)
train_score = pipeline.score(X_train, y_train)
print(f'R2 score (train): {train_score:.5f}')
test_score = pipeline.score(X_test, y_test)
print(f'R2 score (test): {test_score:.5f}')
На этот раз мы получаем следующие баллы:
R2 score (train): 0.60188 R2 score (test): 0.58393
которые аналогичны показателям R², которые мы получили с решением в закрытой форме (они немного ниже, поскольку SGD приближается к глобальному минимуму, но не к самому минимуму).
Ключевые выводы
- В линейной регрессии мы пытаемся найти линейную связь между набором функций и целевой переменной.
- Ключевые предположения заключаются в том, что признаки линейно независимы, а члены ошибки независимы друг от друга и нормально распределены с нулевым средним значением.
- В обычной регрессии методом наименьших квадратов (OLS) мы пытаемся минимизировать сумму квадратов остатков. Функция стоимости МНК выпукла, поэтому имеет единственную точку минимума.
- Двумя основными подходами для нахождения оптимальных параметров модели являются решение в замкнутой форме и градиентный спуск. Градиентный спуск предпочтительнее, когда набор данных велик или когда нам нужно поддерживать онлайн-обучение.
- При использовании градиентного спуска важно нормализовать ваши функции и выбрать соответствующую скорость обучения.
- Мы оцениваем эффективность регрессионных моделей с помощью показателя R², который варьируется от 0 до 1 и измеряет, насколько модель лучше базовой модели, которая всегда предсказывает среднее значение целевого значения, и RMSE, который квадратный корень из среднеквадратичной ошибки.
Решение упражнения с повторяющимися данными
Напомним, что решение в закрытой форме:
Если мы удвоим точки данных, то вместо X и y у нас будет следующая матрица проектирования и целевой вектор:
Матрица A имеет 2n строк и m столбцов, где строки 1, …, n матрицы идентичны строкам n+ 1, …, 2n. Точно так же вектор z имеет 2n строк, где первые n строки идентичны последним n строкам. .
Мы можем легко показать, что
Доказательство:
Точно так же мы имеем (мы можем рассматривать z как матрицу с одним столбцом):
Следовательно, мы можем написать:
Получаем те же веса, что и в исходной модели! т. е. модель регрессии не изменится.
Заключительные примечания
Все изображения, если не указано иное, принадлежат автору.
Вы можете найти примеры кода этой статьи на моем github: https://github.com/roiyeho/medium/tree/main/multiple_linear_regression
Спасибо за прочтение!
