Нелинейная регрессия#
Градиентный спуск#
Напомним формулировку задачу регрессии (семинар 1.3).
Пусть имеем конечный набор точек \(\left\{\left(x_i, y_i\right)\right\}_{i=1}^N\), причём \(x_i \in \mathbb{R}^m, y_i \in \mathbb{R}\). Хотим найти функцию \(\check{y}\), наилучшим образом аппроксимирующий данный набор точек.
«Наилучшим» в плане минимизации какой-либо функции ошибки (Loss). Самой простой и широко-используемой является средний квадрат отклонений (MSE)
В отличие от предыдущих семинаров, сейчас рассмотрим задачу нелинейной регрессии.
Пусть ищем наилучшую функцию среди некоторого класса функций \(\check{y} = f(w, x)\) - где зависимость от вектора весов \(w=(w_1, w_2, ..., w_m)\) в общем случае нелинейна.
Алгоритм не меняется. Пытаемся минимизировать MSE по весам \(w_k\):
где введено обозначение для невязки на каждой точке из выборки \(r_i = y_i - f(w, x_i)\). Как обычно, невязка это «то что хотим» минус «что имеем».
Это выражение можно минимизировать итеративно с помощью любого локального минимизатора выше (или его стохастической версии). Вот, например, обычный градиентный спуск:
\[ w^{(n+1)}_k=w^{(n)}_k-\alpha \cdot \left.\frac{\partial \mathcal{L}(w)}{\partial w_{k}}\right|_{w^{(n)}} =w^{(n)}_k+\frac{2\alpha}{N}\cdot\sum_{i=1}^N r_i\cdot \left. \frac{\partial f(w, x_i)}{\partial w_k} \right|_{w^{(n)}} \]
Также можно воспользоваться итерационным методом нормального уравнения.
Метод нормального уравнения#
В предположении о том, что на каждом шаге мы будем менять веса несильно, разложим модельную функцию в ряд Тейлора по приращению весов на каждом шаге и на каждом \(x_i\):
где введено обозначение \(J_{ik}^{(n)}=\left. \frac{\partial f(w, x_i)}{\partial w_k} \right|_{w^{(n)}}\).
Заметим,
А тогда производная Loss-а:
Далее забьём на константу перед суммой. Положим, что для итерации \(w^{(n+1)}\) данное равенство выполнено (или по крайней мере выполнено приближенно) и разложим \(f(w^{(n+1)},x_i)\) в ряд Тейлора:
Перетасуем слагаемые:
Запишем это равенство в матричном виде (\(r^{(n)}\)- столбец невязок для каждого объекта при \(w^{(n)}\)):
\[ \left(J_{(n)}^TJ_{(n)}\right)\Delta w^{(n+1)} = J_{(n)}^T\cdot r^{(n)} \]
Эта система линейных уравнений называется нормальным уравнением. Решая её, мы получаем наилучшую траекторию для поиска минимума Lossa. Заметим, что в случае линейной регрессии, сходимся к ответу за 1 итерацию.
Пример с нелинейной регрессией окружности#
Напишем функцию, которая по заданным координатам 4 точек на плоскости находит уравнение окружности (координаты центра x0, y0 и радиус R), для которой сумма квадратов расстояний от этих 4 точек до центра окружности является минимальной. Нарисуем график (отметим 4 точки, центр и нарисуем саму окружность).
from scipy.optimize import minimize # в данном методе реализовано много методов локальной минимизации
import numpy as np
def distance_to_circle(t, x0, y0, R, x, y):
# x0, y0, R - параметры окружности
# x, y - точка, до которой определяем минимальное расстояние от окружности
# t - угол на окружности, до точки до которой определяем расстояние
x_na_okr = R*np.cos(t) + x0
y_na_okr = R*np.sin(t) + y0
return ((x_na_okr - x)**2 + (y_na_okr - y)**2)**2
def minimum_distance_to_circle(x0, y0, R, x, y):
# Определяет минимальное расстояние от точки до окружности
dist_min = (minimize(distance_to_circle, 0, (x0, y0, R, x, y)))['fun']
return dist_min
# Реализуем функцию ошибки
def MSE(circle, *args):
x0 = circle[0]
y0 = circle[1]
R = circle[2]
x1, y1, x2, y2, x3, y3, x4, y4 = args
return 1./4*(minimum_distance_to_circle(x0, y0, R, x1, y1)**2
+ minimum_distance_to_circle(x0, y0, R, x2, y2)**2
+ minimum_distance_to_circle(x0, y0, R, x3, y3)**2
+ minimum_distance_to_circle(x0, y0, R, x4, y4)**2)
# Генерируем 4 случайные точки на плоскости
x1 = 10*np.random.random() - 5
y1 = 10*np.random.random() - 5
x2 = 10*np.random.random() - 5
y2 = 10*np.random.random() - 5
x3 = 10*np.random.random() - 5
y3 = 10*np.random.random() - 5
x4 = 10*np.random.random() - 5
y4 = 10*np.random.random() - 5
# Минимизируем MSE по параметрам окружности
x0_start = 0
y0_start = 0
R_start = 1
circle_start = np.array([x0_start, y0_start, R_start])
minimize(MSE, circle_start, (x1, y1, x2, y2, x3, y3, x4, y4))
message: Optimization terminated successfully.
success: True
status: 0
fun: 2.0665597889114427e-07
x: [ 7.465e-01 -1.123e+00 3.941e+00]
nit: 53
jac: [ 6.657e-06 2.268e-06 -8.472e-06]
hess_inv: [[ 1.944e+02 -2.241e+02 -2.884e+01]
[-2.241e+02 2.638e+02 -3.534e+01]
[-2.884e+01 -3.534e+01 8.459e+02]]
nfev: 244
njev: 61
circle = minimize(MSE, circle_start, (x1, y1, x2, y2, x3, y3, x4, y4))['x']
x0 = circle[0]
y0 = circle[1]
R = circle[2]
import matplotlib.pyplot as plt
plt.figure(figsize = (7, 7))
ts = np.linspace(0, 2*np.pi, 200)
xs = x0 + R*np.cos(ts)
ys = y0 + R*np.sin(ts)
plt.plot(xs, ys)
plt.scatter(x0, y0)
plt.scatter(x1, y1, c='r')
plt.scatter(x2, y2, c='r')
plt.scatter(x3, y3, c='r')
plt.scatter(x4, y4, c='r')
<matplotlib.collections.PathCollection at 0x7fb7c3e8e8e0>
Пример регрессии нелинейной функцией на scipy#
Нелинейную регрессию очень просто реализовать, используя метод curve_fit. Как сЛеДуЕт из названия, он подбирает нелинейные параметры функции для удовлетворения близости какому-либо множеству.
Примечание. «Под капотом» он использьует различные методы локальной оптимизации. Методы глобальной оптимизации можно посмотреть здесь.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a * np.exp(-b * x) + c
xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
y_noise = 0.1 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.scatter(xdata, ydata, label='data')
popt, pcov = curve_fit(func, xdata, ydata, p0 = [1, 1, 1])
print(popt)
plt.plot(xdata, func(xdata, *popt), 'r-',
label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
[2.49863354 1.4204748 0.56672885]
[<matplotlib.lines.Line2D at 0x7fb7cea55760>]
