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

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

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

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

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

Сите кодови дадени во статијата се напишани на јазикот питон 2.7 користење Upупитер тетратка. Изворниот код и датотеката со примерок на податоци се објавени на 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
крај{случај}
крај{равенка*}

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

Вчитување и форматирање податоци

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

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

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

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

Ајде да ги користиме најчестите алатки во python и реши го системот на равенки:

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

# для начала добавим столбец с не изменяющимся значением в 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, 
# напишем функцию определения суммы квадратов отклонений также с использованием 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 „Збир на квадратни отстапувања за време на стохастичко спуштање на градиент“

Гледајќи го распоредот, се си доаѓа на свое место и сега сè ќе поправиме.

Што се случи? Се случи следново. Кога по случаен избор избираме месец, тогаш нашиот алгоритам се обидува да ја намали грешката во пресметувањето на приходот за избраниот месец. Потоа избираме уште еден месец и ја повторуваме пресметката, но ја намалуваме грешката за вториот избран месец. Сега запомнете дека првите два месеци значително отстапуваат од линијата на едноставната линеарна регресивна равенка. Ова значи дека кога ќе се избере некој од овие два месеци, со намалување на грешката на секој од нив, нашиот алгоритам сериозно ја зголемува грешката за целиот примерок. Па што да се прави? Одговорот е едноставен: треба да го намалите чекорот на спуштање. На крајот на краиштата, со намалување на чекорот за спуштање, грешката исто така ќе престане да „скока“ нагоре и надолу. Или подобро кажано, грешката „скокање“ нема да престане, но нема да го направи тоа толку брзо :) Ајде да провериме.

Код за извршување на 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)

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

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

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

Код за одредување на времето за пресметување на 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)

Колку подалеку во шумата, толку се потемни облаците: повторно, формулата „само напишана“ го покажува најдобриот резултат. Сето ова сугерира дека мора да има уште посуптилни начини за користење на библиотеката Нуп, кои навистина ги забрзуваат пресметковните операции. Во оваа статија нема да научиме за нив. Ќе има за што да размислувате во слободното време :)

Да резимираме

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

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

И сега вистинското кратко резиме.

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

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

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

Четврто, во нашиот случај, функциите „домашно напишани“ покажаа најдобри временски резултати за пресметките. Ова веројатно се должи на ненајпрофесионалното користење на можностите на библиотеката Нуп. Но, како и да е, следниов заклучок се сугерира сам по себе. Од една страна, понекогаш вреди да се преиспитуваат воспоставените мислења, а од друга страна, не вреди секогаш да се комплицира сè - напротив, понекогаш е поефективен поедноставниот начин на решавање на проблемот. И бидејќи нашата цел беше да анализираме три пристапи за решавање на едноставна линеарна регресивна равенка, употребата на „само напишани“ функции ни беше сосема доволна.

Литература (или нешто слично)

1. Линеарна регресија

http://statistica.ru/theory/osnovy-lineynoy-regressii/

2. Метод на најмали квадрати

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Дериват

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Градиент

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Спуштање на градиент

habr.com/ru/post/471458

habr.com/ru/post/307312

artemarakcheev.com//2017-12-31/linear_regression

6. NumPy библиотека

docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.linalg.solve.html

docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.pinv.html

pythonworld.ru/numpy/2.html

Извор: www.habr.com