Шестаков А.В. Майнор по анализу данных - 02/02/2016
Сегодня мы рассмотрим следущие темы:
1) Методы оптимизации функций в Python и графика с ней связанная
2) Символьные вычисления
3) Решение задачи линейной регрессии тремя способами!
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np
import scipy.optimize as opt
import sympy
plt.style.use('ggplot')
%matplotlib inline
Рассмотрим две функции:
x = np.linspace(-5, 5, 100)
y1 = lambda(x): 0.5*x**2 + 10*np.sin(x) - 2
y2 = lambda(x): x**2 + 0.5**x - 4
fig, ax = plt.subplots(1, 2, figsize=(14,7))
ax[0].plot(x, y1(x))
ax[0].set_ylabel('$y_1$', fontsize=20)
ax[0].set_xlabel('$x$', fontsize=20)
ax[1].plot(x, y2(x))
ax[1].set_ylabel('$y_2$', fontsize=20)
ax[1].set_xlabel('$x$', fontsize=20)
t = [0,5]
alpha = 0.5
f1 = alpha*y1(t[0]) + (1-alpha)*y1(t[1])
f2 = y1(alpha*t[0] + (1-alpha)*t[1])
ax[0].plot(alpha*t[0] + (1-alpha)*t[1], f1, '*r')
ax[0].plot(alpha*t[0] + (1-alpha)*t[1], f2, 'sb')
Функция $f(x)$ называется выпуклой (convex), если для любых двух точек $u$ и $v$ и для любого числа $\alpha \in [0,1]$ выполняется: $$ f\left(\alpha u + [1-\alpha] v\right) \leq \alpha f(u) + (1-\alpha) f(v) $$
Какая из этих функций является выпуклой? И заодно посмотрим, что это за выражение такое.
Чем хороши выпуклые фукнции с точки хрения поиска оптимального значения?
Посмотрим на выпуклую функцию в 3D!
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x,y)
Z = X**2 + 0.25*Y**2
fig = plt.figure(figsize=(14, 7))
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.view_init(40, 25)
ax.plot_surface(X, Y, Z, alpha=0.3,)
# ax.plot_(X, Y, Z)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax = fig.add_subplot(1, 2, 2)
contour = ax.contour(X, Y, Z)
plt.clabel(contour, inline=1, fontsize=10)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
Гладкие функции (непрерывно дифференцируемые, smooth functions) удобны тем, что в любой точке можно посчитать её производую (градиент)
x = np.linspace(-5, 5, 100)
y1 = lambda(x): 10*x**(2)
y2 = lambda(x): abs(x)
fig, ax = plt.subplots(1, 2, figsize=(14,7))
ax[0].plot(x, y1(x))
ax[0].set_ylabel('$y_1$', fontsize=20)
ax[0].set_xlabel('$x$', fontsize=20)
ax[1].plot(x, y2(x))
ax[0].set_ylabel('$y_2$', fontsize=20)
ax[0].set_xlabel('$x$', fontsize=20)
В зависимости от того, какой функции вы смотрите в глаза, нужно выбирать свои методы. О описанием различных ситуаций и методов.
Для вас есть хорошая новость. Как правило мы будем работать с выпуклыми и гладкими функциями. В общем случае это делается так:
# Задаём функцию
def f(x):
return x**2 + np.exp(x) - 4
# Находим минимум
x_min = opt.minimize(f, -3)
x_min
x = np.linspace(-5, 5, 100)
plt.plot(x, f(x))
plt.plot(x_min['x'], x_min['fun'], 'sb')
Символьные вычисления в Python выполняются в модуле sympy
sympy.init_printing(use_latex=True)
a = sympy.Rational(1, 2)
print a
print a + 2
print sympy.pi
print sympy.pi.evalf(n=4)
x = sympy.Symbol('x')
y, z = sympy.symbols('y z')
w = (x + y)**2
w
sympy.expand(w)
sympy.factor(x**2 - y**2)
sympy.simplify(w - 2*x*y)
sympy.simplify(sympy.sin(x)/sympy.cos(x))
sympy.diff(w, x)
sympy.diff(w, x, 2)
sympy.solve(x**2 - 4, x)
res = sympy.solve([x + 2*y - 2,
x + 3*y - 6], x, y)
res
type(res)
При оптимизации функций, для некоторых методлов необходимо расчитать градиент функций.
Используйте sympy
если испытывайте какие-то сложности, затем, дополнительно проверяйте с scipy.optimize.check_grad()
sympy.init_printing(use_latex=False)
Рассмотрим 3 наблюдения:
X = np.array([[1, 1],
[1, 2],
[1, 3]])
y = np.array([1, 2, 2])
plt.scatter(X[:,1], y)
plt.xlabel('$x_1$')
plt.ylabel('$y$')
Можем ли мы провести через них прямую? Почему?
Как бы вы решали эту задачу, если бы нужно было провести прямую через 2 точки?
Если система $$ A x = b $$ не имеет решения, то решайте $$A^\top A x = A^\top b$$
Beta = np.linalg.inv((X.T.dot(X))).dot(X.T).dot(y)
print Beta
y_hat = X.dot(Beta)
plt.scatter(X[:,1], y)
plt.plot(X[:,1], y_hat)
$X$ - признаковое описание наблюдений,<br> $y$ - прогнозируемая величина
Пусть задана функция ошибки (функция потерь) $L(\cdot)$. <br> Нам надо построить такой функционал $f(X)$, который будет выдавать значение наиболее близкие к $y$, иначе говоря: $$L\left(f(X) - y\right) \rightarrow\min $$
Определим функцию потерь, как сумму квадратов разности выдаваемого ответа функционала и реального значения: $$ L(\cdot) = \frac{1}{2n}\sum_{i=1}^n(f(x^{(i)}) - y^{(i)})^2$$
Так как среди всего множества моделей мы выбрали линейную регрессию, то $$f(X) = \beta_0 + \beta_1x_1$$ Подставляем это выражение в $L(\cdot)$ и находим $\beta_0$, $\beta_1$!
Градиентый спукс - это итеративный метод оптимизации функции. Он заключается в постепенном перемещении к точке экспетмума в направлении градиента этой функции в точке.
Пусть дана функция $L(\theta)$, необходимо найти $\hat{\theta}:\ L(\hat{\theta}) = \min L(\theta)$ <br>
Шаги алгоритма:
1. Случайным образом фиксируется начальное состояние
2. Пока алгоритм не сойдется выполняется обновление:
$$ \theta = \theta - \alpha \frac{\partial L}{\partial \theta} $$где $\alpha$ - скорость спуска
Пример: $$L(\theta) = 3\theta^2 + 2\theta - 10$$ Найдем градиент (в данном случае простую производую, так как у нас одна переменная) $$ \frac{\partial L}{\partial \theta} = 6\theta + 2$$
def f(x):
return 3*x**2 + 2*x - 10
def grad(x):
return 6*x + 2
iter_max = 500 # максимальное кол-во итераций
alpha = 0.01 # скорость спуска
old_min = 0
temp_min = 4
precision = 0.001 # требуемая точность
i = 0
mins = [temp_min] # значений на каждой итерации алгоритма
cost = [] # разность с инстинным решением
# Напишите код, реализующий данный алгоритм
#
while abs(temp_min - old_min) > precision and i<iter_max:
## ...
x = np.linspace(-5, 5, 100)
# Рисуем функцию
fig, ax = plt.subplots(1, 2, figsize=(14,7))
ax[0].plot(x, f(x), 'b', linewidth=2)
ax[0].set_ylabel('$ L(\\theta) $', fontsize=20)
ax[0].set_xlabel('$\\theta$', fontsize=20)
# Рисуем пусть градиента
mins = np.array(mins)
idx = np.linspace(0, len(mins)-1, 10).astype(int)
ax[0].plot(mins[idx], f(mins[idx]), '-*')
# Смотрим как уменьшается ошибка с итерацией
ax[1].plot(cost)
ax[1].set_ylabel('$Error$', fontsize=20)
ax[1].set_xlabel('$Iteration$', fontsize=20)
Вернемся к линейной регрессии..
Теперь, когда мы знаем, градиент функции квадратичной ошибки для линейной регрессии $\beta_0$, $\beta_1$ можно использовать метод градиентного спуска и для неё!
Важно! Веса должны обновляться одновременно, а не по очереди!
# Напишите код, реализующий данный алгоритм для линейной регрессии
#
def gradient_descent(X, y, iters, alpha):
costs = []
n = y.shape[0]
Beta = np.random.rand(2)
history = [Beta]
preds = []
for i in xrange(iters):
#....
return history, costs, preds, Beta