Математиката зад линейна и логистична регресия¶
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import mglearn
from IPython.display import display
import warnings
warnings.filterwarnings(action="ignore", module="scipy", message="^internal gelsd")
warnings.filterwarnings(action="ignore", module="sklearn", message="^Objective did not")
%matplotlib inline
Ще се опитаме да формулиране проблема на линейната регресия по-формално. Ще вкараме някаква доза математика за да са по-точни нещата. След това ще се направим имплементация на базата на тази формула.
Нотацията леко ще се различава от тази в scikit-learn за да е консистентна с материала по темата (например курса на Andrew Ng).
Първо ще разгледаме линейна регресия с една променлива. Веднъж разбрали това, нещата се генерализират лесно.
Нотация¶
Ще ползваме нотацията на Andrew Ng.
pd.DataFrame([[1, 2, 3, 6], [5, 6, 7, 18], [8, 9, 10, 27], [11, 12, 13, 36]], columns=['$x_1$', '$x_2$', '$x_3$', '$y$'], index=[1, 2, 3, 4])
- $m$ е броя елементи в training set-а
- $n$ е броя feature-и
- $X$ е $m \times n$ матрица с стойностите на feature-ите за всеки ред
- $y$ е вектор с label-ите
- $x^{(i)}$ е $i$-тия ред от feature-ите. Обърнете внимание, че $(i)$ не е степен.
- $x^{(i)}_j$ е $j$-тия feature ($x_j$) на $i$-тия ред. Отново, $(i)$ не е степен.
- $y^{(i)}$ е $i$-тия target
Колко са $m$, $n$, $x^{(2)}$, $x^{(3)}_1$ и $y^{(1)}$?
$m = 4$, $n = 3$, $x^{(2)} = (5, 6, 7)$, $x^{(3)}_1 = 8$, $y^{(1)} = 6$
Формулировка на проблема¶
Имаме някакви данни, за които ще търсим линейна функция, която ги апроксимира.
np.random.seed(0)
m = 50
points = np.linspace(0, 10, m).reshape(-1, 1) + np.random.randn(m, 2) * 0.5
X = points[:, 0].reshape(-1, 1)
y = points[:, 1]
plt.scatter(points[:, 0], points[:, 1]);
pd.DataFrame(points, columns=['$x_1$', '$y$'], index=range(1, m + 1)).head(20)
Ще наричаме функцията "хипотеза":
$$h_\theta(x) = \theta_0 + \theta_1 x_1$$
Тук:
- $\theta$ е вектор от коефициенти
- $h_\theta$ е хипотезата с определените в $\theta$ коефициенти
- $\theta_i$ е $i$-тия компонент на вектора $\theta$ (още $i$-тия коефициент)
- $\theta_0$ е свободния член във хипотезата
Друга нотация е $\hat y$, което е просто кратък запис на хипотезата са даден ред:
$$\hat y = h_\theta(x)$$
И разбира се:
$$\hat{y}^{(i)} = h_\theta(x^{(i)})$$
Ще е полезно за всяка хипотеза да имаме функция, която оценява колко добре се справя. Тази функция е добре да бъде неотрицателна, като по-малки стойности значат по-добро справяне. Ще я наричаме cost функция. За линейна регресия може да ползваме средната квадратна грешка:
$$ J(\theta) = \frac{1}{2m}\sum_{i=1}^m \Big(\hat{y}^{(i)} - y^{(i)}\Big) ^ 2 = \frac{1}{2m}\sum_{i=1}^m \Big(h_\theta(x_1^{(i)}) - y^{(i)}\Big) ^ 2 $$
Така линейната регресия се свежда до много по-прост и разбираем проблем – да намерим глобален минимум на $J(\theta)$, който да ползваме.
Обърнете внимание, че:
- Функцията е винаги неотрицателна, т.е. оптималния възможен минимум за някакъв dataset е 0.
- Функцията е изпъкнала (всъщност, квадратна, макар и в две измерения), което ще рече, че има един минимум.
- 0 не е гарантиран минимум за всеки dataset. В интерес на истината, може да стигнем до 0 само ако всички точки лежат на една права.
Как да намерим минимум, обаче?
Gradient Descent¶
Това е алгоритъм за намиране на минимум на функция. Не е нужно да знаем каква е функцията символично, а само да можем да я изчислим. Първоначално нека го разгледаме извън контекста на линейната регресия.
Ще трябва да си припомним малко училищен анализ.
def draw_derivative(func, x, bounds):
(a, b) = bounds
epsilon = (b - a) * 0.0001
y = func(x)
dx = (func(x + epsilon) - y) / epsilon
di = y - dx * x
xs = np.linspace(a, b, 400)
ys = func(xs)
span = ys.max() - ys.min()
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(xs, ys)
plt.plot([x], [func(x)], marker='o', color='r')
plt.plot(xs, (xs * dx + di), color='r')
plt.grid()
ax.axhline(color='k')
ax.axvline(color='k')
ax.set_ylim(ys.min() - 0.1 * span, ys.max() + 0.1 * span)
Ще разгледаме $x^2 + 0.2$ и нейната производна в някои точки. Тази функция също е изпъкнала.
Какъв е знака на производната в тази точка ($x = 1$)?
draw_derivative(lambda x: x**2 + 0.2, 1.0, (-2, 2))
Знакът е положителен, защото функцията расте.
Каква е стойността на променливата в тази точка ($x = 0.5$), сравнение с предишната?
draw_derivative(lambda x: x**2 + 0.2, 0.5, (-2, 2))
Стойността е по-малка, понеже функцията расте по-бавно.
Какъв е знакът тук ($x = -1$)?
draw_derivative(lambda x: x**2 + 0.2, -1.0, (-2, 2))
Знакът е отрицателен, защото функцията намалява.
А тук ($x = 0$)?
draw_derivative(lambda x: x**2 + 0.2, 0.0, (-2, 2))
Производната е 0, защото има локален минимум.
И финално, каква е връзката между стойността на производната в тази точка и червената права?
draw_derivative(lambda x: x**2 + 0.2, 1.0, (-2, 2))
Стойността на производната е тангенса на ъгъла между червената права и абцисата.
С други думи, производната в дадена точка ни дава информация в каква посока трябва да се предвижим, ако искаме се приближим към най-близкия локален минимум.
Просто трябва да повторим следната операция (с подходящо избрано число $\alpha$) определен брой пъти $N$:
$$ x := x - \alpha f'(x) $$
За начало може да инициализираме $x$ с произволна точка. При достатъчен брой итерации (и достатъчно малка $\alpha$) ще стигнем до минимума.
from itertools import tee
def pairwise(iterable):
a, b = tee(iterable)
next(b, None)
return zip(a, b)
def draw_gradient_descent(func, start, alpha, iterations, region):
(a, b) = region
epsilon = 0.00001
steps = [start]
x = start
for _ in range(iterations):
dx = (func(x + epsilon) - func(x)) / epsilon
x -= alpha * dx
steps.append(x)
xs = np.linspace(a, b, 400)
ys = func(xs)
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(xs, ys)
plt.grid()
ax.axhline(color='k')
ax.axvline(color='k')
for (xa, xb) in pairwise(steps):
ya = func(xa)
yb = func(xb)
plt.annotate('', xy=(xa, ya), xycoords='data',
xytext=(xb, yb), textcoords='data', size=25,
arrowprops={'arrowstyle': '<-', 'ec': 'r', 'lw': 2})
Ако изберем $\alpha$ и $N$ подходящо, с градиентно спускане може да намерим минимума:
draw_gradient_descent(lambda x: x**2 + 0.2, start=4.0, alpha=0.3, iterations=4, region=(-2.0, 4.1))
Ако $\alpha$ е твърде малко, вероятно няма да ни стигнат итерациите:
draw_gradient_descent(lambda x: x**2 + 0.2, start=4.0, alpha=0.03, iterations=10, region=(-2.0, 4.1))
Ако пък $\alpha$ е твърде голямо, резултатите ни няма да клонят към минимума и градиентното спускане ще го прескача постоянно.
draw_gradient_descent(lambda x: x**2 + 0.2, start=4.0, alpha=1.22, iterations=3, region=(-12.0, 12.1))
Може да имплементираме прост алгоритъм, който взема производната на функция и работи с нея:
f = lambda x: x**2 + 2 * x + 1
df = lambda x: 2*x + 2
def gradient_descent(derivative, start, alpha, iterations):
x = start
for _ in range(iterations):
x -= alpha * derivative(x)
return x
x_min = gradient_descent(df, start=-10.0, alpha=0.1, iterations=1000)
x_min
f(x_min)
Ако не знаехме производната, можехме да се оправим и по селския начин:
def gradient_descent(func, start, alpha, iterations):
epsilon = 0.0000001
x = start
for _ in range(iterations):
dx = (func(x + epsilon) - func(x)) / epsilon
x -= alpha * dx
return x
x_min = gradient_descent(f, start=-10.0, alpha=0.01, iterations=1000)
x_min
Резултата, разбира се, е по-неточен, защото така производната е по-неточна.
Повече измерения¶
При повече измерения вземаме частните производни по всяко и създаваме вектор от тях. Този вектор ни е градиента, по който се движим:
$$ x_j := x_j - \alpha \frac{\partial}{\partial x_j}f(x) $$
Ако не помните какво е частна производна – третираме всички други променливи, освен $x_j$ като константи и дифиренцираме.
Gradient Descent върху линейна регресия¶
Връщайки се на линейната регресия, трябва да оптимизираме следната функция.
$$ J(\theta) = \frac{1}{2m}\sum_{i=1}^m \Big(h_\theta(x_1^{(i)}) - y^{(i)}\Big) ^ 2 $$
Тя е изпъкнала и неотрицателна. Трябва да й сметнем производната. Преди това, нека я променим малко.
За да може да работим в вектори, нека да въведем константен feature $x_0 = 1$, така че да имаме $x = (x_0, x_1) = (1, x_1)$.
По този начин работим с вектори с еднаква форма:
$$ x = (x_0, x_1) $$ $$ \theta = (\theta_0, \theta_1) $$
Тогава функцията е:
$$ J(\theta) = \frac{1}{2m}\sum_{i=1}^m \Big(h_\theta(x^{(i)}) - y^{(i)}\Big) ^ 2 $$
И:
$$ h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 $$
Трябва да сметнем:
$$ \frac{\partial}{\partial \theta_j} J(\theta) = \frac{\partial}{\partial \theta_j} \frac{1}{2m}\sum_{i=1}^m \Big(h_\theta(x^{(i)}) - y^{(i)}\Big) ^ 2 $$
Да опростим тая сметка, нека да я видим само с един член от сумата:
$$ \begin{align*} \frac{\partial}{\partial \theta_j}J(\theta) &= \frac{\partial}{\partial \theta_j}\frac{1}{2}\big(h_\theta(x) - y)\big)^2 \\ &= 2 \cdot \frac{1}{2}\big(h_\theta(x) - y \big) \cdot \frac{\partial}{\partial \theta_j}\big(h_\theta(x) - y\big) \\ &= \big(h_\theta(x) - y\big) \cdot \frac{\partial}{\partial \theta_j} \Bigg( \sum_{i=0}^n \theta_i x_i - y \Bigg) \\ &= \big(h_\theta(x) - y\big) x_j \end{align*} $$
Частната производна по $\theta_j$ е:
$$ \frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}\sum_{i=1}^m\Big(\big(h_\theta(x^{(i)}) - y^{(i)}\big)x_j^{(i)}\Big) $$
За градиентно спускане трябва да обновяваме коефициентите така:
$$ \begin{align*} \theta_0 &:= \theta_0 - \alpha \frac{1}{m}\sum_{i=1}^m\big(h_\theta(x^{(i)}) - y^{(i)}\big) \\ \theta_1 &:= \theta_1 - \alpha \frac{1}{m}\sum_{i=1}^m\Big(\big(h_\theta(x^{(i)}) - y^{(i)}\big)x_1^{(i)}\Big) \end{align*} $$
Нека се опитаме да имплементираме тази сметка в NumPy с векторизация. Ще даваме стъпка по стъпка, за да е по-ясно.
Най-лесно ще е да илюстрираме това с dataset.
plt.rc('figure', figsize=(12, 8))
np.random.seed(0)
m = 20
line = np.linspace(0, 10, m)
points = np.random.randn(m, 2) * 0.5 + line.reshape(-1, 1)
X = np.hstack([np.ones((m, 1)), points[:, 0].reshape(-1, 1)])
y = points[:, 1]
plt.scatter(points[:, 0], points[:, 1])
plt.ylim(0, 10);
X
y
Първо ще ползваме scikit-learn да натренираме една линейна регресия, за да видим какъв резултат тъсрим:
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(points[:, :-1], points[:, 1])
best_theta = np.array([model.intercept_, *model.coef_])
best_theta
Ако я начертаем, изглежда убедително:
plt.scatter(points[:, 0], points[:, 1])
plt.plot(line, best_theta[1] * line + best_theta[0], color='b')
plt.ylim(0, 10);
Ние, обаче, ще започнем с по-лоша $\theta$:
theta = np.array([-1., 1.8])
plt.scatter(points[:, 0], points[:, 1])
plt.plot(line, best_theta[1] * line + best_theta[0], color='b')
plt.plot(line, theta[1] * line + theta[0], color='r')
plt.ylim(0, 10);
Сега ще стъпим през една итерация на gradient descent, опитвайки се да го имплементираме с NumPy.
Ще започнем с хипотезата:
$$ \frac{1}{m}\sum_{i=1}^m\Big(\big(\bbox[yellow]{h_\theta(x^{(i)})} - y^{(i)}\big)x_j^{(i)}\Big) $$
Понеже $h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 = \theta^T x$, може да го пренапишем така:
$$ \frac{1}{m}\sum_{i=1}^m\Big(\big(\bbox[yellow]{\theta^T x^{(i)}} - y^{(i)}\big)x_j^{(i)}\Big) $$
Резултата ще бъде NumPy масив с $m$ елемента, всеки от който съдържа предвидената стойност за ред от $X$:
X @ theta
Ще продължим с изваждането на $y$:
$$ \frac{1}{m}\sum_{i=1}^m\Big(\bbox[yellow]{\big(\theta^T x^{(i)} - y^{(i)}\big)}x_j^{(i)}\Big) $$
Резултата вече трябва да бъде грешката за всеки от елементите:
X @ theta - y
Нека добавим и умножението по $x_j$:
$$ \frac{1}{m}\sum_{i=1}^m\Big(\bbox[yellow]{\big(\theta^T x^{(i)} - y^{(i)}\big)x_j^{(i)}}\Big) $$
Тук ще ползваме broadcasting. Резултата ще бъде $n \times m$ матрица, където всеки ред ще съдържа резултатите от това произведение за определен feature. Първия е за $x_0$, втория е за $x_1$:
(X @ theta - y) * X.T
Нека добавим и сумата:
$$ \frac{1}{m}\bbox[yellow]{\sum_{i=1}^m\Big(\big(\theta^T x^{(i)} - y^{(i)}\big)x_j^{(i)}\Big)} $$
Резултата ще бъде масив с градиентите за всеки един от feature-ите. Да не забравяме да сумираме по правилната ос:
((X @ theta - y) * X.T).sum(1)
Финално, трябва да разделим на m:
$$ \bbox[yellow]{\frac{1}{m}\sum_{i=1}^m\Big(\big(\theta^T x^{(i)} - y^{(i)}\big)x_j^{(i)}\Big)} $$
((X @ theta - y) * X.T).sum(1) / m
Може да ползваме този градиент да получим по-добро приближение на най-добрата линейна регресия:
gradient = ((X @ theta - y) * X.T).sum(1) / m
theta1 = theta - 0.02 * gradient
theta1
Като го начертаем, виждаме, че новата линия е по-добра апроксимация
plt.scatter(points[:, 0], points[:, 1])
plt.plot(line, theta[1] * line + theta[0], c='r')
plt.plot(line, theta1[1] * line + theta1[0], c='orange')
plt.ylim(-1, 11);
Ако итерираме още веднъж ще си приближим повече:
theta2 = theta1 - 0.02 * (((X @ theta1 - y) * X.T).sum(1) / m)
theta2
plt.scatter(points[:, 0], points[:, 1])
plt.plot(line, theta[1] * line + theta[0], c='r')
plt.plot(line, theta1[1] * line + theta1[0], c='orange')
plt.plot(line, theta2[1] * line + theta2[0], c='cyan')
plt.ylim(-1, 11);
theta3 = theta2 - 0.02 * (((X @ theta2 - y) * X.T).sum(1) / m)
plt.scatter(points[:, 0], points[:, 1])
plt.plot(line, theta[1] * line + theta[0], c='r')
plt.plot(line, theta1[1] * line + theta1[0], c='orange')
plt.plot(line, theta2[1] * line + theta2[0], c='cyan')
plt.plot(line, theta3[1] * line + theta3[0], c='green')
plt.ylim(-1, 11);
Все още има голяма разлика между това, което правим ние и това, което прави линейната регресия в scikit-learn.
theta3
best_theta
Нека просто пуснем 10000 итерации и да видим какво ще стане:
result = theta
for _ in range(10000):
gradient = ((X @ result - y) * X.T).sum(1) / m
result -= 0.01 * gradient
result
Това го намери и scikit-learn:
best_theta
Нека пробваме да изчислим и $R^2$ score-а:
$$ \operatorname{R^2} = 1 - \frac{\sum_{i=1}^n \big( \hat{y}_i - y_i \big)^2}{\sum_{i=1}^n \big( \bar{y} - y_i \big)^2} $$
Където $\bar y$ е средно-аритметичното на y.
Да започнем с числителя:
$$ \sum_{i=1}^n \big( \bbox[yellow]{\hat{y}_i} - y_i \big)^2 $$
$\hat y$ е просто предвижданията на нашата функция, което е $X\theta^T$. Резултата трябва да е масив от $m$ елемента със предположението за всеки ред от $X$:
X @ theta
Да започнем с числителя:
$$ \sum_{i=1}^n \big( \bbox[yellow]{\hat{y}_i - y_i} \big)^2 $$
Нека извадим $y$ за да получим грешката:
X @ theta - y
Трябва да повдигнем резултата на квадрат:
$$ \sum_{i=1}^n \bbox[yellow]{\big( \hat{y}_i - y_i \big)^2} $$
(X @ theta - y) ** 2
Финално, трябва да го сумираме, за да получим скалар:
$$ \operatorname{R^2} = 1 - \frac{\bbox[yellow]{\sum_{i=1}^n \big( \hat{y}_i - y_i \big)^2}}{\sum_{i=1}^n \big( \bar{y} - y_i \big)^2} $$
((X @ theta - y) ** 2).sum()
Знаменателя е тривиален и ще го направим наведнъж:
$$ \operatorname{R^2} = 1 - \frac{\sum_{i=1}^n \big( \hat{y}_i - y_i \big)^2}{\bbox[yellow]{\sum_{i=1}^n \big( \bar{y} - y_i \big)^2}} $$
((y.mean() - y) ** 2).sum()
Фималното изчисление е ясно:
$$ \operatorname{R^2} = 1 - \frac{\sum_{i=1}^n \big( \hat{y}_i - y_i \big)^2}{\sum_{i=1}^n \big( \bar{y} - y_i \big)^2} $$
1 - ((X @ theta - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()
Какво ли намери нашия scikit-learn по-рано?
model.score(points[:, :-1], points[:, 1])
Същото! Колко яко!
Estimator¶
Май събрахме достатъчно код да си направим собствен estimator.
Ще го даваме на YOLO и направо ще се опитаме да си имплементираме един:
from sklearn.base import BaseEstimator
class PoorMansLinearRegression(BaseEstimator):
def __init__(self, alpha=0.01, iterations=2000):
self.alpha = alpha
self.iterations = iterations
def fit(self, X, y):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
theta = np.random.rand(n + 1)
for _ in range(self.iterations):
gradient = ((X @ theta - y) * X.T).sum(1) / m
theta -= self.alpha * gradient
self.theta_ = theta
return self
def predict(self, X):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
return X @ self.theta_
def score(self, X, y):
X = np.hstack([np.ones((X.shape[0], 1)), X])
return 1 - ((X @ self.theta_ - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()
Отвъд малки дребни детайли, тук няма нищо интересно. Нека си вземем dataset-а с който почнахме, но да го направим малко по-подходящ за scikit-learn (без bias колоната):
X = points[:, :-1]
y = points[:, 1]
X
y
Да тренираме нашия регресор!
poor = PoorMansLinearRegression(alpha=0.02, iterations=10000)
poor.fit(X, y)
print(poor.predict(X))
poor.score(X, y)
Какво ли дава тоя на scikit-learn?
model = LinearRegression().fit(X, y)
print(model.predict(X))
model.score(X, y)
Бижу!
Я да го пробваме с "по-истински dataset":
from sklearn.datasets import load_boston
from sklearn.preprocessing import minmax_scale
from sklearn.model_selection import train_test_split
boston = load_boston()
X = minmax_scale(boston.data)
y = boston.target
X_train, X_test, y_train, y_test = train_test_split(X, y)
proper = LinearRegression().fit(X_train, y_train)
proper.fit(X_train, y_train)
proper.score(X_test, y_test)
Нека сега си натренираме с нашия estimator:
poor = PoorMansLinearRegression(alpha=0.5, iterations=50000)
poor.fit(X_train, y_train)
print("Our score: {}".format(poor.score(X_test, y_test)))
print("sklearn's score: {}".format(proper.score(X_test, y_test)))
print("Bigest error wrt sklearn: {0:.14f}".format(np.abs(proper.predict(X_test) - poor.predict(X_test)).max()))
Normal equation method¶
Всъщност, има и друг подход с който може да решим проблема на линейната регресия – може просто да решим уравнението:
$$\frac{\partial}{\partial \theta_j}J(\theta) = 0 \space \space \space \space \text{(for every j)}$$
Ако приложим малко анализ, може да намерим символично решение. Ще ви спестим математиката, понеже включва добро количество векторен анализ, и ще ви покажем направо резултата:
$$\theta = (X^T X)^{-1}X^Ty$$
Обърнете внимане, че тази операция включва обръщане на матрица, което създава няколко проблема.
Първо, нека си възстановим dataset-а от по-рано с точките, за да го пробваме.
np.random.seed(0)
m = 20
line = np.linspace(0, 10, m)
points = np.random.randn(m, 2) * 0.5 + line.reshape(-1, 1)
plt.scatter(points[:, 0], points[:, 1])
plt.ylim(0, 10);
Нека си направим X и y:
X = np.hstack([np.ones((m, 1)), points[:, :-1]])
y = points[:, 1]
X
y
Сега нека просто сметнем уравнението:
$$\theta = (X^T X)^{-1}X^Ty$$
theta = np.linalg.pinv(X.T @ X) @ X.T @ y
theta
LinearRegression намира същите коефициенти:
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(points[:, :-1], points[:, 1])
np.array([model.intercept_, *model.coef_])
Въпреки това, този метод има няколко проблема:
- Намирането на обратна матрица е бавна операция – $O(n^3)$. При твърде голяма матрица $X$ няма шанс да я изчакате.
- Ако имате линейни зависимости между feature-ите, матрицата ще е необратима.
Регуляризация¶
Как изглежда регуляризацията при линейната регресия? Стандартната (l2) регуляризация просто ползва малко по-различна cost функция с допълнителен член:
$$ J(\theta) = \frac{1}{2m} \Bigg[ \sum_{i=1}^m \Big(h_\theta(x^{(i)}) - y^{(i)}\Big) ^ 2 + \lambda \sum_{j=1}^n\theta_j^2 \Bigg] $$
По този начин се опитваме да държим коефициентите на линейната регресия малки. $\lambda$ е регуляризационния параметър (C или alpha в scikit-learn).
Частната производна и градиента изглеждат малко по-различно:
$$ \theta_j := \theta_j - \alpha \Bigg[ \Bigg (\frac{1}{m}\sum_{i=1}^m\Big(\big(h_\theta(x^{(i)}) - y^{(i)}\big)x_j^{(i)}\Big) \Bigg) + \frac{\lambda}{m}\theta_j \Bigg] $$
Нека натренираме един Ridge, за да има с какво да сравняваме:
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
X = minmax_scale(boston.data, (-1, 1))
y = boston.target
grid = {'alpha': [0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0]}
search = GridSearchCV(Ridge(), grid)
search.fit(X, y)
search.best_estimator_
ridge = Ridge(alpha=10).fit(X, y)
ridge.score(X, y)
Сега ще си направим наша регуляризация:
from sklearn.base import RegressorMixin
class Regularized(BaseEstimator, RegressorMixin):
def __init__(self, alpha=0.01, l=0, iterations=20000):
self.alpha = alpha
self.iterations = iterations
self.l = l
def fit(self, X, y):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
theta = np.random.rand(n + 1)
for _ in range(self.iterations):
regularization = self.l / m * theta
regularization[0] = 0
gradient = ((X @ theta - y) * X.T).sum(1) / m + regularization
theta -= self.alpha * gradient
self.theta_ = theta
return self
def predict(self, X):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
return X @ self.theta_
Може да ползваме този регресор с GridSearchCV:
grid = {'l': [0.1, 1.0, 10.0, 100.0]}
search = GridSearchCV(Regularized(alpha=0.1, iterations=20000), grid)
search.fit(X, y)
search.best_estimator_
Сравнявайки резултатите с Ridge, получаваме нещо сходно:
poor = Regularized(alpha=0.1, iterations=20000, l=10.0).fit(X, y)
poor.score(X, y)
ridge.score(X, y)
Логистична регресия¶
Логистичната регресия строи ползва хипотеза като линейната, но прекаране през сигмоид:
$$ h_\theta(x) = g(\theta^Tx) = g(\theta_0 + \theta_1 x_1 + \theta_2 x_2 + \ldots + \theta_n x_n) $$
Където:
$$ g(z) = \frac{1}{1 + e^{-z}} $$
Ще си дефинираме функцията за сигмоида:
def sigmoid(z):
return 1 / (1 + np.exp(-z))
_x = np.linspace(-10, 10)
plt.plot(_x, [sigmoid(i) for i in _x])
plt.grid()
За cost функция използваме:
$$ \begin{align*} & J(\theta) = \frac{1}{m} \sum_{i=1}^{m}\mathrm{Cost} \big( h_\theta(x^{(i)}), y^{(i)}) \big) \\ & \mathrm{Cost}(h_\theta(x), y) = -log(1 - h_\theta(x)) & \text{if y = 0} \\ & \mathrm{Cost}(h_\theta(x), y) = -log(h_\theta(x)) & \text{if y = 1} \end{align*} $$
Може да я запишем и така:
$$ J(\theta) = - \frac{1}{m} \sum_{i=1}^{m} \Big[ y^{(i)}\log(h_\theta(x^{(i)})) + (1 - y^{(i)})\log(1 - h_\theta(x^{(i)})) \Big] $$
Или векторно:
$$ \begin{align*} & h = g(X\theta) \\ & J(\theta) = \frac{1}{m} \cdot \Big( -y^T \log(h) - (1 - y)^T\log(1 - h) \Big) \end{align*} $$
Градиента е:
$$ \theta_j := \theta_j - \frac{\alpha}{m}\sum_{i=1}^{m}\Big(h_\theta(x^{(i)}) - y^{(i)}\Big)x^{(i)}_j $$
Нека се опитаме да имплементираме и това. Първо ще започнем с някакъв dataset:
from sklearn.datasets import make_classification
X, y = make_classification(n_samples=20, n_features=2, n_redundant=0, random_state=42)
mglearn.discrete_scatter(X[:, 0], X[:, 1], y);
X = np.hstack([np.ones((X.shape[0], 1)), X])
(m, n) = X.shape
X
Ще започнем с някаква хипотеза:
theta = np.random.rand(n)
Градиента ще се опитаме да го напишем наведнъж:
$$ \theta_j := \theta_j - \frac{\alpha}{m}\sum_{i=1}^{m}\Big(h_\theta(x^{(i)}) - y^{(i)}\Big)x^{(i)}_j $$
((sigmoid(X @ theta) - y).reshape(-1, 1) * X).sum(0)
Да пуснем няколко итерации:
theta = np.random.rand(n)
for _ in range(50000):
theta -= 0.3 * ((sigmoid(X @ theta) - y).reshape(-1, 1) * X).sum(0) / m
theta
Да видим какво ще направи логистичната регресия.
from sklearn.linear_model import LogisticRegression
log = LogisticRegression(C=1000000000000).fit(X[:, 1:], y)
np.array([log.intercept_[0], *log.coef_[0]])
Тук резултатите се различават малко повече, защото има регуляризация. Голямо C значи малка регуляризация (но все пак, не и 0).
Ще си напишем логистична регресия с регуляризация. Ето го градиента:
$$ \theta_j := \theta_j - \frac{\alpha}{m} \Bigg[ \sum_{i=1}^{m}\Big(h_\theta(x^{(i)}) - y^{(i)}\Big)x^{(i)}_j + \lambda\theta_j \Bigg] $$
Абсолютно същото като линейната.
from sklearn.base import ClassifierMixin
class ItsSomethingLogistic(BaseEstimator, ClassifierMixin):
def __init__(self, alpha=0.01, l=1.0, iterations=20000):
self.alpha = alpha
self.iterations = iterations
self.l = l
def fit(self, X, y):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
theta = np.random.rand(n + 1)
for _ in range(self.iterations):
regularization = self.l / m * theta
regularization[0] = 0
gradient = ((sigmoid(X @ theta) - y).reshape(-1, 1) * X).sum(0) / m + regularization
theta -= self.alpha * gradient
self.theta_ = theta
return self
def predict(self, X):
(m, n) = X.shape
X = np.hstack([np.ones((m, 1)), X])
return np.round(sigmoid(X @ self.theta_))
Може да тестваме нашия класификатор на cancer dataset-а:
from sklearn.datasets import load_breast_cancer
cancer = load_breast_cancer()
X = minmax_scale(cancer.data, (-1, 1))
y = cancer.target
Ще ни трябват повече итерации, смела $\alpha$ и никаква регуляризация:
np.random.seed(0)
ours = ItsSomethingLogistic(l=0.0, alpha=20.0, iterations=500000)
ours.fit(X, y)
ours.score(X, y)
Да видим какъв резултат вади scikit-learn:
LogisticRegression(C=10000).fit(X, y).score(X, y)
Чудесно! Вадим същото, макар и много по-бавно.