🥇Решавање једначине једноставне линеарне регресије

У чланку се разматра неколико начина за одређивање математичке једначине једноставне (упарене) регресионе линије.

Све методе решавања једначине о којима се овде говори заснивају се на методи најмањих квадрата. Означимо методе на следећи начин:

Аналитичко решење
Градиент Десцент
Стохастички градијентни пад

За сваки метод решавања једначине праве, чланак даје различите функције, које се углавном деле на оне које су написане без коришћења библиотеке НумПи и оне које користе за прорачуне НумПи. Верује се да је вешта употреба НумПи смањиће трошкове рачунара.

Сав код дат у чланку је написан на језику питхон 2.7 Користећи Јупитер Нотебоок. Изворни код и датотека са узорцима података се објављују на Гитхуб

Чланак је више намењен и почетницима и онима који су већ постепено почели да савладавају проучавање веома широког дела вештачке интелигенције - машинског учења.

За илустрацију материјала користимо врло једноставан пример.

Пример услова

Имамо пет вредности које карактеришу зависност Y из X (Табела бр. 1):

Табела бр. 1 „Пример услова”

Претпоставићемо да су вредности је месец у години, и — приход овог месеца. Другим речима, приход зависи од месеца у години, и - једини знак од кога зависи приход.

Пример је тако-тако, како са становишта условне зависности прихода од месеца у години, тако и са становишта броја вредности - врло их је мало. Међутим, такво поједностављење ће омогућити, како кажу, да се објасни, не увек са лакоћом, материјал који почетници усвајају. А такође, једноставност бројева ће омогућити онима који желе да реше пример на папиру без значајних трошкова рада.

Претпоставимо да се зависност дата у примеру може прилично добро апроксимирати математичком једначином једноставне (упарене) регресионе линије облика:

где је месец у коме је приход примљен, — приход који одговара месецу, и су коефицијенти регресије процењене линије.

Имајте на уму да коефицијент често се назива нагиб или градијент процењене линије; представља износ којим се када се промени .

Очигледно, наш задатак у примеру је да изаберемо такве коефицијенте у једначини и , при чему су одступања наших израчунатих вредности прихода по месецима од тачних одговора, тј. вредности представљене у узорку ће бити минималне.

Метод најмањег квадрата

Према методи најмањих квадрата, одступање треба израчунати квадрирањем. Ова техника вам омогућава да избегнете међусобно отказивање одступања ако имају супротне знаке. На пример, ако је у једном случају одступање +5 (плус пет), а у другом -5 (минус пет), тада ће се збир одступања међусобно поништавати и износити 0 (нула). Могуће је не квадратирати одступање, већ користити својство модула и тада ће сва одступања бити позитивна и акумулирати. Нећемо се детаљно задржавати на овој тачки, већ једноставно назначити да је због погодности прорачуна уобичајено да се одступање квадратира.

Овако изгледа формула којом ћемо одредити најмањи збир квадрата одступања (грешке):

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

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

је индекс узорка (број месеца у коме је утврђено одступање)

Хајде да издиференцирамо функцију, дефинишемо парцијалне диференцијалне једначине и будимо спремни да пређемо на аналитичко решење. Али прво, хајде да направимо кратак излет о томе шта је диференцијација и присетимо се геометријског значења изведенице.

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

Диференцијација је операција проналажења извода функције.

За шта се користи дериват? Извод функције карактерише брзину промене функције и говори нам њен правац. Ако је извод у датој тачки позитиван, онда се функција повећава; у супротном, функција опада. А што је већа вредност апсолутног извода, то је већа брзина промене вредности функције, као и стрмији нагиб графика функције.

На пример, под условима Декартовог координатног система, вредност извода у тачки М(0,0) је једнака +25 значи да у датој тачки, када се вредност помери десно по конвенционалној јединици, вредност повећава за 25 конвенционалних јединица. На графикону то изгледа као прилично стрм пораст вредности из дате тачке.

Други пример. Вредност деривата је једнака -КСНУМКС значи да када се расељава по једној конвенционалној јединици, вредност смањује се за само 0,1 конвенционалну јединицу. Истовремено, на графику функције можемо уочити једва приметан нагиб надоле. Поводећи аналогију са планином, као да се веома полако спуштамо низ благу падину са планине, за разлику од претходног примера где смо морали да се пењемо на веома стрме врхове :)

