Множественная регрессия и полиномиальная регрессия

 

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

Напомним, мы начали с одномерной линейной регрессии, где имели входные данные {(x1, y1), (x2, y2), …, (xN, yN)}.

Определение многомерной задачи и ее решение

Сейчас же, поскольку рассматривается сразу несколько факторов, то и наше xi становится вектором – он ещё называется вектором признаков. Количество элементов вектора называется размерностью и, как правило, обозначается буквой D.  Наша модель имеет вид

ŷ=wTx.

Поскольку нам необходимо умножать x на транспонированный вектор параметров w, вектор w также должен иметь размерность D.

Заметьте, что мы всегда можем включить свободный член b в вектор параметров w. Это связано с тем, что мы можем переименовать b в w0 и добавить к нему x0, установив значение x0 = 1:

ŷ = b + w1x1 + … + wDxD,

ŷ = w0 + w1x1 + … + wDxD,

ŷ = w0x0 + w1x1 + … + wDxD = w’Tx’, x0 = 1.

Это эквивалентно добавлению колонки единиц в нашу матрицу данных X (имеющую первоначальную размерность NxD).

Почему наша матрица данных имеет размерность NxD? Не забывайте, N – это количество опытов (наблюдений), а D – количество входных факторов. Если мы возьмём одну строку из X – что означает, что мы взяли результат только одного наблюдения, – мы получим вектор-строку xi размерности 1xD, который и есть вектором признаков. Но поскольку в линейной алгебре векторы принято рассматривать как векторы-столбцы размерностью Dx1, в случае одного опыта нам приходится преобразовать формулу так:

ŷ=wTx.

Если же мы хотим вычислить значение исходящей переменной ŷ для всех наблюдений N одновременно, мы вынуждены записать

 \rightarrow Y_N_X_1 = X_N_X_D W_D_X_1.

Это связано с тем, что вектор-столбец w имеет размерность Dx1, а матрица данных X – размерность NxD, а как вы знаете, для корректного перемножения матриц необходимо, чтобы их соответствующие размерности совпадали. Тогда в качестве результата мы получим правильный вектор-столбец  размерности Nx1. Конечно, сначала это выглядит несколько странно, поскольку наше w  теперь находится справа.

Чтобы прояснить эту идею, давайте рассмотрим простой пример. Пусть имеем N=4, то есть количество проделанных наблюдений равно четырём, а D=3, то есть количество изучаемых входящих факторов равно трём. Таким образом, наша матрица данных X имеет размерность 4×3. Пусть также w – вектор с тремя элементами 1, 2, 3:

X =\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\\1&1&0\end{bmatrix}_{4x3} ; \;w = \begin{bmatrix}1\\2\\3\end{bmatrix}_{3x1}

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

y_i = w^Tx_i = [1\;\; 2 \;\; 3] \; \begin{bmatrix}0\\0\\1\end{bmatrix} = 1*0 + 2*0 + 3*1 = 3 = x^T_i w.

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

Но такие вычисления малоэффективны. Нам бы хотелось максимально использовать возможности NumPy, проделав при этом как можно меньше операций. Не забывайте, NumPy специально оптимизирован для численных расчётов, поэтому нам можно не вычислять y1, y2, y3 и y4 по отдельности. Мы хотим рассчитать их сразу все и чтобы все четыре y были в одном массиве размерностью четыре. Мы можем это сделать, сразу перемножив всю матрицу X на w.

 y_{4x1} = \begin{bmatrix}3\\2\\1\\3\end{bmatrix}_{4x1} = \begin{bmatrix}0&0&1\\0&1&0\\1&0&0\\1&1&0\end{bmatrix}_{4x3} \begin{bmatrix}1\\2\\3\end{bmatrix}_{3x1}

Это может вначале показаться непонятным, но не забывайте, в матрице данных X все наши xi размещаются горизонтально в одной строке, а если мы рассматриваем xi отдельно, он берётся как вектор-столбец, поэтому данные идут по вертикали. Просто помните всегда «золотое правило»: при перемножении матриц их внутренние размерности должны совпадать.

