本文讨论了确定简单（配对）回归线的数学方程的几种方法。
这里讨论的所有求解方程的方法都基于最小二乘法。 让我们将这些方法表示如下：
- 解析解
- 梯度下降
- 随机梯度下降
对于每种求解直线方程的方法，文章都提供了各种函数，主要分为不使用库编写的函数 NumPy的 以及那些用于计算的 NumPy的。 相信熟练运用 NumPy的 将降低计算成本。
文章中给出的所有代码都是用该语言编写的 python 2.7 同 Jupyter笔记本。 源代码和包含示例数据的文件发布在
本文更针对初学者和那些已经逐渐开始掌握人工智能非常广泛的部分（机器学习）的研究的人。
为了说明本材料，我们使用一个非常简单的例子。
条件示例
我们有五个表征依赖的价值观 Y 从 X （表 1）：
表 1“示例条件”
我们假设这些值 是一年中的月份，并且 ——本月收入。 换句话说，收入取决于一年中的月份，并且 - 收入所依赖的唯一标志。
这个例子马马虎虎，无论是从收入对一年中月份的条件依赖的角度来看，还是从值的数量来看——它们都很少。 然而，正如他们所说，这种简化将有可能解释初学者所吸收的材料，但并不总是那么容易。 而且数字的简单性将使那些希望在纸上解决示例的人无需花费大量的劳动力。
让我们假设示例中给出的相关性可以通过以下形式的简单（成对）回归线的数学方程很好地近似：
哪里 是收到收入的月份， — 对应月份的收入， и 是估计线的回归系数。
注意系数 通常称为估计线的斜率或坡度； 代表的金额 当它改变时 .
显然，我们在示例中的任务是在方程中选择这样的系数 и ，其中我们按月计算的收入值与真实答案的偏差，即样本中呈现的值将是最小的。
最小二乘法
根据最小二乘法，应通过对其进行平方来计算偏差。 如果偏差具有相反的符号，则此技术可让您避免相互抵消。 例如，如果在一种情况下，偏差是 +5 （加五），以及在其他 -5 （负五），那么偏差之和将相互抵消并等于 0（零）。 可以不对偏差进行平方，而是利用模数的性质，然后所有偏差都将是正值并且会累加。 关于这一点，我们不做详细阐述，只是简单说明，为了计算方便，习惯上对偏差进行平方。
这就是我们用来确定最小偏差平方和（误差）的公式：
哪里 是真实答案的近似函数（即我们计算的收入），
是真实答案（样本中提供的收入），
是样本索引（确定偏差的月份数）
让我们对函数进行微分，定义偏微分方程，并准备好继续求解解析解。 但首先，让我们简单了解一下微分是什么，并记住导数的几何意义。
差异化
微分是求函数导数的运算。
衍生品有什么用？ 函数的导数表征了函数的变化率并告诉我们它的方向。 如果给定点的导数为正，则函数增大；否则，函数减小。 并且绝对导数的值越大，函数值的变化率就越高，函数图形的斜率也就越陡。
例如，在直角坐标系条件下，M(0,0)点处的导数值等于 +25 意味着在给定点，当值发生变化时 向右按常规单位，值 增加25个常规单位。 在图表上，值看起来相当陡峭的上升 从给定点。
另一个例子。 导数值相等 -0,1 意味着当移位时 每一个常规单位，价值 仅减少0,1个常规单位。 同时，在函数图上，我们可以观察到几乎不明显的向下斜率。 与一座山进行类比，就好像我们从一座山上慢慢地走下缓坡，不像前面的例子，我们必须爬上非常陡峭的山峰:)
因此，对函数求导后 按赔率 и ，我们定义一阶偏微分方程。 确定方程后，我们将得到一个由两个方程组成的系统，通过求解该方程组，我们将能够选择系数的值 и ，给定点处相应导数的值变化非常非常小，而在解析解的情况下根本不变化。 换句话说，找到的系数处的误差函数将达到最小值，因为这些点处的偏导数值将等于零。
因此，根据微分规则，关于系数的一阶偏导方程 将采用以下形式：
一阶偏导方程关于 将采用以下形式：
结果，我们得到了一个具有相当简单解析解的方程组：
开始{equation *}
开始{案例}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
结束{案例}
结束{equation *}
在解方程之前，我们先预加载，检查加载是否正确，并格式化数据。
加载和格式化数据
应该注意的是，由于对于解析解以及随后的梯度和随机梯度下降，我们将使用两种变体的代码：使用库 NumPy的 如果不使用它，那么我们将需要适当的数据格式（参见代码）。
数据加载和处理代码
# импортируем все нужные нам библиотеки
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
import pylab as pl
import random
# графики отобразим в Jupyter
%matplotlib inline
# укажем размер графиков
from pylab import rcParams
rcParams['figure.figsize'] = 12, 6
# отключим предупреждения Anaconda
import warnings
warnings.simplefilter('ignore')
# загрузим значения
table_zero = pd.read_csv('data_example.txt', header=0, sep='t')
# посмотрим информацию о таблице и на саму таблицу
print table_zero.info()
print '********************************************'
print table_zero
print '********************************************'
# подготовим данные без использования NumPy
x_us = []
[x_us.append(float(i)) for i in table_zero['x']]
print x_us
print type(x_us)
print '********************************************'
y_us = []
[y_us.append(float(i)) for i in table_zero['y']]
print y_us
print type(y_us)
print '********************************************'
# подготовим данные с использованием NumPy
x_np = table_zero[['x']].values
print x_np
print type(x_np)
print x_np.shape
print '********************************************'
y_np = table_zero[['y']].values
print y_np
print type(y_np)
print y_np.shape
print '********************************************'可视化
现在，在我们首先加载数据，其次检查加载的正确性并最后格式化数据之后，我们将进行第一次可视化。 经常使用的方法是 配对图 文库 海生。 在我们的示例中，由于数量有限，使用该库没有意义 海生。 我们将使用常规库 Matplotlib 只需查看散点图即可。
散点图代码
print 'График №1 "Зависимость выручки от месяца года"'
plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()图1“收入对月份的依赖性”
解析解
让我们使用最常用的工具 蟒蛇 并求解方程组：
开始{equation *}
开始{案例}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
结束{案例}
结束{equation *}
根据克莱默法则 我们将找到一般行列式以及以下行列式 和 ，然后将行列式除以 一般行列式 - 求系数 ，类似地我们求出系数 .
解析解代码
# определим функцию для расчета коэффициентов a и b по правилу Крамера
def Kramer_method (x,y):
# сумма значений (все месяца)
sx = sum(x)
# сумма истинных ответов (выручка за весь период)
sy = sum(y)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x[i]*y[i]) for i in range(len(x))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x[i]**2) for i in range(len(x))]
sx_sq = sum(list_x_sq)
# количество значений
n = len(x)
# общий определитель
det = sx_sq*n - sx*sx
# определитель по a
det_a = sx_sq*sy - sx*sxy
# искомый параметр a
a = (det_a / det)
# определитель по b
det_b = sxy*n - sy*sx
# искомый параметр b
b = (det_b / det)
# контрольные значения (прооверка)
check1 = (n*b + a*sx - sy)
check2 = (b*sx + a*sx_sq - sxy)
return [round(a,4), round(b,4)]
# запустим функцию и запишем правильные ответы
ab_us = Kramer_method(x_us,y_us)
a_us = ab_us[0]
b_us = ab_us[1]
print ' 33[1m' + ' 33[4m' + "Оптимальные значения коэффициентов a и b:" + ' 33[0m'
print 'a =', a_us
print 'b =', b_us
print
# определим функцию для подсчета суммы квадратов ошибок
def errors_sq_Kramer_method(answers,x,y):
list_errors_sq = []
for i in range(len(x)):
err = (answers[0] + answers[1]*x[i] - y[i])**2
list_errors_sq.append(err)
return sum(list_errors_sq)
# запустим функцию и запишем значение ошибки
error_sq = errors_sq_Kramer_method(ab_us,x_us,y_us)
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений" + ' 33[0m'
print error_sq
print
# замерим время расчета
# print ' 33[1m' + ' 33[4m' + "Время выполнения расчета суммы квадратов отклонений:" + ' 33[0m'
# % timeit error_sq = errors_sq_Kramer_method(ab,x_us,y_us)这是我们得到的：
因此，系数的值已经找到，偏差平方和也已经确定。 让我们根据找到的系数在散射直方图上画一条直线。
回归线代码
# определим функцию для формирования массива рассчетных значений выручки
def sales_count(ab,x,y):
line_answers = []
[line_answers.append(ab[0]+ab[1]*x[i]) for i in range(len(x))]
return line_answers
# построим графики
print 'Грфик№2 "Правильные и расчетные ответы"'
plt.plot(x_us,y_us,'o',color='green',markersize=16, label = '$True$ $answers$')
plt.plot(x_us, sales_count(ab_us,x_us,y_us), color='red',lw=4,
label='$Function: a + bx,$ $where$ $a='+str(round(ab_us[0],2))+',$ $b='+str(round(ab_us[1],2))+'$')
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.legend(loc=1, prop={'size': 16})
plt.show()图 2“正确和计算出的答案”
您可以查看每个月的偏差图。 在我们的例子中，我们不会从中获得任何重要的实用价值，但我们将满足我们的好奇心，即简单的线性回归方程如何很好地描述收入对一年中月份的依赖性。
偏差图代码
# определим функцию для формирования массива отклонений в процентах
def error_per_month(ab,x,y):
sales_c = sales_count(ab,x,y)
errors_percent = []
for i in range(len(x)):
errors_percent.append(100*(sales_c[i]-y[i])/y[i])
return errors_percent
# построим график
print 'График№3 "Отклонения по-месячно, %"'
plt.gca().bar(x_us, error_per_month(ab_us,x_us,y_us), color='brown')
plt.xlabel('Months', size=16)
plt.ylabel('Calculation error, %', size=16)
plt.show()第 3 号图表“偏差，%”
并不完美，但我们完成了任务。
让我们编写一个函数来确定系数 и 使用图书馆 NumPy的，更准确地说，我们将编写两个函数：一个使用伪逆矩阵（在实践中不推荐，因为该过程计算复杂且不稳定），另一个使用矩阵方程。
分析解决方案代码 (NumPy)
# для начала добавим столбец с не изменяющимся значением в 1.
# Данный столбец нужен для того, чтобы не обрабатывать отдельно коэффицент a
vector_1 = np.ones((x_np.shape[0],1))
x_np = table_zero[['x']].values # на всякий случай приведем в первичный формат вектор x_np
x_np = np.hstack((vector_1,x_np))
# проверим то, что все сделали правильно
print vector_1[0:3]
print x_np[0:3]
print '***************************************'
print
# напишем функцию, которая определяет значения коэффициентов a и b с использованием псевдообратной матрицы
def pseudoinverse_matrix(X, y):
# задаем явный формат матрицы признаков
X = np.matrix(X)
# определяем транспонированную матрицу
XT = X.T
# определяем квадратную матрицу
XTX = XT*X
# определяем псевдообратную матрицу
inv = np.linalg.pinv(XTX)
# задаем явный формат матрицы ответов
y = np.matrix(y)
# находим вектор весов
return (inv*XT)*y
# запустим функцию
ab_np = pseudoinverse_matrix(x_np, y_np)
print ab_np
print '***************************************'
print
# напишем функцию, которая использует для решения матричное уравнение
def matrix_equation(X,y):
a = np.dot(X.T, X)
b = np.dot(X.T, y)
return np.linalg.solve(a, b)
# запустим функцию
ab_np = matrix_equation(x_np,y_np)
print ab_np让我们比较一下确定系数所花费的时间 и ，按照所提出的 3 种方法。
计算计算时间的代码
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов без использования библиотеки NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием псевдообратной матрицы:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием матричного уравнения:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)
对于少量数据，“自写”函数出现在前面，它使用 Cramer 方法找到系数。
现在您可以继续使用其他方法来查找系数 и .
梯度下降
首先，我们来定义什么是渐变。 简单来说，梯度就是表示函数最大增长方向的一段。 类比爬山，坡度面向的地方就是登顶山顶最陡的地方。 以山为例，我们记得实际上我们需要最陡的下降才能尽快到达低地，即最小值 - 函数不增加或减少的地方。 此时导数将等于零。 因此，我们不需要梯度，而是反梯度。 要找到反梯度，您只需将梯度乘以 -1 （减一）。
让我们注意这样一个事实：一个函数可以有多个最小值，并且使用下面提出的算法下降到其中一个最小值时，我们将无法找到另一个最小值，该最小值可能低于找到的最小值。 放心吧，这对我们没有威胁！ 在我们的例子中，我们正在处理一个最小值，因为我们的函数 图上是一条正抛物线。 正如我们从学校数学课程中应该知道的那样，抛物线只有一个最小值。
当我们发现为什么需要梯度，并且梯度是一个段，即具有给定坐标的向量，它们是完全相同的系数之后 и 我们可以实现梯度下降。
在开始之前，我建议先阅读几句话有关下降算法的内容：
- 我们以伪随机方式确定系数的坐标 и 。 在我们的示例中，我们将定义接近于零的系数。 这是一种常见的做法，但每种情况可能都有自己的做法。
- 从坐标 减去该点的一阶偏导数的值 。 因此，如果导数为正，则函数增大。 因此，通过减去导数的值，我们将向增长的相反方向移动，即向下降的方向移动。 如果导数为负，则此时的函数会减小，通过减去导数的值，我们将向下降方向移动。
- 我们用坐标进行类似的操作 ：减去该点的偏导数值 .
- 为了不跳过最小值而飞入深空，需要设置下降方向的步长。 一般来说，您可以写一篇完整的文章来介绍如何正确设置步长以及如何在下降过程中更改步长以减少计算成本。 但现在我们面前的任务略有不同，我们将使用“戳”的科学方法，或者正如他们通常所说的那样，凭经验确定步长。
- 一旦我们从给定的坐标开始 и 减去导数的值，我们得到新的坐标 и 。 我们已经根据计算出的坐标进行下一步（减法）。 如此循环一次又一次地开始，直到达到所需的收敛。
全部！ 现在我们准备出发寻找马里亚纳海沟最深的峡谷。 让我们开始吧。
梯度下降的代码
# напишем функцию градиентного спуска без использования библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = sum(x_us)
# сумма истинных ответов (выручка за весь период)
sy = sum(y_us)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x_us[i]*y_us[i]) for i in range(len(x_us))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x_us[i]**2) for i in range(len(x_us))]
sx_sq = sum(list_x_sq)
# количество значений
num = len(x_us)
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = [a,b]
errors.append(errors_sq_Kramer_method(ab,x_us,y_us))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
我们潜入马里亚纳海沟的最底部，在那里我们发现了所有相同的系数值 и ，这正是我们所期望的。
让我们再潜水一次，只是这一次，我们的深海飞行器将充满其他技术，即图书馆 NumPy的.
梯度下降代码 (NumPy)
# перед тем определить функцию для градиентного спуска с использованием библиотеки NumPy,
# напишем функцию определения суммы квадратов отклонений также с использованием NumPy
def error_square_numpy(ab,x_np,y_np):
y_pred = np.dot(x_np,ab)
error = y_pred - y_np
return sum((error)**2)
# напишем функцию градиентного спуска с использованием библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = float(sum(x_np[:,1]))
# сумма истинных ответов (выручка за весь период)
sy = float(sum(y_np))
# сумма произведения значений на истинные ответы
sxy = x_np*y_np
sxy = float(sum(sxy[:,1]))
# сумма квадратов значений
sx_sq = float(sum(x_np[:,1]**2))
# количество значений
num = float(x_np.shape[0])
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = np.array([[a],[b]])
errors.append(error_square_numpy(ab,x_np,y_np))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
系数值 и 不可改变的。
我们来看看梯度下降过程中误差如何变化，即偏差平方和如何随着每一步变化。
绘制偏差平方和的代码
print 'График№4 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()图 4“梯度下降期间偏差平方和”
在图表上我们看到，每一步误差都会减小，经过一定次数的迭代后，我们观察到一条几乎水平的线。
最后，我们来估计一下代码执行时间的差异：
确定梯度下降计算时间的代码
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска без использования библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска с использованием библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
也许我们做错了什么，但这又是一个简单的“自制”函数，不使用该库 NumPy的 优于使用库的函数的计算时间 NumPy的.
但我们并没有停滞不前，而是正在研究另一种令人兴奋的方法来求解简单的线性回归方程。 见面！
随机梯度下降
为了快速理解随机梯度下降的运行原理，最好确定它与普通梯度下降的区别。 我们，在梯度下降的情况下，在导数方程中 и 使用样本中可用的所有特征和真实答案的值的总和（即所有特征的总和 и ）。 在随机梯度下降中，我们不会使用样本中存在的所有值，而是伪随机地选择所谓的样本索引并使用其值。
例如，如果索引确定为数字 3（三），则我们取值 и ，然后我们将这些值代入导数方程并确定新的坐标。 然后，确定了坐标后，我们再次伪随机确定样本索引，将索引对应的值代入偏微分方程中，以新的方式确定坐标 и ETC。 直到收敛变成绿色。 乍一看，这似乎根本行不通，但确实行得通。 确实，值得注意的是，误差并不会随着每一步而减少，但肯定有一种趋势。
随机梯度下降与传统梯度下降相比有哪些优点？ 如果我们的样本量非常大并且以数万个值进行测量，那么处理随机的数千个值而不是整个样本要容易得多。 这就是随机梯度下降发挥作用的地方。 当然，就我们而言，我们不会注意到太大的差异。
让我们看一下代码。
随机梯度下降的代码
# определим функцию стох.град.шага
def stoch_grad_step_usual(vector_init, x_us, ind, y_us, l):
# выбираем значение икс, которое соответствует случайному значению параметра ind
# (см.ф-цию stoch_grad_descent_usual)
x = x_us[ind]
# рассчитывыаем значение y (выручку), которая соответствует выбранному значению x
y_pred = vector_init[0] + vector_init[1]*x_us[ind]
# вычисляем ошибку расчетной выручки относительно представленной в выборке
error = y_pred - y_us[ind]
# определяем первую координату градиента ab
grad_a = error
# определяем вторую координату ab
grad_b = x_us[ind]*error
# вычисляем новый вектор коэффициентов
vector_new = [vector_init[0]-l*grad_a, vector_init[1]-l*grad_b]
return vector_new
# определим функцию стох.град.спуска
def stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800):
# для самого начала работы функции зададим начальные значения коэффициентов
vector_init = [float(random.uniform(-0.5, 0.5)), float(random.uniform(-0.5, 0.5))]
errors = []
# запустим цикл спуска
# цикл расчитан на определенное количество шагов (steps)
for i in range(steps):
ind = random.choice(range(len(x_us)))
new_vector = stoch_grad_step_usual(vector_init, x_us, ind, y_us, l)
vector_init = new_vector
errors.append(errors_sq_Kramer_method(vector_init,x_us,y_us))
return (vector_init),(errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
我们仔细观察这些系数，不禁问自己“这怎么可能？” 我们得到了其他系数值 и 。 也许随机梯度下降已经为方程找到了更优化的参数？ 很不幸的是，不行。 只需查看偏差平方和即可发现，使用新的系数值，误差会更大。 我们并不急于绝望。 让我们构建一个误差变化图。
绘制随机梯度下降中偏差平方和的代码
print 'График №5 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()图 5“随机梯度下降过程中的偏差平方和”
看看日程安排，一切都已就位，现在我们将解决所有问题。
所以发生了什么事？ 发生了以下情况。 当我们随机选择一个月时，我们的算法会针对所选月份寻求减少计算收入的误差。 然后我们选择另一个月份并重复计算，但我们减少了第二个选定月份的误差。 现在请记住，前两个月明显偏离简单线性回归方程的直线。 这意味着，当选择这两个月中的任何一个月时，通过减少其中每个月的误差，我们的算法会严重增加整个样本的误差。 那么该怎么办？ 答案很简单：你需要减少下降步长。 毕竟，通过减少下降步长，误差也会停止上下“跳跃”。 或者更确切地说，“跳跃”错误不会停止，但不会那么快发生:)让我们检查一下。
以较小增量运行 SGD 的代码
# запустим функцию, уменьшив шаг в 100 раз и увеличив количество шагов соответсвующе
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print 'График №6 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
图 6“随机梯度下降期间的偏差平方和（80 万步）”
系数有所改善，但仍不理想。 假设，这可以通过这种方式纠正。 例如，我们在最后 1000 次迭代中选择误差最小的系数值。 确实，为此我们还必须写下系数本身的值。 我们不会这样做，而是要注意日程安排。 看起来很平滑，误差似乎均匀减少。 其实这不是真的。 让我们看一下前 1000 次迭代，并将它们与最后一次进行比较。
SGD 图表代码（前 1000 步）
print 'График №7 "Сумма квадратов отклонений по-шагово. Первые 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])),
list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
print 'График №7 "Сумма квадратов отклонений по-шагово. Последние 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])),
list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()图 7“偏差平方和 SGD（前 1000 步）”
图 8“偏差平方和 SGD（最后 1000 步）”
在下降的一开始，我们观察到误差相当均匀且急剧下降。 在最后一次迭代中，我们看到误差在 1,475 的值附近徘徊，在某些时刻甚至等于这个最佳值，但随后它仍然上升......我再说一遍，你可以写下 的值系数 и ，然后选择误差最小的那些。 然而，我们遇到了一个更严重的问题：我们必须执行 80 万步（参见代码）才能获得接近最优的值。 而这已经与相对于梯度下降用随机梯度下降节省计算时间的想法相矛盾了。 哪些方面可以纠正和改进？ 不难注意到，在第一次迭代中，我们自信地向下走，因此，我们应该在第一次迭代中留下一个大的步骤，并在前进时减少步骤。 我们不会在本文中这样做 - 它已经太长了。 有兴趣的可以自己想想怎么做，并不难：）
现在让我们使用该库执行随机梯度下降 NumPy的 （我们不要被我们之前确定的石头绊倒）
随机梯度下降代码 (NumPy)
# для начала напишем функцию градиентного шага
def stoch_grad_step_numpy(vector_init, X, ind, y, l):
x = X[ind]
y_pred = np.dot(x,vector_init)
err = y_pred - y[ind]
grad_a = err
grad_b = x[1]*err
return vector_init - l*np.array([grad_a, grad_b])
# определим функцию стохастического градиентного спуска
def stoch_grad_descent_numpy(X, y, l=0.1, steps = 800):
vector_init = np.array([[np.random.randint(X.shape[0])], [np.random.randint(X.shape[0])]])
errors = []
for i in range(steps):
ind = np.random.randint(X.shape[0])
new_vector = stoch_grad_step_numpy(vector_init, X, ind, y, l)
vector_init = new_vector
errors.append(error_square_numpy(vector_init,X,y))
return (vector_init), (errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print
结果与不使用下降时的值几乎相同 NumPy的。 然而，这是合乎逻辑的。
让我们看看随机梯度下降花了多长时间。
确定SGD计算时间的代码（80万步）
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска без использования библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска с использованием библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
越深入森林，云层就越暗：“自写”公式再次显示了最好的结果。 所有这些都表明必须有更微妙的方式来使用该库 NumPy的，这确实加快了计算操作。 在本文中我们不会了解它们。 闲暇时会有一些事情可以思考:)
我们总结一下
在进行总结之前，我想回答一个很可能是我们亲爱的读者提出的问题。 事实上，如果我们手中有如此强大而简单的设备，为什么我们需要在山上走上走下（大部分是下山）才能找到珍贵的低地？分析解决方案的形式，它立即将我们传送到正确的地方？
这个问题的答案就在表面上。 现在我们看一个非常简单的例子，其中真正的答案是 取决于一个标志 。 这种情况在生活中并不常见，所以让我们想象一下我们有 2 个、30 个、50 个或更多的标志。 让我们为每个属性添加数千甚至数万个值。 在这种情况下，解析解可能无法经受住考验而失败。 反过来，梯度下降及其变化将缓慢但肯定地使我们更接近目标——函数的最小值。 不要担心速度 - 我们可能会寻找允许我们设置和调节步长（即速度）的方法。
现在是实际的简要总结。
首先，我希望本文中提供的材料能够帮助刚入门的“数据科学家”了解如何求解简单（而不仅仅是）线性回归方程。
其次，我们研究了求解方程的几种方法。 现在，我们可以根据情况选择最适合解决问题的方案。
第三，我们看到了附加设置的威力，即梯度下降步长。 这个参数不能忽略。 如上所述，为了减少计算成本，应该在下降过程中改变步长。
第四，在我们的例子中，“自制”函数显示了最佳的计算时间结果。 这可能是由于没有最专业地使用图书馆的功能 NumPy的。 但无论如何，以下结论不言自明。 一方面，有时值得质疑既定的观点，另一方面，并不总是值得让一切变得复杂——相反，有时解决问题的更简单的方法更有效。 由于我们的目标是分析解决简单线性回归方程的三种方法，因此使用“自编写”函数对我们来说就足够了。
文学（或类似的东西）
1. 线性回归
2.最小二乘法
3. 衍生品
4。 梯度
5.梯度下降
6.NumPy库
来源： habr.com