Дакле, након диференцирања функције по изгледу и , дефинишемо парцијалне диференцијалне једначине 1. реда. Након одређивања једначина, добићемо систем од две једначине, решавањем којих ћемо моћи да изаберемо такве вредности коефицијената и , за које се вредности одговарајућих извода у датим тачкама мењају за веома, веома мали износ, а у случају аналитичког решења уопште се не мењају. Другим речима, функција грешке код пронађених коефицијената ће достићи минимум, пошто ће вредности парцијалних извода у овим тачкама бити једнаке нули.

Дакле, према правилима диференцијације, једначина парцијалног деривата 1. реда у односу на коефицијент ће имати облик:

Једначина парцијалног извода 1. реда у односу на ће имати облик:

Као резултат, добили смо систем једначина који има прилично једноставно аналитичко решење:

почети{једначина*}
почети{случајеви}
на + бсумлимитс_{и=1}^нк_и — сумлимитс_{и=1}^ни_и = 0

сумлимитс_{и=1}^нк_и(а +бсумлимитс_{и=1}^нк_и — сумлимитс_{и=1}^ни_и) = 0
крај{цасес}
крај{једначина*}

Пре него што решимо једначину, хајде да унапред учитамо, проверимо да ли је учитавање исправно и форматирамо податке.

Учитавање и форматирање података

Треба напоменути да ћемо због чињенице да ћемо за аналитичко решење, а потом и за градијент и стохастички градијентни спуст, користити код у две варијације: коришћење библиотеке НумПи и без употребе, биће нам потребно одговарајуће форматирање података (погледајте код).

Код за учитавање и обраду података

# импортируем все нужные нам библиотеки
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 „Зависност прихода од месеца у години“

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

Хајде да користимо најчешће алате у питон и реши систем једначина:

почети{једначина*}
почети{случајеви}
на + бсумлимитс_{и=1}^нк_и — сумлимитс_{и=1}^ни_и = 0

сумлимитс_{и=1}^нк_и(а +бсумлимитс_{и=1}^нк_и — сумлимитс_{и=1}^ни_и) = 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 „Одступања, %“

Није савршено, али смо завршили свој задатак.

Напишимо функцију која, да одредимо коефицијенте и користи библиотеку НумПи, тачније, написаћемо две функције: једну користећи псеудоинверзну матрицу (не препоручује се у пракси, пошто је процес рачунарски сложен и нестабилан), другу користећи матричну једначину.

Код аналитичког решења (НумПи)

# для начала добавим столбец с не изменяющимся значением в 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)

Са малом количином података, напред излази „самонаписана“ функција која проналази коефицијенте користећи Црамеров метод.

Сада можете прећи на друге начине проналажења коефицијената и .

Градиент Десцент

Прво, хајде да дефинишемо шта је градијент. Једноставно речено, градијент је сегмент који указује на правац максималног раста функције. По аналогији са пењањем на планину, где је успон тамо где је најстрмији успон на врх планине. Развијајући пример са планином, сетимо се да нам је у ствари потребан најстрмији спуст да бисмо што брже стигли до низије, односно минимума – места где се функција не повећава нити смањује. У овом тренутку извод ће бити једнак нули. Дакле, не треба нам градијент, већ антиградијент. Да бисте пронашли антиградијент, само треба да помножите градијент са -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
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)

Можда радимо нешто погрешно, али опет то је једноставна „кућно написана“ функција која не користи библиотеку НумПи надмашује време израчунавања функције која користи библиотеку НумПи.

Али ми не стојимо на месту, већ се крећемо ка проучавању још једног узбудљивог начина решавања једноставне једначине линеарне регресије. Сусрет!

Стохастички градијентни пад

Да би се брзо разумео принцип рада стохастичког градијентног спуштања, боље је утврдити његове разлике од обичног градијентног спуштања. Ми, у случају градијентног спуштања, у једначинама извода од и користио збир вредности свих карактеристика и тачних одговора доступних у узорку (односно збира свих и ). У стохастичком градијентном спуштању, нећемо користити све вредности присутне у узорку, већ уместо тога, псеудо-случајно бирамо такозвани индекс узорка и користимо његове вредности.

На пример, ако се одреди да је индекс број 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 „Збир одступања на квадрат током стохастичког градијента спуштања“

Гледајући распоред, све долази на своје место и сада ћемо све поправити.

И шта се десило? Десило се следеће. Када насумично одаберемо месец, онда наш алгоритам за изабрани месец настоји да смањи грешку у израчунавању прихода. Затим бирамо други месец и понављамо обрачун, али смањујемо грешку за други изабрани месец. Сада запамтите да прва два месеца значајно одступају од линије једноставне једначине линеарне регресије. То значи да када се изабере било који од ова два месеца, смањењем грешке сваког од њих, наш алгоритам озбиљно повећава грешку за цео узорак. Па ста да радим? Одговор је једноставан: морате смањити корак спуштања. На крају крајева, смањењем корака спуштања, грешка ће такође престати да "скаче" горе-доле. Тачније, грешка „скакање“ неће престати, али неће то учинити тако брзо :) Хајде да проверимо.

