<jupyter_start><jupyter_text>Линейная регрессия Курс "Модели и методы интеллектуального анализа данных"<NAME>., доцент кафедры ИТЗИ УрГУПС (Екатеринбург)<EMAIL>  Оглавление ВведениеМодель линейной регрессииЭффект переобученияПример 1. Построение линии регрессии "дата"-"температура" Пример 2: СсылкиЗадание для самостоятельной работы      Введение  Задача регрессии (предсказания) относится к типовым задачам машинного обучения с учителем. Задача заключается в том, чтобы научиться предсказывать значения целевой (объясняемой) переменной на основе знаний о значениях этой переменной при заданных значениях признаков (объясняющих переменных). Данная задача весьма часто встречается на практике. Назад к Оглавлению  Модель линейной регрессии  Линейная регрессия (Linear regression) — модель зависимости целевой (объясняемой) переменной от одной или нескольких других независимых (объясняющих) переменных (которые называются факторами или регрессорами) с использованием линейной функции $\vec{y}=A\vec{x}+\vec{B}$, где $\vec{y}$ - вектор целевых значений, $A$ - матрица весов для имеющихся значений признаков, $\vec{B}$ - вектор значений смещений, $\vec{x}$ - имеющиеся значения признаков.Матрица $A$ имеет размерность $m$ на $n$, где $m$ - количество признаков, $n$ - количество наблюдений признаков. Например, зависимость$$РОСТ = A \cdot ВЕС + B$$линейно связывает рост и вес человека, позволяя по значению веса предположить значение роста, для этого определяется константа (коэффициент, вес) $a$ и смещение $b$.Пример более сложной зависимости:$$ЦЕНА = A \cdot КАЧЕСТВО + B \cdot ВРЕМЯ + C \cdot ЗАТРАТЫ + D$$Также более сложные зависимости можно строить с помощью конструирования новых признаков с использованием существующих (рассматривая, например, степени известных признаков):$$y(x)=a x^2 + b x + c,$$в которой объясняющими переменными (признаками)  являются $x$ и $x^2$, а искомых коэффициентов три: $a$, $b$ и $c$. В примере ниже показано приближение множества точек моделью $y(x) = Ax + B$ (линейная зависимость, один признак $x$) и моделью $y(x)=Ax^2 + Bx + C$ (квадратичная зависимость, два признака $x^2$ и $x$).<jupyter_code>xs = np.linspace(0, 5, 20)
ys = xs**2 + np.random.randn(20)
plt.scatter(xs, ys, c='b')
plt.plot(xs, 5*xs, c='r', label='Приближение $y=Ax+B$')
plt.plot(xs, xs**2+xs, c='g', label='Приближение $y=Ax^2+Bx+c$')
plt.grid(True)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>В общем случае, если задан вектор признаков $\vec{x}=(x_1, \cdots, x_n)$ и вектор весов $\vec{w}=(w_0, w_1, \cdots, w_n)$, то их линейной комбинацией является выражение$$y = w_0 + x_1 w_1 + \cdots + w_n x_n,$$или, при добавлении единицы в вектор признаков $\vec{x}=(1, x_1, \cdots, x_n)$ в более коротком виде$$y=w^T x.$$ В самом простом случае (одна объясняющая переменная, один признак) модель линейной регрессии представляет собой линейную функцию $y(x)=Ax+b$. Графически это означает, что мы пытаемся построить прямую линию $y(x) = Ax+B$, которая проходит "как можно ближе" к таблично заданным точкам $(x_1,y_1), \cdots,(x_n,y_n)$.<jupyter_code>xs = np.linspace(0, 5, 10)
ys = xs**2 + np.random.randn(10)
plt.scatter(xs, ys, c='b')
plt.plot(xs, 5*xs-2, c='r')
for x, y in zip(xs, ys):
    plt.plot([x,x], [5*x-2,y], c='g')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Параметр, характеризующий ошибку модели, очевидно связан с длинной зеленых отрезков на рисунке (т.е. расстояний между известными заданными точками $y_i$ и значениями, предсказанными моделью $Ax_i+B$). Поэтому в качестве меры ошибки берется сумма среднеквадратичных отклонений, для наиболее оптимальной модели эта мера должна быть минимальна $$F(A,B)=\sqrt{\sum_{i=1}^{n}{\left[ Ax_i + B -y_i\right]^2}}\to \min$$ Решение этой задачи оптимизации можно получить в виде точной формулы. Для этого используется тот факт, что экстремум (минимум) достигается в точке, в которой частные производные целевой функции $F(A,B)$ равны нулю: $$\frac{\partial F(A,B)}{\partial A}=0,~\frac{\partial F(A,B)}{\partial B}=0.$$Из этих двух уравнений получается система линейных уравнений для нахождения коэффициентов $A$ и $B$:$$A \cdot n + B \cdot \sum_{i=1}^{n}{x_i}=\sum_{i=1}^{n}{y_i}$$$$A \sum_{i=1}^{n}{x_i} + B \cdot \sum_{i=1}^{n}{x_i ^2}=\sum_{i=1}^{n}{x_i y_i}.$$Из этой системы в явном виде находятся коэффициенты регрессии $A$ и $B$, что и дает решение задачи. Однако у такого подхода есть существенный минус - при решении СЛАУ методом Крамера или методом обратной матрицы необходимо вычислять определитель, что неэффективно для очень больших матриц (когда много признаков) или для матриц, определитель которых близок к нулю.Поэтому в большинстве практических задач более эффективными являются приближенные численные методы, которые позволяют найти решение с заданной точностью с меньшими ресурсами (например, градиентный метод). Назад к Оглавлению  Переобучение и регуляризация<jupyter_code>def f(x):
    r = []
    for elem in x:
        r.append(round(elem,2))
    return(r)