Итак, вернёмся к постановке задачи. Как же нам вычислить w?Тут мы встречаемся с уже знакомыми нам вычислениями. Функция ошибок не изменилась, да и переменные y и ŷ по-прежнему остаются скалярными величинами. Но поскольку задача усложнилась, теперь это выражение имеет несколько иной вид:

E =\sum_{i=1} ^ {N} (y_{i} -\widehat y_i)^2=\sum_{i=1}^{N} (y_{i}-w^{T}x_{i})^2.

Но мы по-прежнему можем взять производные по всем элементам w . Назовём их wj, j=1..D.

\frac{dE}{dw_j} =\sum_{i=1}^{N}2(y_i - w^T x_i)(-\frac{d(w^T x_i)}{d w_j})=\sum_{i=1}^{N}2(y_i - w^T x_i)(- x_i_j).

Если вы сразу не поняли, как получилось такое решение, возможно, вам будет легче, если перепишете уравнение в скалярной форме. В любом случае приостановитесь на минуту, чтобы убедиться, что вы поняли решение. Наше X является матрицей, поэтому её элементы можно пронумеровать двумя индексами xij, где i показывает номер наблюдения, а j – номер фактора (признака).

Как вы помните, следующее действие – приравнять производную к нулю. Поскольку j имеет размерность D и пробегает все значения от 1 до D, то мы имеем D уравнений с D неизвестными.

\frac{dE}{dw_j}=\sum_{i=1}^{N}2(y_i-w^T x_i)(-x_i_j)=0.

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

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

\sum _ {i=1} ^ {N}y_i(- x_i_j) - \sum_{i=1}^ {N} w^Tx_i(- x_i_j) = 0,

 \sum_ {i=1} ^ {N} w^T x_i x_i_j = \sum_ {i=1} ^ {N} y_i x_i_j.

Это уже кое-что, но далеко не всё, так как нам нужно ещё выделить w. Прежде всего мы можем вынести w из суммы, так как суммирование идёт по i, а w не зависит от i.

 w^T \sum_{i=1}^{N} x_i x_i_j = \sum_{i=1}^{N} y_i x_i_j.

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

w^T \sum_{i=1} ^ {N} x_i x_i_1 = \sum_{i=1}^{N} y_i x_i_1

 w^T \sum_{i=1}^{N} x_i x_i_2 = \sum_{i=1}^{N} y_i x_i_2

 w^T \sum_{i=1}^{N} x_i x_i_D = \sum_{i=1}^{N} y_i x_i_D.

Напомним определение скалярного произведения:

a^T b = \sum_{i=1} ^ {N} a_i b_i.

Поэтому мы можем записать в матричной форме

w^T(X^TX)=y^TX

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

Поскольку нам нужно получить просто величину w, а не транспонированную, транспонируем обе части уравнения.

w^T(X^TX)=y^TX

[w^T(X^TX)]^T=[y^TX]^T

(X^TX)w=X^Ty.

Обратите внимание, что, например, при линейной зависимости Ax=b решение равно произведению b на величину, обратную А:

Ax=b\Rightarrow x\Rightarrow = A^{-1}b.

Соответственно, всё, что нам нужно сделать, – произвести ту же операцию для выделения w.

(X^TX)w=X^Ty

w=(X^TX)^{-1} X^Ty.

Это и есть решение для множественной линейной регрессии.

Как и в случае одномерной линейной регрессии, NumPy позволяет ускорить вычисления по сравнению с прямыми расчётами. Операции по нахождению решения линейной регрессии настолько распространены, что в NumPy для этого есть специальная функция solve. Поэтому всё, что нам нужно, – это вызвать её в коде и задать значения коэффициентов A и b.

Ax=b\Rightarrow x=np.linalg.solve(A, b)

w=(X^TX)^{-1} X^Ty = \; \Rightarrow w = np.linalg.solve(X^TX, X^Ty).

Помните, что эта статья – альтернативная ,которая рассматривалась в первой статье.

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

