<jupyter_start><jupyter_text>Линейная регрессияПри изучении задач машинного обучения удобно отталкиваться от реальных данных и реальных бизнес проблем. В качестве задачи для обсуждения метода линейной регрессии мы рассмотрим сильно учененный набор данных из конкурса kaggle [House Sales in King County, USA](https://www.kaggle.com/harlfoxem/housesalesprediction).Перед нами стоит задача по предсказанию цены за дом в зависимости от различных его параметров. В этом наборе данных мы для простоты специально оставим из всех парамтеров только жилую площадь в квадратных футах.  Загрузка данныхВоспользуемся для заргузки данных из CSV файла библиотекой pandas.<jupyter_code>import pandas as pd

housesales = pd.read_csv('../datasets/kc_house_data_reduced.csv')
housesales['price'] = housesales['price']/100000
housesales['sqft_living'] = housesales['sqft_living']/1000<jupyter_output><empty_output><jupyter_text>Исследование данныхСледующим шагом будет исследование данных, с которыми нам предстоит иметь дело. Самым мощным инструментов для этого является визуализация. Посмотим как цена дома зависит от его жилой площади<jupyter_code>%pylab inline
import seaborn as sns  # для красивого стиля графиков
#from matplotlib.ticker import FuncFormatter # для форматирования разметки по оси y
ax = housesales.plot.scatter(x='sqft_living', y='price')
ax.set_xlabel('Жилая площадь, 1k футы$^2$')
ax.set_ylabel('Цена, 100k $')<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Мы видим, что точек не очень много и между жилой площадью и ценой определенно нет простой зависимости. При этом вполне логично, что есть характерная закономерность: чем больше жилая площадь, тем больше цена. В качестве более сильного предположения можно выдвинуть следующее. Если площадь увеличивается в 2 раза, то и цена увеличивается примерно в 2 раза. Это означает, что мы можем попробовать описать наблюдаемую *зависимость* как *прямо пропорциональную*, и, что еще важнее, как *линейную*.  Постановка задачиВажнейшим элементом *любого* исследования данных является постановка формальной задачи, которую нам предстоит решить. В данном примере она напрашивается естественным образом: *по заданной площади оценить стоимость дома*.Например, для площади 2000 квадратных футов цена примерно может составить 500,000$.<jupyter_code>ax = housesales.plot.scatter(x='sqft_living', y='price')
ax.set_xlabel('Жилая площадь, 1k футы$^2$')
ax.set_ylabel('Цена, 100k $')
ax.axvline(x=2, ymax=0.45, color='green', linestyle='--')
ax.axhline(y=5, xmax=0.5, color='red', linestyle='--')
scatter(2, 5, color='red', marker='x', s=40, linewidth=3)<jupyter_output><empty_output><jupyter_text>Признаки и целевая переменнаяПредставим наши данные в табличном виде<jupyter_code>housesales.head(10)<jupyter_output><empty_output><jupyter_text>Каждоя строка состоит из двух чисел: площадь дома и его цена. Обозначим число строка как $n$<jupyter_code>n = len(housesales)
n<jupyter_output><empty_output><jupyter_text>Будем называть признаками те параметры, которые нам будут в задаче даны и целевыми переменными те, которые требуется определить. В данном случае у нас один признак - площадь дома и одна целевая переменная - цена. Каждую строку таблицы будем называть объектом.Обозначим за $x^{(i)}$ - значение признака для $i$-го объекта, а за $y^{(i)}$ - значение целевой переменной для $i$-го объекта. Индекс $i$ пробегает значения от $1$ до $n$.Пару $(x^{(i)}, y^{(i)})$ будем называть элементом обучающей выборки. A весь набор таких пар обучающей выборкой и обозначать ее как $(X, Y)$.Следующим важным элементом является функция $h_\theta(x)$, которую называют гипотезой или моделью. В данной задаче, если на вход этой функции подать площадь дома, то на выходе должны быть вычислена примерная стоимость. Можно грубо сказать, что все методы машинного обучения направленны на решение одного вопроса: как найти функцию $h$?Суть любой задачи машинного обучения сводится к схеме приведенной ниже **TODO**: вставить схему  Линейная регрессияВ методе линейной регресии предлагается представить гипотезу в виде$$h_\theta(x) = \theta_0 + \theta_1 x.$$Здесь $\theta_0$ и $\theta_1$ - неизвестные вещественнозначные параметры, которые нам предстоит определить (выучить).Рассмотрим несколько примеров для различных параметров $\theta_i$. Для простоты попложим $\theta_0= 0$.<jupyter_code>ax = housesales.plot.scatter(x='sqft_living', y='price')
ax.set_xlabel('Жилая площадь, 1k футы$^2$')
ax.set_ylabel('Цена, 100k $')

xx = np.linspace(1, 3, 2)
plot(xx, 1 * xx, label='$\\theta_1 = 1$')
plot(xx, 2 * xx, label='$\\theta_1 = 2$')
plot(xx, 3 * xx, label='$\\theta_1 = 3$')
ax.legend(loc=4)<jupyter_output><empty_output><jupyter_text>Функция потерьЧтобы выбрать параметры $\theta_i$ нам необходим числовой критерий $L(\theta_0, \theta_1)$. Допустим мы хотим придумать такой критерий, чтобы чем его значение меньше, тем лучшие значения параметров гипотезы мы выбрали. При этом чтобы этот критейри не принимал отрицательных значений. Это означает, что если критерий равен нулю, то мы считаем, что выбрали наилучшие из возможных значений параметров.В методе линейно регресси в качестве такого критерий предлагается рассмотреть среднее значение квадрата ошибки.$$L(\theta_0, \theta_1) = \frac{1}{2n} \sum\limits_{i=1}^{n}\left(h_\theta(x^{(i)}) - y^{(i)}\right)^2. $$$L$ также называют функцией потерь.  Интуиция за функцией потерьЧтобы пролить свет на выбор имеено такого вида функции $L$ рассмотрим упрощенный пример задачи регрессии. Обучающая выборка представлена на рисунке, а гипотеза $h$ имеет упрощенный вид$$h_\theta(x) = \theta_1 x$$<jupyter_code># Обучающая выборка
x = np.array([1, 2, 3])
y = np.array([1, 2, 3])

# Будем строить график h(x) и L(theta)
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

# Отметим на левом графике обучающую выборку с помощью зеленых крестиков
ax1.scatter(x, y, s=100, linewidth=3, marker='x', color='green')

# И построим график гипотезы при трех различных значениях theta
theta_1 = 0.75
theta_2 = 1
theta_3 = 1.3

xx = np.linspace(0.5, 3.5, 10)
ax1.plot(xx, theta_1*xx, color='red', label='$\\theta_1={}$'.format(theta_1))
ax1.plot(xx, theta_2*xx, color='blue', label='$\\theta_1={}$'.format(theta_2))
ax1.plot(xx, theta_3*xx, color='green', label='$\\theta_1={}$'.format(theta_3))
ax1.legend(loc=4)
ax1.set_xlim((0, 3.5))
ax1.set_ylim((0, 3.5))
ax1.set_xlabel('$x$')
ax1.set_ylabel('$y$')

# Определим функцию L для вычисления функции потерь
def L(x, y, theta_1):
    n = len(x)
    return 1/(2*n)*np.sum((theta_1*x - y)**2)

# И построим график функции потерь в зависимости от различных значений параметра theta_1
ax2.scatter(theta_1, L(x, y, theta_1), color='red', s=100)
ax2.scatter(theta_2, L(x, y, theta_2), color='blue', s=100)
ax2.scatter(theta_3, L(x, y, theta_3), color='green', s=100)
theta = np.linspace(0.6, 1.4, 20)
L_theta = [L(x, y, theta_1) for theta_1 in theta]
_ = ax2.plot(theta, L_theta)
ax2.set_xlabel('$\\theta$')
ax2.set_ylabel('$L(\\theta)$')<jupyter_output><empty_output><jupyter_text>Как можно увидеть из правого графика, благодря тому, что функция $L$ имеет квадратичных вид, то у нее есть ровно один минимум, которых автоматически является глобальным. В данном примере этот минимум достугается в точке $\theta_1 = 1$.  Метод градиентного спука Будучи уверены, что функция $L$ имеет ровно один минимум мы можем сформулировать задачу поиска наилучшей гипотезы как задачу оптимизации, в которой необходимо минимизировать функцию потерь $L$ по параметрам гипотезы $\theta_0, \theta_1$.\begin{equation}\theta_0^*, \theta_1^* = \underset{\theta_0, \theta_1 \in \mathbb{R}}{\mathrm{argmin}} L(\theta_0, \theta_1).\end{equation}В данной поставноке мы должны минимизировать функцию сразу по двум переменным. Построим на основе нашей обучающей выборки поверхность $L(\theta_0, \theta_1)$.<jupyter_code>def L(x, y, theta_0, theta_1):
    n = len(x)
    h = theta_0 + theta_1*x
    return 1/(2*n)*np.sum((h - y)**2)

X = housesales['sqft_living']
y = housesales['price']

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

theta_0 = np.linspace(-30, 30, 100)
theta_1 = np.linspace(-10, 10, 100)
L_theta = np.array([[L(X, y, theta_0i, theta_1j) for theta_0i in theta_0] for theta_1j in theta_1])
theta_0v, theta_1v = np.meshgrid(theta_0, theta_1)

ax.plot_surface(theta_0v, theta_1v, L_theta)
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.title('Поверхность функции потерь $L(\\theta_0, \\theta_1)$')<jupyter_output><empty_output><jupyter_text>Линии уровеняДля визуального анализа график поверхности не всегда является удобным способом представления функции потерь. Для случая функции 2-х переменных удобно рассмотреть ее графическое представление в виде линий уровня - уривых, на которых функция принимает одинаковое значение.<jupyter_code>levels = [2, 5, 10, 50, 100]
CS = plt.contour(theta_0v, theta_1v, L_theta, levels, colors='k')
plt.clabel(CS, inline=1, fontsize=10, colors='k')
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.title('Линии уровня функции потерь $L(\\theta_0, \\theta_1)$')<jupyter_output><empty_output><jupyter_text>Из приведенного выше графика хорошо видно, что у функции $L$ есть всего одна точка экстремума, которая соответсвует глобальному минимум. Важный вопрос, ответ на который нам предстоит рассмотреть, заключается в том, как найти эту точку минимума. Т.е. решить задачу оптимизации или по терминологии машинного обучения - обучить модель линейной регресии.В качестве основного метода решения мы будет использовать градиентный спуск. Напомним, что под градиентом функции в точке понимается вектор, который указывает направление наискорейшего роста функции. Компоненты вектора градиента определяеются как частные производные.$$grad~ L = \left(\frac{\partial L}{\partial \theta_0}, \frac{\partial L}{\partial \theta_1}\right).$$Для пояснение того, как метод градиентного спуска позволяет найти минимум функции, мы вернемся к одномерному случаю, когда зафиксируем $\theta_0 = 0$ и будем рассматривать зависимость функции потерь $L$ только от $\theta_1$. График этой зависимости для задачи оценки стоимост жилья приведен на рисунке ниже.Суть метода градиентного спуска заключается в следующем. На каждом шаге мы обновляем оценку для аргумента минимума $\theta_1^{(k)}$, начиная с некоторого начального значения $\theta_1^{(0)}$ по правилу$$\theta_1^{(k + 1)} = \theta_1^{(k)} - \alpha \frac{\partial L}{\partial \theta_1}(\theta_1^{(k)}). $$Здесь $\alpha$ - параметр, который разывают скоростью обучения (learning rate). Удачный выбор значения $\alpha$ может полохительно сказаться на числе шагов алгоритма, которые необходимы для приближения к значения минимума.Попробуем с помощью визуализации пояснить суть работы алгоритма.Пусть $\theta_1^{(0)} = 8$. Поскольку производная в точке описывает тангенс угла наклона касательной к функции в этой точке, то как видно из рисунка, производная будет иметь положительный знак$$\frac{\partial L}{\partial \theta_1}(\theta_1^{(0)}) > 0.$$А следовательно, производная взятая с отрицательным знаком будет направлена в сторону уменьшения значения функции $L$ и это действительно то направление, в котором следует искать минимум.<jupyter_code>theta_1 = np.linspace(-5, 9, 100)
L_theta_1 = [L(X, y, 0, theta_1i) for theta_1i in theta_1]
plot(theta_1, L_theta_1)
theta1_0 = 8
scatter(theta1_0, L(X, y, 0, theta1_0), s=50)

def gradL(X, y, theta_0, theta_1):
    n = len(X)
    h = theta_0 + theta_1*X.values
    return 1/(n)*np.sum(X.transpose()@(h - y.values))

arrow(theta1_0, L(X, y, 0, theta1_0), -1, -gradL(X, y, 0, theta1_0), linewidth=2, head_width=0.5, head_length=2, fc='k', ec='k')
vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
vlines(x=theta1_0 - 1, ymin=-20, ymax=L(X, y, 0, theta1_0 -1), color='green', alpha=0.3, linestyles='--')
ylim(-20, 100)
xlabel('$\\theta_1$')
ylabel('$L(0, \\theta_1)$')
text(x=theta1_0-0.45, y = -35, text='$\\theta_1^{(0)}$', s=50)
text(x=theta1_0-1.45, y = -35, text='$\\theta_1^{(1)}$', s=50)<jupyter_output><empty_output><jupyter_text>Рассмотрим другую ситуацию. Пусть $\theta_1^{(0)} = -2$, тогда производная в этой точке будет иметь отрицательный знак$$\frac{\partial L}{\partial \theta_1}(\theta_1^{(0)}) < 0.$$А следовательно взятая с обратным знаком будет подталкивать значение $\theta_1^{(i)}$ вправо, что в таком случае и требуется.<jupyter_code>plot(theta_1, L_theta_1)
theta1_0 = -2
scatter(theta1_0, L(X, y, 0, theta1_0), s=50)

def gradL(X, y, theta_0, theta_1):
    n = len(X)
    h = theta_0 + theta_1*X.values
    return 1/(n)*np.sum(X.transpose()@(h - y.values))

arrow(theta1_0, L(X, y, 0, theta1_0), 1, gradL(X, y, 0, theta1_0), linewidth=2, head_width=0.5, head_length=2, fc='k', ec='k')
vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
vlines(x=theta1_0 + 1, ymin=-20, ymax=L(X, y, 0, theta1_0 +1), color='green', alpha=0.3, linestyles='--')
ylim(-20, 100)
xlabel('$\\theta_1$')
ylabel('$L(0, \\theta_1)$')
text(x=theta1_0-0.45, y = -35, text='$\\theta_1^{(0)}$', s=50)
text(x=theta1_0+0.55, y = -35, text='$\\theta_1^{(1)}$', s=50)<jupyter_output><empty_output><jupyter_text>Повторяя шаг алгорима многократно мы постеменно будем приближаться к точке минимума функции как это показано на рисунке ниже.<jupyter_code>plot(theta_1, L_theta_1)
theta1_0 = -2
scatter(theta1_0, L(X, y, 0, theta1_0), s=50)

alpha = 0.1
for i in range(8):
    arrow(theta1_0, L(X, y, 0, theta1_0),
          -alpha*gradL(X, y, 0, theta1_0), -alpha*gradL(X, y, 0, theta1_0)**2,
          linewidth=2, head_width=0.5, head_length=2, fc='k', ec='k')
    vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
    text(x=theta1_0-0.45, y = -35, text='$\\theta_1^{(%d)}$' % i, s=50)
    theta1_0 += -alpha*gradL(X, y, 0, theta1_0)
vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
ylim(-20, 100)
xlabel('$\\theta_1$')
ylabel('$L(0, \\theta_1)$')
title('$\\alpha = {}$'.format(alpha))<jupyter_output><empty_output><jupyter_text>Отметим следующее наблюдение, которое хорошо видно на графике выше. Чем ближе мы к минимуму, тем меньше изменяется значение параметра $\theta_1^{(i)}$ и это при том, что скорость обучения постоянна $\alpha = const$. Объяснить этот эффект можно следующим образом. В точке экстремума касательная к функции горизонтальна и проивзодная равна нулю. Следовательно, чем ближе мы к точке оптимума, тем меньше обсолютное значение производной и тем меньше изменяется параметр $\theta_1^{(i)}$. В некоторых случаях это может создавать проблемы с необходимостью использовать очень большое число итераций для достижения результата. В этом случае следует увеличть скорость обучения обучения $\alpha$.<jupyter_code>plot(theta_1, L_theta_1)
theta1_0 = -2
scatter(theta1_0, L(X, y, 0, theta1_0), s=50)

alpha = 0.2
for i in range(3):
    arrow(theta1_0, L(X, y, 0, theta1_0),
          -alpha*gradL(X, y, 0, theta1_0), L(X, y, 0, theta1_0 - alpha*gradL(X, y, 0, theta1_0)) - L(X, y, 0, theta1_0),
          linewidth=2, head_width=0.5, head_length=2, fc='k', ec='k')
    vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
    text(x=theta1_0-0.45, y = -35, text='$\\theta_1^{(%d)}$' % i, s=50)
    theta1_0 += -alpha*gradL(X, y, 0, theta1_0)
vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
ylim(-20, 100)
xlabel('$\\theta_1$')
ylabel('$L(0, \\theta_1)$')
title('$\\alpha = {}$'.format(alpha))<jupyter_output><empty_output><jupyter_text>При это следует учитывать, что слишком большое значение $\alpha$ может приводить к тому, что мы будем постоянно перепрыгивать через локальный минимум и алгоритм может не сойтись к искомой точке.<jupyter_code>plot(theta_1, L_theta_1)
theta1_0 = -2
scatter(theta1_0, L(X, y, 0, theta1_0), s=50)

alpha = 0.63
for i in range(6):
    arrow(theta1_0, L(X, y, 0, theta1_0),
          -alpha*gradL(X, y, 0, theta1_0), L(X, y, 0, theta1_0 - alpha*gradL(X, y, 0, theta1_0)) - L(X, y, 0, theta1_0),
          linewidth=2, head_width=1, head_length=0.1, fc='k', ec='k')
    vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
    text(x=theta1_0-0.45, y = -35, text='$\\theta_1^{(%d)}$' % i, s=50)
    theta1_0 += -alpha*gradL(X, y, 0, theta1_0)
vlines(x=theta1_0, ymin=-20, ymax=L(X, y, 0, theta1_0), color='green', alpha=0.3, linestyles='--')
ylim(-20, 100)
xlabel('$\\theta_1$')
ylabel('$L(0, \\theta_1)$')
title('$\\alpha = {}$'.format(alpha))<jupyter_output><empty_output><jupyter_text>Рекомендуется использовать следующие значения для скорости обучения в порядке их возрастания.$$\alpha = 0.01, 0.03, 0.1, 0.3, 1.$$ Вернемся к двумерному случаю и проиллюстрируем работу метода градиентного спуска с помощью линий уровня.<jupyter_code>CS = plt.contour(theta_0v, theta_1v, L_theta, levels, colors='k')
plt.clabel(CS, inline=1, fontsize=10, colors='k')
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.title('Линии уровня функции потерь $L(\\theta_0, \\theta_1)$')

def gradL(X, y, theta_0, theta_1):
    n = len(X)
    h = theta_0 + theta_1*X.values
    return (1/(n)*np.sum(h - y.values),
            1/(n)*np.sum(X.transpose()@(h - y.values)))

theta0_i = 20
theta1_i = 5

alpha = 0.1
for i in range(1000):
    upd_theta0_i, upd_theta1_i = gradL(X, y, theta0_i, theta1_i)
    if i < 10 or i % 50 == 0:
        scatter(theta0_i, theta1_i, s=20)
    theta0_i -= alpha*upd_theta0_i
    theta1_i -= alpha*upd_theta1_i<jupyter_output><empty_output><jupyter_text>Выражения для градиентовДля реализации метода нам необходимо в явном виде вычислить частные производные функции потерь.\begin{eqnarray}\frac{\partial L}{\partial \theta_0} & = & \frac{1}{n}\sum\limits_{i=1}^n (h_\theta(x^{(i)}) - y^{(i)}) \\\frac{\partial L}{\partial \theta_1} & = & \frac{1}{n}\sum\limits_{i=1}^n (h_\theta(x^{(i)}) - y^{(i)})x^{(i)} \\\end{eqnarray}Здесь использовано правило вычисления производной для составной функции.Используя данное представление мы уже нашли выше оценку для минимума функции $L$. Координаты этой точки равны<jupyter_code>print("theta_0* = {:g}, theta_1* = {:g}".format(theta0_i, theta1_i))<jupyter_output>theta_0* = 0.431547, theta_1* = 2.16905<jupyter_text>Функция потерь при данном значении параметров принимает значение<jupyter_code>L(X, y, theta0_i, theta1_i)<jupyter_output><empty_output><jupyter_text>Построим график для полученной гипотезы<jupyter_code>ax = housesales.plot.scatter(x='sqft_living', y='price')
ax.set_xlabel('Жилая площадь, 1k футы$^2$')
ax.set_ylabel('Цена, 100k $')

xx = np.linspace(1, 3, 2)
plot(xx, theta0_i + theta1_i * xx)<jupyter_output><empty_output><jupyter_text>И вернемся к вопросу, который мы можем задать этой модели. Какова примерная стоимость дома площадью 2000 квадратных футов?<jupyter_code>print("Примерная стоимость: {:.2g} * 100k$".format((theta0_i + theta1_i * 2)))<jupyter_output>Примерная стоимость: 4.8 * 100k$<jupyter_text>Основные результатыПоздравляю! С помощью метода линейной регресии мы смогли построить систему машинного обучения, которая на основе предоставленных данных делает оценку для стоимости жилья. Для закрепления материала приведем еще раз основные понятия* задача регрессии (regression problem) - задача машинного обучения, в которой по заданому входу нужно предсказать значение из множества вещественных чисел.* объекты (objects)- сущности, для которых алгоритм делает предсказания* признаки (features) - числовое представление объектов* матрица объекты-признаки  - матрица по строкам которой отложены объекты, а по столбцам - признаки* целевая  переменная (target) - значение, которое необходимо предсказать* обучающая выборка (training set)- множество пар (объект, значение целевой переменной), которые используются для выбора модели* линейная регрессия (linear regression) - метод решения задачи регрессии* гипотеза (hypothesis) - вид функции, которая используется для нахождения целевой переменной* параметры гипотезы (hypothesis parameters) - неизвестные величины, которые необходимо найти в ходе обучения* функция потерь (loss function) - штраф за неправильный результат* градиентный спуск (gradient descent) - метод решения задачи оптимизации, который базируется на использовании градиента функции* срокость обучения (learning rate) - основной параметр метода градиентного спуска, который отвечает за скорость сходимостиДалее мы рассмотрим как динейную регресси можно обобщить на случай произвольной размерности и что можно сделать с явно нелинейными зависимостями.  Многомерная линейная регрессия В случае, если целевай переменная зависит от нескольких признаков, например, площадь дома и цисло ванных комнат, линейную регрессию можно обобщить на случай многмерного вектора признаков.Рассмотрим задачу прогноза цены жилья в зависимости от следующих параметров: * количество спальных комнат * количество ванных комнат * жилая площадь * общая площадь * количество этажей * год постройки<jupyter_code>housesales = pd.read_csv('../datasets/house_price/train.csv.gz')[['SalePrice',
                                                              'LotArea',
                                                              'LotFrontage',
                                                              '1stFlrSF',
                                                              'OverallQual',
                                                              'GarageArea',
                                                              'YearBuilt',
                                                              'TotalBsmtSF',
                                                              '2ndFlrSF',]].dropna()
housesales['SalePrice'] = housesales['SalePrice']/100000
housesales.head()<jupyter_output><empty_output><jupyter_text>Предствим матрицу объекты-признаки как$$X = \begin{pmatrix}x_{11} & x_{12} & \cdots & x_{1m} \\x_{21} & x_{22} & \cdots & x_{2m}  \\\vdots & \vdots & \ddots & \vdots \\x_{n1} & x_{n2} & \cdots & x_{nm} \\\end{pmatrix}$$Здесь $n$ - число объектов и $m$ - число признаков.Вектор целевых переменных $y$ имеет ту же форму, что и для одного признака$$y = \begin{pmatrix}y_1 \\y_2 \\\vdots \\y_n\end{pmatrix}$$Для многомерного случая гипотеза $h_\theta(x)$ пример следующий вид$$h_\theta(x) = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_mx_m.$$Из приведенного веше выражения видно, что слагаемое $\theta_0$ отличается от всех остальных, которые представляют собой произведения различных компанент вектора параметров $\theta$ и вектора признаков $x = (x_1, x_2 , \dots, x_m)$. Чтобы привести выражение для гипотезы к более простому виду мы применим следующий приём. Расширим матрицу объекты-признаки дополнительным столбцом слева, состоящим из одних единиц$$X = \begin{pmatrix}1 & x_{11} & x_{12} & \cdots & x_{1m} \\1 & x_{21} & x_{22} & \cdots & x_{2m}  \\\vdots & \vdots & & \ddots & \vdots \\1 & x_{n1} & x_{n2} & \cdots & x_{nm} \\\end{pmatrix}$$и используем нулевой индекс для обозначения этого столбца ($x_{i0} = 1, \quad i = \overline{1, n}$).Тогда выражение для гипотезы $h_\theta(x)$ можно переписать в виде$$h_\theta(x) = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_mx_m = \theta_0x_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_mx_m = \sum\limits_{i = 0}^{m}\theta_ix_i = \theta^Tx.$$Т.е. представить в виде скалярного произведения ветора параметров $\theta$ и вектора признаков $x$. Здесь и далее под матрицей объектов-признаков $X$ будем понимать расширенную матрицу с единичным первым столбцом.В качестве функции потерь $L$ как и прежде рассмотрим средний квадрат ошибки$$L(\theta) = \frac{1}{2n}\sum\limits_{i=1}^{n}\left(h_\theta(x^{(i)}) - y^{(i)}\right)^2.$$Здесь с помощью $x^{(i)} = (x_{i0}, x_{i1}, \dots, x_{im})^T$ обозначена вектор, составленый из строки матрицы $X$.  Векторная форма градиентного спускаРассмотрим как метод градиентного спуска можно обобщить на многомерный случай. Для этого вычислим градиент фукнци потерь $J$ по всем компонентам вектора параметров $\theta$.\begin{eqnarray}\frac{\partial L}{\partial \theta_0} & = & \frac{1}{n}\sum\limits_{i=1}^n (h_\theta(x^{(i)}) - y^{(i)}) \\\frac{\partial L}{\partial \theta_j} & = & \frac{1}{n}\sum\limits_{i=1}^n (h_\theta(x^{(i)}) - y^{(i)})x_{ij}, \quad j = \overline{1, m} \\\end{eqnarray}Или, используя прием с дополнительным единичным столбцом ($x_{i0} = 1$), можно переписать выражением для компонент вектора градиента в единообразном стиле\begin{eqnarray}\frac{\partial L}{\partial \theta_j} & = & \frac{1}{n}\sum\limits_{i=1}^n (h_\theta(x^{(i)}) - y^{(i)})x_{ij}, \quad j = \overline{0, m} \\\end{eqnarray}Наконец, в векторном виде выражение для градиента пример вид$$grad~ L(\theta) = \frac{1}{n} X^T(X\theta - y).$$Уравнение для обновления параметров $\theta$ для шага градиентного спуска можно также представить в вектрном виде$$\theta^{(k+1)} = \theta^{(k)} - \alpha ~grad ~ L(\theta^{(k)}).$$  Кривая обученияВ многомерном случае не представляется возможным визуализировать работу метода градиентного спуска с помощью линий уровня, поскольку простанство параметров в общес случае имеет размерность $m+1$. На помощь нам приходит другое представление, которое показывается изменение значений функции потерь в точке $\theta^{(k)}$ в зависимости от номера шага $k$. Такой график называется кривой обучения.  Проведем предобработку данных<jupyter_code>n = len(housesales)
m = len(housesales.drop('SalePrice', axis=1).columns)
y = housesales['SalePrice'].values.reshape((n, 1))
X = housesales.drop('SalePrice', axis=1).values.reshape((n, m))
X = np.hstack((np.ones((n, 1)), X))<jupyter_output><empty_output><jupyter_text>Градиентный спуск<jupyter_code>def grad(y, X, theta):
    n = y.shape[0]
    return 1/n * X.transpose() @ (X @ theta - y)

def L(y, X, theta):
    n = y.shape[0]
    return 1/(2*n)*np.sum(np.power(X @ theta - y, 2))

def fit(y, X, theta_0, alpha=0.001, nsteps = 100):
    theta = np.copy(theta_0)
    loss = [L(y, X, theta)]
    for i in range(nsteps):
        theta -= alpha*grad(y, X, theta)
        loss.append(L(y, X, theta))
    return loss, theta
theta_0 = np.zeros((m + 1, 1))
loss_history, theta_star = fit(y, X, theta_0, alpha=1e-10, nsteps=10000)
plt.plot(loss_history)
plt.xlabel('$k$')
plt.ylabel('$L(\\theta^{(k)})$')
_ = plt.title('Кривая обучения')
L(y, X, theta_star)<jupyter_output><empty_output><jupyter_text>Стоит обратить внимание, что наблюдается очень медленная сходимость метода градиентного спука для такой постановки задачи. Попробуем решить эту проблему с помощью предобработчки данных, а конкретно - нормализации признаков.  Нормирование признаков Поскольку различные признаки могу иметь отличные размерности, то, для устойчисвости численной процедуры метода градиентного спуска, рекомендуется провести их нормализацию. Для этого каждый признак можно по отдельности нормировать, вычтя из него среднее значение и разделив на средне-квадратичное отклонение$$    U^{(j)} = \frac{X^{(j)} - M[X^{(j)}]}{\sqrt{D[X^{(j)}]}}, \quad j = \overline{1, m}.$$При этом первый столбце из одних единиц остается неизменным.$$    U^{(0)} = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1\end{pmatrix}$$<jupyter_code>U = np.ones((n, m + 1))
for j in range(1, m + 1):
    U[:, j] = (X[:, j] - np.mean(X[:, j]))/np.std(X[:, j])
    
loss_history, theta_star = fit(y, U, theta_0, alpha=1e-1, nsteps=50)
plt.plot(loss_history)
plt.xlabel('$k$')
plt.ylabel('$L(\\theta^{(k)})$')
_ = plt.title('Кривая обучения')<jupyter_output><empty_output><jupyter_text>Из графика видно, что при нормализации признаков удалось достичь того же значения функции потерь сделав в 10 раз меньше шагов. Значение функции потерь для нашей задачи составило<jupyter_code>L(y, U, theta_star)<jupyter_output><empty_output><jupyter_text>Как же можно улучшить этот результат оставаясь в рамках метода линейно регрессии?   Нелинейные зависимости Проведем анализ нашей задачи. Для этого построим зависимость целевой переменной от каждого из признаков по отдельности.<jupyter_code>sns.pairplot(housesales, y_vars="SalePrice",
                 x_vars=["LotArea", "LotFrontage", "1stFlrSF", "OverallQual"],
            )
sns.pairplot(housesales, y_vars="SalePrice",
                 x_vars=["GarageArea", "YearBuilt", "TotalBsmtSF", "2ndFlrSF"],
            )<jupyter_output><empty_output><jupyter_text>По характерному виду некоторых зависимостей можно предположить, что целевая переменная нелинейно зависит от некоторых признаков. Чтобы учеть подобные зависимости мы модем обогатить числовые признаки их квадратами и произведениями различных признаков, получив, таким образом, новую матрицу объекты-признаки.$$X = \begin{pmatrix}1 & x_{11} & x_{12} & \cdots & x_{1m} & x_{11}^2 & x_{11}x_{12} & \cdots & x_{1m}^2\\1 & x_{21} & x_{22} & \cdots & x_{2m} & x_{21}^2 & x_{21}x_{22} & \cdots & x_{2m}^2 \\\vdots & \vdots & & \ddots & \vdots \\1 & x_{n1} & x_{n2} & \cdots & x_{nm} & x_{n1}^2 & x_{n1}x_{n2} & \cdots & x_{nm}^2\\\end{pmatrix}$$После чего используем тот же метод линейно регрессии.<jupyter_code>n = len(housesales)
m = len(housesales.drop('SalePrice', axis=1).columns)
y = housesales['SalePrice'].values.reshape((n, 1))

tmpX = np.ones((n, m + 1))
X = np.zeros((n, (m + 1)**2))
tmpX[:, 1:] = housesales.drop('SalePrice', axis=1).values.reshape((n, m))
for i in range(n):
    X_i = tmpX[i, :] .reshape(1, -1)
    X[i, :] = (X_i.T @ X_i).reshape(-1)
U = np.ones(X.shape)
for j in range(1, 2*m + 1):
    U[:, j] = (X[:, j] - np.mean(X[:, j]))/np.std(X[:, j])<jupyter_output><empty_output><jupyter_text>Нормальное уравнение<jupyter_code>theta_0 = np.zeros((U.shape[1], 1))
loss_history, theta_star = fit(y, U, theta_0, alpha=1e-2, nsteps=500)
plt.plot(loss_history)
plt.ylim((0, loss_history[10]))
plt.xlabel('$k$')
plt.ylabel('$L(\\theta^{(k)})$')
_ = plt.title('Кривая обучения')<jupyter_output><empty_output><jupyter_text>Значение функции потерь после процедуры градиентного спуска равно<jupyter_code>L(y, U, theta_star)<jupyter_output><empty_output><jupyter_text>что немного улучшает результат, полученный без использования квадратичных признаков.  Нормальное уравнение Значение оптимальных параметров линейной регресси может быть получено аналитически в матричном виде и вычисленно без применения метода градиентного спуска\begin{equation}\theta = (X^TX)^{-1}X^T y.\end{equation}Вывод данного уравениея оставим в качестве самостоятельного задания.Данное выражение называется **нормальным уравнением** линейной регрессии.Его основная вычислительная сложность связана с необходимостью обращения матрицы, которая может быть плохо обусловлена. Это может сказаться на точности вычислений. Метод граидентного спуска наоборот обладает тем свойством, что при удачном выборе скрости обучения будет с каждым шагом сходиться к оптимальному значения параметров $\theta$.Для устойчивого обращения плохо обусловленной матрицы мы можем использовать алгоритм обращения Мура-Пенроуза (numpy.linalg.pinv)<jupyter_code>theta_n = np.linalg.pinv(X.T @ X) @ X.T @ y<jupyter_output><empty_output><jupyter_text>Значение функции потерь при параметрах, вычисленных с помощью нормального уравнения, равно<jupyter_code>L(y, X, theta_n)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1370161497052038, 'fraction_numerical': 0.031241989233529863, 'mean_word_length': 5.0529480217222655, 'pattern_counts': {'":': 0, '<': 120, '<?xml version=': 0, '>': 119, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '26642', 'license': 'Apache-2.0', 'n_tokens_mistral': 13408, 'n_tokens_neox': 14516, 'n_words': 3095, 'path': 'notebooks/LinearRegression.ipynb', 'repo_name': 'it-sec-std/ml-intro', 'size': '31208'}