xs = np.sort(f(10*np.random.random(25)))
ys = xs**2 + np.random.randn(25)*10
fig, ax = plt.subplots(4, 4, figsize=(20,10))
for pos, degree in enumerate(range(1,17)):
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(xs.reshape(-1,1), ys)
    result = model.predict(xs.reshape(-1,1))
    ax[pos//4][pos%4].plot(xs, result, c='r')
    ax[pos//4][pos%4].plot(xs, xs**2, c='g')
    ax[pos//4][pos%4].scatter(xs, ys, c='b', alpha=0.3)
    ax[pos//4][pos%4].set_title("n={}".format(degree))
model = make_pipeline(PolynomialFeatures(50), LinearRegression())
model.fit(xs.reshape(-1,1), ys)
model.score(xs.reshape(-1,1), ys)<jupyter_output><empty_output><jupyter_text>Регуляризация используется для борьбы с переобучением. При $L_2$-регуляризации минимизируется функция:$$M(X, \vec{y}, \vec{w}) = L (X, \vec{y}, \vec{w}) + \lambda |\vec{w}|^2$$<jupyter_code>from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso<jupyter_output><empty_output><jupyter_text>Назад к Оглавлению  Пример 1. Построение линии регрессии "дата"-"температура"  Зададим массив значений $ts$ (значение переменной времени, моменты измерений) и $Ts$ (температура).<jupyter_code>ts = np.linspace(0, 10, 15)
Ts = ts + np.random.randn(15)
df = pd.DataFrame({"Время": ts, "Температура": Ts})
display(df.T)
df.plot()
plt.grid(True)
plt.show()
plt.boxplot(df.Температура, vert=False)
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Отобразим точки<jupyter_code>plt.plot(df['Время'], df['Температура'], 'ro')
plt.xlabel("Время")
plt.ylabel("Температура")
plt.grid(True)
plt.show()<jupyter_output>/usr/lib/python3/dist-packages/matplotlib/cbook/__init__.py:1377: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  x[:, None]
/usr/lib/python3/dist-packages/matplotlib/axes/_base.py:237: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  x = x[:, np.newaxis]
/usr/lib/python3/dist-packages/matplotlib/axes/_base.py:239: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  y = y[:, np.newaxis]<jupyter_text>Будем искать приближение в виде линейной функции $y(x)=Ax+b$.Как рассмотрено выше, для нахождения коэффициентов $A$ и $B$ решаем систему линейных уравнений:$$A \cdot n + B \cdot \sum_{i=1}^{n}{x_i}=\sum_{i=1}^{n}{y_i}$$$$A \sum_{i=1}^{n}{x_i} + B \cdot \sum_{i=1}^{n}{x_i ^2}=\sum_{i=1}^{n}{x_i y_i}.$$ Решим СЛАУ методом обратной матрицы<jupyter_code>xs = df.Время
ys = df.Температура
M = np.matrix([[len(xs), sum(xs)], [sum(xs), sum([i**2 for i in xs])]])
N = np.array([sum(ys), np.dot(xs, ys)])
#import scipy.linalg as spla
invM = spla.inv(M)
S = np.dot(N, invM)
print("Коэффициент A: {}, коэффициент B:{}.".format(round(S[1], 2), round(S[0], 2)))<jupyter_output>Коэффициент A: 0.99, коэффициент B:0.12.<jupyter_text>Построим график функции $y(x)=Ax+B$ совместно с заданными точками.Также построим предсказываемое значение для $Ж=15$ (желтый ромб в правом верхнем углу).<jupyter_code>Ж = 15
plt.plot(xs, ys,'ro')
plt.plot(xs, [S[1]*i+S[0] for i in xs], marker="o")
plt.plot(Ж, S[1]*Ж+S[0], marker='d')
plt.grid(True)
plt.show()
print("Прогнозное значение температуры: {:.2f}".format(S[1]*Ж+S[0]))<jupyter_output><empty_output><jupyter_text>Оценим среднее квадратичное отклонение.<jupyter_code>vErr=0
for i in range(len(xs)):
    for j in range(len(ys)):
        if i==j: vErr += (S[1]*xs[j] + S[0] - ys[j])**2
print("Среднеквадратичное отклонение: {:.2f}".format(math.sqrt(vErr)))<jupyter_output>Среднеквадратичное отклонение: 4.09<jupyter_text>Добавим дополнительный признак - квадрат уже имеющегося признака $x$. Для этого построим приближение в виде квадратичной функции $y(x)=Ax^2+Bx+C$ При этом функция ошибки имеет вид: $$F(A,B,C)=\sqrt{\sum_{i=1}^{n}{\left[ Ax_i^2 + Bx_i + C -y_i\right]^2}}\to \min$$. Также как и для предыдущего примера решение задачи может быть найдено аналитически (в виде явной формулы), однако, как указано выше, этот подход не всегда эффективен. Для иллюстрации применения численных методов минимум функции найдем методом градиентного спуска. Градиентные методы основаны на том свойстве, что вектор градиента         $$\nabla F(x,y) = \left\{ F_x, F_y \right\}$$ направлен в сторону максимально быстрого роста функции. Алгоритм работы градиентного спуска:    Выбираем точку $(x_0,y_0)$ - заданное начальное приближение    Делаем шаг с заданной длиной $h$ (например $h=0.01$) в направлении вектора антиградиента (поскольку разыскиваем минимум)    Оцениваем расстояние между соседними точками - если оно меньше заданной точности $\epsilon$         (например, $\epsilon=0.01$), то останавливаем работу алгоритма и последнее найденное значение         $(x_n, y_n)$ объявляем решением. Реализация градиентного спуска в Python$$x_{n+1}=x_n - h \left. \frac{\frac{\partial F(x,y)}{\partial x}}{\left| \nabla \vec{F(x,y)}\right|} \right|_{(x_n,y_n)}y_{n+1}=y_n - h \left. \frac{\frac{\partial F(x,y)}{\partial y}}{\left| \nabla \vec{F(x,y)}\right|} \right|_{(x_n,y_n)}$$<jupyter_code>def f(A,B,C):
    return math.sqrt(sum([(A*(i**2)+B*i+C-j)**2 for i,j in zip(xs,ys)]))

def grad_step(pf, pA, pB, pC, h):
    d=0.001
    l = math.sqrt(
        (pf(pA+d, pB  , pC  ) - pf(pA,pB,pC))**2 +
        (pf(pA  , pB+d, pC  ) - pf(pA,pB,pC))**2 +
        (pf(pA  , pB  , pC+d) - pf(pA,pB,pC))**2
    )
    rA = pA - h*((pf(pA+d, pB  , pC  ) - pf(pA,pB,pC))/l)
    rB = pB - h*((pf(pA  , pB+d, pC  ) - pf(pA,pB,pC))/l)
    rC = pC - h*((pf(pA  , pB  , pC+d) - pf(pA,pB,pC))/l)             
    return round(rA,3), round(rB,3), round(rC,3)<jupyter_output><empty_output><jupyter_text>Применим метод градиентного спуска для отыскания коэффициентов $A$, $B$, $C$ в цикле из 10000 шагов, отобразим график функции $y(x)=Ax^2+Bx+C$ для каждого 1000-го вычисления, очевидно графики приближаются к заданным точкам.<jupyter_code>%%time
m = []
rA, rB, rC = 0, 0, 0
for _ in range(2000):
    rA, rB, rC = grad_step(f, rA, rB, rC, 0.01)
    m.append((rA, rB, rC))
plt.plot(xs,ys,'ro')
for A, B, C in m[50::400]:
    plt.plot(xs, [A*i**2+B*i+C for i in xs])
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Отобразим итоговый результат и используем построенную функцию для предсказания значения при $x=5$.<jupyter_code>Ж = 15
plt.plot(xs, ys, 'ro')
plt.plot(xs, [rA*(i**2)+rB*i+rC for i in xs], marker="o")
plt.plot(Ж, rA*(Ж**2)+rB*Ж+rC, marker='d')
plt.grid(True)
plt.show()
print("Прогнозное значение температуры: {:.2f}".format(rA*(5**2)+rB*5+rC))<jupyter_output><empty_output><jupyter_text>Вычислим функцию ошибки.<jupyter_code>vErr=0
for i in range(len(xs)):
    for j in range(len(ys)):
        if i==j: vErr+=(rA*(xs[j]**2)+rB*xs[j]+rC-ys[j])**2
print("Среднеквадратичное отклонение равно: {:.2f}".format(math.sqrt(vErr)))<jupyter_output>Среднеквадратичное отклонение равно: 4.25<jupyter_text>Для использования модели линейной регрессии в библиотеке sklearn есть специальная модель LinearRegression().<jupyter_code>xs = df.Время.values
x_train = np.array([[i] for i in xs])
ys = df.Температура.values
linreg = LinearRegression()
linreg.fit(x_train, ys)
print(linreg.coef_, linreg.intercept_)
plt.plot(xs, ys,'ro')
plt.plot(xs, [linreg.coef_*i+linreg.intercept_ for i in xs])
plt.xlabel("Время")
plt.ylabel("Температура")
plt.title("Предсказание")
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Для получения приближения с использованием произведений признака $x$ (3,4,5 степеней) можно использовать функцию библиотеки sklearn.preprocessing.PolynomialFeatures. При этом можно убрать лишние признаки с помощью sklearn.linear_model.Ridge. Для объединения нескольких действий в один поток (пайплайн, pipeline) используется функция sklearn.pipeline.make_pipeline.<jupyter_code>xs_test = np.array([1, 2, 3, 4, 5]).reshape(-1,1)
r = PolynomialFeatures(3).fit_transform(xs_test)
print(r)
for count, degree in enumerate([3, 4, 5]):
    model = make_pipeline(PolynomialFeatures(degree), Ridge())
    model.fit(xs.reshape(-1,1), ys)
    y_plot = model.predict(np.array(xs_test).reshape(-1,1))
    plt.scatter(xs_test, y_plot)

plt.plot(xs, ys)    
plt.show()
xs = np.array(range(1,11)).reshape(-1,1)
ys1 = np.array([i for i in xs])
ys2 = np.array([1/i for i in xs])
ys3 = np.array([5/i for i in ys2])

linreg1 = LinearRegression()
linreg1.fit(xs, ys1)
print(linreg1.score(xs, ys1))

linreg2 = LinearRegression()
linreg2.fit(xs, ys2)
print(linreg2.score(xs, ys2))

linreg3 = LinearRegression()
linreg3.fit(xs, ys3)
print(linreg3.score(xs, ys3))
fig, ax = plt.subplots(1,2)
ax[0].scatter(xs, ys1)
ax[0].scatter(xs, linreg1.predict(np.array(xs).reshape(-1,1))+0.1)
ax[1].scatter(xs, ys2)
ax[1].scatter(xs, linreg2.predict(np.array(xs).reshape(-1,1)))

plt.show()<jupyter_output><empty_output><jupyter_text>Назад к Оглавлению  Пример 2. Анализ стоимости квартир<jupyter_code>from sklearn.datasets import load_boston
data = load_boston()
dir(data)<jupyter_output><empty_output><jupyter_text>Данный датасет содержит информацию о стоимости квартир в Бостоне с учетом признаков (перечент признаков 'feature_names', описание можно посмотреть в 'DESCR'). Всего имеются данные о 506 квартирах.<jupyter_code>print(data.DESCR)
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

model = LinearRegression()
X_train, X_test, y_train, y_test = train_test_split(data.data, data.target, test_size=0.3, shuffle=True)
model.fit(X_train, y_train)
model.score(X_test, y_test)
model.predict(data.data[12].reshape(1,-1)), data.target[12]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.152733875849082, 'fraction_numerical': 0.0192346492702939, 'mean_word_length': 5.71028880866426, 'pattern_counts': {'":': 2, '<': 67, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '458656', 'license': 'CC0-1.0', 'n_tokens_mistral': 6617, 'n_tokens_neox': 7241, 'n_words': 1432, 'path': 'lessons/lec-DA-LinearRegression.ipynb', 'repo_name': 'valeriaustenko/Data-Science-Course-USURT', 'size': '14898'}