Код за покретање СГД-а са мањим корацима

# запустим функцию, уменьшив шаг в 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 итерација и упоредимо их са последњим.

Код за СГД графикон (првих 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 „Збир квадратних девијација СГД (првих 1000 корака)”

Графикон бр. 8 „Збир квадратних девијација СГД (последњих 1000 корака)”

На самом почетку спуштања примећујемо прилично равномерно и стрмо смањење грешке. У последњим итерацијама видимо да грешка иде около и око вредности од 1,475 и у неким тренуцима се чак изједначи са овом оптималном вредношћу, али онда и даље расте... Понављам, можете да запишете вредности коефицијенти и , а затим изаберите оне за које је грешка минимална. Међутим, имали смо озбиљнији проблем: морали смо да предузмемо 80 хиљада корака (погледајте код) да бисмо вредности приближили оптималним. А ово је већ у супротности са идејом да се уштеди време рачунања са стохастичким градијентом спуштања у односу на градијентни пад. Шта се може исправити и побољшати? Није тешко приметити да у првим итерацијама самопоуздано идемо наниже и, стога, треба оставити велики корак у првим итерацијама и смањивати корак како идемо напред. Нећемо то радити у овом чланку - већ је предугачак. Они који желе могу сами да смисле како то да ураде, није тешко :)

Хајде сада да извршимо стохастички градијентни спуст користећи библиотеку НумПи (и да се не спотичемо о камење које смо раније идентификовали)

Код за стохастички градијентни спуштање (НумПи)

# для начала напишем функцию градиентного шага
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

Испоставило се да су вредности скоро исте као при спуштању без употребе НумПи. Међутим, ово је логично.

Хајде да сазнамо колико нам је требало стохастички градијентни спуст.

Код за одређивање времена израчунавања СГД (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)

Што даље у шуму, то су облаци тамнији: опет, „самонаписана“ формула показује најбољи резултат. Све ово сугерише да морају постојати још суптилнији начини коришћења библиотеке НумПи, што заиста убрзава рачунске операције. У овом чланку нећемо научити о њима. Имаћете о чему да размишљате у слободно време :)

Резимирајте

Пре него што сумирам, желео бих да одговорим на питање које је највероватније настало од нашег драгог читаоца. Чему, у ствари, такво „мучење“ спуштањима, зашто морамо да шетамо горе-доле (углавном доле) да бисмо пронашли драгоцену низију, ако у рукама имамо тако моћну и једноставну справу, у облик аналитичког решења, које нас моментално телепортује на право место?

Одговор на ово питање лежи на површини. Сада смо погледали веома једноставан пример, у коме је прави одговор зависи од једног знака . Ово се не виђа често у животу, па замислимо да имамо 2, 30, 50 или више знакова. Додајмо овоме хиљаде, или чак десетине хиљада вредности за сваки атрибут. У овом случају, аналитичко решење можда неће издржати тест и неће успети. Заузврат, градијентни спуст и његове варијације ће нас полако али сигурно довести ближе циљу - минимуму функције. И не брините о брзини - вероватно ћемо погледати начине који ће нам омогућити да подесимо и регулишемо дужину корака (тј. брзину).

А сада прави кратак резиме.

Прво, надам се да ће материјал представљен у чланку помоћи почетницима „научника података“ у разумевању како да реше једноставне (и не само) једначине линеарне регресије.

Друго, погледали смо неколико начина да решимо једначину. Сада, у зависности од ситуације, можемо изабрати ону која је најпогоднија за решавање проблема.

Треће, видели смо моћ додатних подешавања, односно дужину корака спуштања у градијенту. Овај параметар се не може занемарити. Као што је горе наведено, да би се смањили трошкови прорачуна, дужину корака треба променити током спуштања.

Четврто, у нашем случају, „кућно писане“ функције су показале најбоље временске резултате за прорачуне. Ово је вероватно због не професионалне употребе библиотечких могућности НумПи. Али како год било, намеће се следећи закључак. С једне стране, понекад је вредно преиспитивати устаљена мишљења, а са друге стране, не вреди увек све компликовати – напротив, понекад је ефикаснији једноставнији начин решавања проблема. А пошто је наш циљ био да анализирамо три приступа решавању једноставне линеарне регресионе једначине, употреба „самонаписаних“ функција нам је била сасвим довољна.