🥇Решете просто уравнение на линейна регресия

Статията обсъжда няколко начина за определяне на математическото уравнение на проста (двойка) регресионна линия.

Всички методи за решаване на разглежданото тук уравнение се основават на метода на най-малките квадрати. Ние обозначаваме методите, както следва:

Аналитично решение
градиентно спускане
Стохастичен градиентен спускане

За всеки от начините за решаване на уравнението на права линия статията представя различни функции, които са разделени главно на тези, които са написани без използване на библиотеката numpy и тези, които се използват за изчисления numpy. Смята се, че умелото използване numpy ще намали цената на изчисленията.

Целият код в тази статия е написан на python 2.7 с Джупиър Бележник. Изходният код и примерният файл с данни са достъпни на Github

Статията е по-фокусирана както върху начинаещите, така и върху тези, които вече постепенно са започнали да овладяват изучаването на много обширен раздел в изкуствения интелект - машинното обучение.

Нека използваме един много прост пример, за да илюстрираме материала.

Примерни условия

Имаме пет ценности, които характеризират пристрастяването Y от X (Маса 1):

Таблица № 1 "Условия на примера"

Ще приемем, че стойностите е месецът от годината и - печалба този месец. С други думи, приходите зависят от месеца в годината и - единственият знак, от който зависят приходите.

Примерът е така, както по отношение на условната зависимост на приходите от месеца на годината, така и по отношение на броя на стойностите - има много малко от тях. Подобно опростяване обаче ще позволи, както се казва на пръстите, да обясни, не винаги с лекота, материала, усвоен от начинаещи. Освен това простотата на числата ще позволи на тези, които желаят да решат примера на „хартия“ без значителни разходи за труд.

Да предположим, че зависимостта, дадена в примера, може да бъде апроксимирана доста добре чрез математическото уравнение на линията на проста (двойка) регресия на формата:

където е месецът, през който са получени приходите, - приходи, съответстващи на месеца, и са регресионните коефициенти на оценената линия.

Имайте предвид, че коефициентът често се нарича наклон или градиент на прогнозната линия; е сумата, с която когато се промени .

Очевидно нашата задача в примера е да изберем такива коефициенти в уравнението и , при което отклоненията на прогнозните ни стойности на приходите по месеци от верните отговори, т.е. стойностите, представени в извадката, ще бъдат минимални.

Метод на най-малките квадрати

Според метода на най-малките квадрати отклонението трябва да се изчисли чрез повдигане на квадрат. Такава техника позволява да се избегне взаимното изплащане на отклонения, ако те имат противоположни знаци. Например, ако в един случай, отклонението е +5 (плюс пет), а в другата -5 (минус пет), тогава сумата на отклоненията ще се анулира взаимно и ще бъде 0 (нула). Възможно е отклонението да не се квадратизира, а да се използва свойството модул и тогава всички отклонения да са положителни и да се натрупват. Няма да се спираме подробно на тази точка, а просто посочваме, че за удобство на изчисленията е обичайно отклонението да се квадратира.

Ето как изглежда формулата, с помощта на която ще определим най-малката сума на квадратите на отклонения (грешки):

където е функция на приближаване на верните отговори (т.е. приходите, изчислени от нас),

са верните отговори (приходите, предоставени в извадката),

е индексът на извадката (номерът на месеца, в който е определено отклонението)

Нека диференцираме функцията, дефинираме частичните диференциални уравнения и сме готови да преминем към аналитично решение. Но първо, нека направим кратко отклонение относно това какво е диференциране и да си припомним геометричното значение на производната.

Диференциация

Диференцирането е операция за намиране на производната на функция.

За какво е производната? Производната на функция характеризира скоростта на изменение на функцията и ни показва нейната посока. Ако производната в дадена точка е положителна, тогава функцията е нарастваща; в противен случай функцията е намаляваща. И колкото по-голяма е стойността на производната по модул, толкова по-висока е скоростта на промяна на стойностите на функцията, както и толкова по-стръмен е наклонът на графиката на функцията.

Например при условията на декартова координатна система стойността на производната в точка M(0,0) е равна на +25 означава, че в дадена точка, когато стойността е изместена вдясно с условна единица, стойност се увеличава с 25 условни единици. На графиката изглежда като доста стръмен ъгъл на покачване на стойностите от дадена точка.

Друг пример. Стойността на производната е -0,1 означава, че при превключване за една условна единица, стойността намалява само с 0,1 условна единица. В същото време на графиката на функцията можем да наблюдаваме едва забележим наклон надолу. Правейки аналогия с планина, все едно много бавно се спускаме по лек склон от планина, за разлика от предишния пример, където трябваше да изкачваме много стръмни върхове :)