Решение множественной линейной регрессии с помощью матриц

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

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

 J = \sum_{i=1}^{N}(t_i - y_i)^2,

где t – матрица фактических результатов размерности Nx1, y – матрица прогнозируемых данных размерности Nx1.

y = Xw,

где X – матрица размерности NxD, w – вектор размерности Dx1, поэтому при их перемножении получим результат размерности Nx1.

Множественная регрессия и полиномиальная регрессия

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

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

Сначала раскроем выражение с J.

J=(t-y)^T(t-y) = t^Tt-t^Ty-y^Tt+y^Ty = t^Tt-t^TXw-(Xw)^Tt+(Xw)^T(Xw).

Далее находим производную и приравниваем её к нулю.

 \frac{\partial J}{\partial w}= - 2 X^T t + 2X^T Xw = 0

X^TXw=X^Tt

И находим наше решение.

w=(X^tX)=X^Tt.

Решение множественной регрессии в коде языка Python

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

У нас уже есть данные. Таким образом, наша матрица X имеет размерность два. Количество наблюдений равно 100. Следовательно, имеем D =2, N =100.

Создадим новый файл и назовём его lr_2d.py. Далее импортируем необходимые библиотеки. Обратите внимание, что кроме уже привычных NumPy и Matplotlib, нам также понадобится Mplot3d для построения трёхмерной диаграммы.

Matplotlib — библиотека на языке программирования Python для визуализации данных двумерной (2D) графикой (3D графика также поддерживается). Получаемые изображения могут быть использованы в качестве иллюстраций в публикациях.

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

Прежде всего, как и ранее, мы загружаем данные. Создаём пустые списки X и Y и загружаем данные из файла. Как вы помните, у нас также свободный член, характеризующий точку пересечения с осью ординат, поэтому создадим также наше x0 и устанавливаем для него значение 1.

X=[]

Y=[]

for line in open(‘data_2d.csv’):

x1, x2, y = line.split(‘,’)

X.append([float(x1), float(x2), 1])

Y.append(float(y))

Вы уже видели, насколько полезна библиотека NumPy, поэтому преобразуем наши X и Y в массивы NumPy.

X = np.array(X)

Y = np.array(Y)

Теперь изобразим наши данные графически, чтобы понять, как они выглядят. Это будет чуть сложнее по сравнению с предыдущими двухмерными изображениями, поскольку теперь у нас три измерения. Будьте внимательны! Если вы работаете под Mac, это может не работать. Я покажу вам, как это выглядит у меня под Linux.

fig = plt.figure()

ax = fig.add_subplot(111, projection = ‘3d’)

ax.scatter(X[:,0], X[:, 1], Y)

plt.show()

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

w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, Y)

Yhat = np.dot(X, w)

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

d1 = Y – Yhat

d2 = Y – Y.mean()

r2 = 1 – d1.dot(d1) / d2.dot(d2)

print ‘’r-squared:’’, r2

И запускаем программу. R2 у нас получилось очень неплохое – 0,998. Это потому, что мы так подобрали данные.

Полиномиальная регрессия – расширение линейной регрессии с кодом на языке Python

Давайте поговорим о полиномиальной регрессии.

Итак, короткие пояснения перед тем, как начнём писать код. Даже если у нас есть многочлен, скажем, квадратный или кубический, или какая-либо другая функция кривой, мы по-прежнему можем использовать линейную регрессию, несмотря на её название. Дело в том, что линейность регрессии означает, что коэффициенты модели являются линейными, но она ничего не говорит о том, каким должен быть сам x. Если вы возьмёте квадратичную функцию или функцию вида ax2+bx+c, вы увидите, что результат y линейно зависит от коэффициентов a, b и c, а потому к ним также применима линейная регрессия.

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

Итак, давайте создадим новый файл, который назовём lr_poly.py.

Как и прежде, импортируем библиотеки NumPy и Matplotlib.

import numpy as np

import matplotlib.pyplot as plt

Если вы посмотрите на данные, находящиеся в файле data_poly.csv, вы увидите всего два столбца данных, как и в файле data_1d.csv. Итак, загружаем наши данные.