Така след обособяване на функцията по коефициенти и , дефинираме уравненията на частни производни от 1-ви ред. След като дефинираме уравненията, ще получим система от две уравнения, решавайки които можем да изберем такива стойности на коефициентите и , при което стойностите на съответните производни в дадени точки се променят с много, много малка сума, а в случай на аналитично решение не се променят изобщо. С други думи, функцията на грешката при намерените коефициенти ще достигне минимум, тъй като стойностите на частичните производни в тези точки ще бъдат равни на нула.

И така, съгласно правилата за диференциране, уравнението на частната производна от 1-ви ред по отношение на коефициента ще приеме формата:

Уравнение с частна производна от 1-ви ред по отношение на ще приеме формата:

В резултат на това получихме система от уравнения, която има доста просто аналитично решение:

начало{уравнение*}
начало{случаи}
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
край {случаи}
край{уравнение*}

Преди да решим уравнението, нека предварително заредим, проверим коректността на зареждането и форматираме данните.

Зареждане и форматиране на данни

Трябва да се отбележи, че поради факта, че за аналитичното решение и в бъдеще за градиентно и стохастично градиентно спускане, ще използваме кода в два варианта: използване на библиотеката 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 '********************************************'

Визуализация

Сега, след като, първо, заредихме данните, второ, проверихме коректността на зареждането и накрая форматирахме данните, ще извършим първата визуализация. Често се използва този метод парцел библиотека Морскородено. В нашия пример, поради ограничения брой, няма смисъл да се използва библиотеката Морскородено. Ще използваме обичайната библиотека Матплотлиб и гледайте само диаграмата на разсейване.

Код на точечната диаграма

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 "Зависимост на приходите от месеца на годината"

Аналитично решение

Използваме най-често срещаните инструменти в питон и реши системата от уравнения:

начало{уравнение*}
начало{случаи}
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
край {случаи}
край{уравнение*}

Според правилото на Креймър намерете обща детерминанта, както и детерминанти по и , след което, разделяйки детерминантата на към обща детерминанта - намерете коефициента , по подобен начин намираме коефициента .

Код на аналитично решение

# определим функцию для расчета коэффициентов 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 (минус едно).

Нека обърнем внимание на факта, че една функция може да има няколко минимума и след като се спуснем в един от тях според предложения по-долу алгоритъм, няма да можем да намерим друг минимум, който може да е по-нисък от намерения. Спокойно, не сме в опасност! В нашия случай имаме работа с един минимум, тъй като нашата функция на графиката е обикновена парабола. И както всички трябва да знаем много добре от училищния курс по математика, параболата има само един минимум.

След като разбрахме защо имаме нужда от градиент и също, че градиентът е сегмент, тоест вектор с дадени координати, които са същите коефициенти и можем да приложим градиентно спускане.

Преди да започнете, предлагам да прочетете само няколко изречения за алгоритъма за спускане:

Определяме псевдослучайно координатите на коефициентите и . В нашия пример ще дефинираме коефициенти близо до нула. Това е обичайна практика, но всеки случай може да има своя собствена практика.
От координата извадете стойността на частната производна от 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 (три), тогава ние вземаме стойностите и , след това заместваме стойностите в уравненията на производните и определяме нови координати. След това, след като определихме координатите, ние отново определяме примерния индекс по псевдослучаен начин, заместваме стойностите, съответстващи на индекса, в частичните диференциални уравнения и определяме координатите по нов начин и и т.н. до зелена конвергенция. На пръв поглед може да изглежда, че изобщо може да работи, но работи. Вярно е, че трябва да се отбележи, че грешката не намалява с всяка стъпка, но със сигурност има тенденция.

Какви са предимствата на стохастичното градиентно спускане пред конвенционалното градиентно спускане? Ако размерът на нашата извадка е много голям и се измерва в десетки хиляди стойности, тогава е много по-лесно да обработим, да речем, произволни хиляда от тях, отколкото цялата извадка. Това е мястото, където стохастичното градиентно спускане започва. В нашия случай, разбира се, няма да забележим голяма разлика.

Нека да разгледаме кода.

Код за стохастичен градиентен спускане

# определим функцию стох.град.шага
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 „Сума от квадратни отклонения за стохастичен градиент надолу (80k стъпки)“

Коефициентите се подобриха, но все още не са идеални. Хипотетично това може да се коригира по този начин. Избираме, например, на последните 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. Но както и да е, заключението се налага следното. От една страна, понякога си струва да поставите под въпрос установените мнения, а от друга страна, не винаги си струва да усложнявате всичко - напротив, понякога по-простият начин за решаване на проблем е по-ефективен. И тъй като целта ни беше да анализираме три подхода за решаване на просто уравнение на линейна регресия, използването на „самоизписани“ функции беше достатъчно за нас.