X=[]

Y=[]

for line in open(‘data_poly.csv’):

x, y = line.split(‘,’)

Разница в том, что у нас теперь будет x, являющееся просто числом, и X, включающее, кроме x, ещё и x2.

x = float(x)

X.append([1, x, x*x])

Y.append(float(y))

Превращаем, как и раньше, в NumPy, массивы:

X = np.array(X)

Y = np.array(Y)

Теперь изобразим это в графическом виде.

plt.scatter(X[:, 1], Y)

plt.show

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

w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, Y)

Yhat = np.dot(X, w)

И вновь изобразим это графически.

plt.scatter(X[:, 1], Y)

plt.plot(X[:, 1], Yhat)

plt.show()

Получается нечто непонятное. Такое иногда бывает, когда X неотсортирован, – тогда каждая пара точек соединяется линией. Исправим это, переписав последнюю часть кода.

plt.scatter(X[:, 1], Y)

plt.plot(sorted(X[:, 1]), sorted(Yhat))

plt.show()

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

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

d1 = Y – Yhat

d2 = Y – Y.mean()

r2 = 1 – d1.dot(d1) / d2.dot(d2)

print ‘’r-squared:’’, r2

Коэффициент детерминации равен 0,999 – очень хорошее значение.

Прогнозирование систолического артериального давления по возрасту и весу

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

Просто скопируйте и вставьте этот адрес в свой браузер и загрузите Excel-файл. Мы будем использовать библиотеку pandas для чтения Access-файлов, но нам потребуется загрузить внешнюю библиотеку. Вы должны будете запустить

sudo pip install xlrd

что позволит использовать функцию pd dot для чтения Excel-файла.

Загрузив данные, вы увидите, что у нас есть три переменных X1, X2 и X3 для каждого пациента. Первая колонка X1 – это систолическое кровяное давление. Это и будет наше Y, значение которого мы будем прогнозировать. X2 – это наш входной фактор, характеризующий возраст в годах – мы будем анализировать кровяное давление как функцию от возраста. Кроме того, мы хотим проанализировать зависимость давления от веса в фунтах – это наше X3. Мы увидим, что между всеми этими данными есть линейная зависимость.

Начнём с загрузки необходимых библиотек. Импортируем matplotlib как plt, numpy как np и pandas как pd.

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

Используем также функцию чтения Excel-файла из библиотеки pandas. Файл называется mlr02.xls.

df = pd.read_excel(‘mlr02.xls’)

Вы также можете взять этот файл в Riehle.

Первое, что мы сделаем – преобразуем данные в метрическую систему. Неплохо также сразу и отобразить их графически. По оси x будем откладывать прожитые годы, а по оси y – артериальное давление.

X = df.as_matrix()

plt.scatter(X[:, 1], X[:, 0])

plt.show()

То же самое делаем и со вторым фактором, но теперь по оси x будем откладывать вес, а по оси y – давление.

plt.scatter(X[:, 2], X[:, 0])

plt.show()

Мы построим три линейные регрессии. Первая – зависимость исходящей переменной Y только от фактора X2, вторая – зависимость только от фактора X3, и третья – зависимость от обоих факторов, а также вычислим коэффициент детерминации R2. Используем для этого оператор цикла.

df[‘ones’] = 1

Y = df[‘X1’]

X = df[[‘X2’, ‘X3’, ‘ones’]]

X2only = df[[‘X2’, ‘ones’]]

X3only = df[[‘X3’, ‘ones’]]

def get_r2(X, Y):

w = np.linalg.solve( X.T.dot(X), X.T.dot(Y) )

Yhat = X.dot(w)

d1 = Y – Yhat

d2 = Y – Y.mean()

r2 = 1 – d1.dot(d1) / d2.dot(d2)

return r2

print ‘’r2 for x2 only:’’, get_r2(X2only, Y)

print ‘’r2 for x3 only:’’, get_r2(X3only, Y)

print ‘’r2 for both:’’, get_r2(X, Y)

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

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

Понравилась статья? Поделить с друзьями:
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: