🥇Lösa ekvationen för enkel linjär regression

Artikeln diskuterar flera sätt att bestämma den matematiska ekvationen för en enkel (parad) regressionslinje.

Alla metoder för att lösa ekvationen som diskuteras här är baserade på minsta kvadratmetoden. Låt oss beteckna metoderna enligt följande:

Analytisk lösning
Gradient Descent
Stokastisk gradientnedstigning

För varje metod för att lösa ekvationen för en rät linje ger artikeln olika funktioner, som huvudsakligen är uppdelade i de som är skrivna utan att använda biblioteket numpy och de som används för beräkningar numpy. Man tror att skicklig användning numpy kommer att minska datorkostnaderna.

All kod som ges i artikeln är skriven på språket python 2.7 med Jupyter Notebook. Källkoden och filen med exempeldata läggs ut på Github

Artikeln riktar sig mer till både nybörjare och de som redan så smått börjat bemästra studiet av ett mycket brett avsnitt inom artificiell intelligens – maskininlärning.

För att illustrera materialet använder vi ett mycket enkelt exempel.

Exempel på villkor

Vi har fem värden som kännetecknar beroende Y från X (Tabell nr 1):

Tabell nr 1 "Exempel på villkor"

Vi kommer att anta att värdena är månaden på året, och — intäkter denna månad. Med andra ord beror intäkterna på månaden på året, och - det enda tecken som intäkterna beror på.

Exemplet är så som så, både ur synvinkeln av inkomstens villkorliga beroende av månaden på året och ur synvinkeln av antalet värden - det finns väldigt få av dem. Men en sådan förenkling kommer att göra det möjligt, som de säger, att förklara, inte alltid med lätthet, materialet som förvärvats av nybörjare. Och även siffrornas enkelhet kommer att tillåta dem som vill lösa exemplet på papper utan betydande arbetskostnader.

Låt oss anta att beroendet som ges i exemplet kan approximeras ganska väl genom den matematiska ekvationen av en enkel (parad) regressionslinje av formen:

där är den månad då intäkterna erhölls, — inkomst motsvarande månaden, и är regressionskoefficienterna för den uppskattade linjen.

Observera att koefficienten kallas ofta lutningen eller gradienten för den uppskattade linjen; representerar det belopp med vilket när det ändras .

Uppenbarligen är vår uppgift i exemplet att välja sådana koefficienter i ekvationen и , där avvikelserna för våra beräknade intäktsvärden per månad från de sanna svaren, dvs. värden som presenteras i provet kommer att vara minimala.

Minsta kvadratiska metod

Enligt minsta kvadratmetoden ska avvikelsen beräknas genom att kvadrera den. Denna teknik låter dig undvika ömsesidig upphävande av avvikelser om de har motsatta tecken. Till exempel, om i ett fall avvikelsen är +5 (plus fem), och i den andra -5 (minus fem), då kommer summan av avvikelserna att ta bort varandra och uppgå till 0 (noll). Du behöver inte kvadrera avvikelsen utan använd modulegenskapen och då blir alla avvikelser positiva och ackumuleras. Vi kommer inte att uppehålla oss i denna punkt i detalj, utan bara indikera att det för beräkningarnas bekvämlighet är vanligt att kvadrera avvikelsen.

Så här ser formeln ut med vilken vi kommer att bestämma den minsta summan av kvadrerade avvikelser (fel):

där är en funktion av approximation av sanna svar (det vill säga intäkterna vi beräknade),

är de sanna svaren (intäkter som anges i urvalet),

är provindexet (nummer för den månad då avvikelsen bestäms)

Låt oss differentiera funktionen, definiera de partiella differentialekvationerna och vara redo att gå vidare till den analytiska lösningen. Men först, låt oss ta en kort utflykt om vad differentiering är och komma ihåg den geometriska betydelsen av derivatan.

Differentiering

Differentiering är operationen att hitta derivatan av en funktion.

Vad används derivatan till? En funktions derivata kännetecknar funktionens förändringshastighet och berättar för oss dess riktning. Om derivatan vid en given punkt är positiv ökar funktionen, annars minskar funktionen. Och ju större den absoluta derivatans värde är, desto högre är förändringshastigheten för funktionsvärdena, liksom desto brantare lutning på funktionsgrafen.

Till exempel, under villkoren för ett kartesiskt koordinatsystem, är värdet på derivatan vid punkten M(0,0) lika med +25 betyder att vid en given punkt, när värdet ändras till höger av en konventionell enhet, värde ökar med 25 konventionella enheter. På grafen ser det ut som en ganska brant värdestegring från en given punkt.

Ett annat exempel. Det derivata värdet är lika -0,1 betyder att när de förflyttas per en konventionell enhet, värde minskar med endast 0,1 konventionell enhet. Samtidigt kan vi på funktionens graf observera en knappt märkbar nedåtlutning. Om vi drar en analogi med ett berg, är det som om vi mycket långsamt går ner för en svag sluttning från ett berg, till skillnad från föregående exempel, där vi var tvungna att bestiga väldigt branta toppar :)

Alltså efter differentiering av funktionen efter odds и , definierar vi 1:a ordningens partiella differentialekvationer. Efter att ha bestämt ekvationerna kommer vi att få ett system med två ekvationer, genom att lösa vilka vi kommer att kunna välja sådana värden av koefficienterna и , för vilka värdena för motsvarande derivator vid givna punkter ändras med en mycket, mycket liten mängd, och i fallet med en analytisk lösning inte förändras alls. Med andra ord kommer felfunktionen vid de hittade koefficienterna att nå ett minimum, eftersom värdena för partiella derivator vid dessa punkter kommer att vara lika med noll.

Så, enligt reglerna för differentiering, den partiella derivatekvationen av första ordningen med avseende på koefficienten kommer att ha formen:

1:a ordningens partiella derivatekvation med avseende på kommer att ha formen:

Som ett resultat fick vi ett ekvationssystem som har en ganska enkel analytisk lösning:

börja{ekvation*}
börja{cases}
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
slut{cases}
slut{ekvation*}

Innan vi löser ekvationen, låt oss ladda i förväg, kontrollera att laddningen är korrekt och formatera data.

Laddar och formaterar data

Det bör noteras att på grund av det faktum att för den analytiska lösningen, och därefter för gradient och stokastisk gradientnedstigning, kommer vi att använda koden i två varianter: med hjälp av biblioteket numpy och utan att använda det kommer vi att behöva lämplig dataformatering (se kod).

Dataladdnings- och bearbetningskod

# импортируем все нужные нам библиотеки
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 '********************************************'

Visualisering

Nu, efter att vi för det första har laddat data, för det andra kontrollerat korrektheten av inläsningen och slutligen formaterat data, kommer vi att utföra den första visualiseringen. Metoden som ofta används för detta är parplot bibliotek sjöfödd. I vårt exempel, på grund av det begränsade antalet, är det ingen idé att använda biblioteket sjöfödd. Vi kommer att använda det vanliga biblioteket matplotlib och titta bara på scatterplotten.

Scatterplot-kod

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()

Diagram nr 1 "Beroende av intäkter på månaden på året"

Analytisk lösning

Låt oss använda de vanligaste verktygen i pytonorm och lösa ekvationssystemet:

börja{ekvation*}
börja{cases}
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
slut{cases}
slut{ekvation*}

Enligt Cramers regel vi kommer att hitta den allmänna determinanten, såväl som determinanter av och , varefter determinanten divideras med till den allmänna determinanten - hitta koefficienten , på samma sätt finner vi koefficienten .

Kod för analytisk lösning

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

Här är vad vi fick:

Så, värdena på koefficienterna har hittats, summan av de kvadratiska avvikelserna har fastställts. Låt oss rita en rak linje på spridningshistogrammet i enlighet med de hittade koefficienterna.

Regressionslinjekod

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

Diagram nr 2 "Rätta och beräknade svar"

Du kan titta på avvikelsegrafen för varje månad. I vårt fall kommer vi inte att härleda något väsentligt praktiskt värde av det, men vi kommer att tillfredsställa vår nyfikenhet om hur väl den enkla linjära regressionsekvationen karakteriserar inkomstens beroende av månaden på året.

Avvikelsediagramkod

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

Diagram nr 3 "Avvikelser, %"

Inte perfekt, men vi klarade vår uppgift.

Låt oss skriva en funktion som för att bestämma koefficienterna и använder biblioteket numpy, mer exakt kommer vi att skriva två funktioner: den ena använder en pseudoinvers matris (rekommenderas inte i praktiken, eftersom processen är beräkningsmässigt komplex och instabil), den andra använder en matrisekvation.

Analytisk lösningskod (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

Låt oss jämföra tiden som ägnas åt att bestämma koefficienterna и , i enlighet med de 3 presenterade metoderna.

Kod för beräkning av beräkningstid

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)

Med en liten mängd data kommer en "självskriven" funktion fram, som hittar koefficienterna med Cramers metod.

Nu kan du gå vidare till andra sätt att hitta koefficienter и .

Gradient Descent

Låt oss först definiera vad en gradient är. Enkelt uttryckt är gradienten ett segment som anger riktningen för maximal tillväxt för en funktion. I analogi med att bestiga ett berg, där gradienten är där den brantaste stigningen till toppen av berget är. När vi utvecklar exemplet med berget kommer vi ihåg att vi faktiskt behöver den brantaste nedstigningen för att nå låglandet så snabbt som möjligt, det vill säga minimum - platsen där funktionen inte ökar eller minskar. Vid denna tidpunkt kommer derivatan att vara lika med noll. Därför behöver vi inte en gradient, utan en antigradient. För att hitta antigradienten behöver du bara multiplicera gradienten med -1 (minus ett).

Låt oss vara uppmärksamma på det faktum att en funktion kan ha flera minima, och efter att ha gått ner i en av dem med den algoritm som föreslås nedan, kommer vi inte att kunna hitta ett annat minimum, som kan vara lägre än det som hittades. Låt oss slappna av, detta är inget hot mot oss! I vårt fall har vi att göra med ett enda minimum, eftersom vår funktion på grafen finns en vanlig parabel. Och som vi alla borde veta mycket väl från vår skolmatematikkurs, har en parabel bara ett minimum.

Efter att vi fick reda på varför vi behövde en gradient, och även att gradienten är ett segment, det vill säga en vektor med givna koordinater, som är exakt samma koefficienter и vi kan implementera gradientnedstigning.

Innan du börjar föreslår jag att du bara läser några meningar om nedstigningsalgoritmen:

Vi bestämmer på ett pseudoslumpmässigt sätt koefficienternas koordinater и . I vårt exempel kommer vi att definiera koefficienter nära noll. Detta är en vanlig praxis, men varje fall kan ha sin egen praxis.
Från koordinat subtrahera värdet av 1:a ordningens partiella derivata vid punkten . Så om derivatan är positiv ökar funktionen. Därför, genom att subtrahera värdet på derivatan, kommer vi att röra oss i motsatt riktning av tillväxten, det vill säga i nedstigningsriktningen. Om derivatan är negativ, så minskar funktionen vid denna punkt och genom att subtrahera värdet på derivatan rör vi oss i nedstigningsriktningen.
Vi genomför en liknande operation med koordinaten : subtrahera värdet av den partiella derivatan vid punkten .
För att inte hoppa över miniminivån och flyga ut i rymden är det nödvändigt att ställa in stegstorleken i nedstigningsriktningen. I allmänhet kan du skriva en hel artikel om hur du ställer in steget korrekt och hur du ändrar det under nedstigningsprocessen för att minska beräkningskostnaderna. Men nu har vi en lite annan uppgift framför oss, och vi kommer att fastställa stegstorleken med hjälp av den vetenskapliga metoden "poke" eller, som man säger i vanligt språkbruk, empiriskt.
När vi är från de givna koordinaterna и subtrahera värdena på derivatorna får vi nya koordinater и . Vi tar nästa steg (subtraktion), redan från de beräknade koordinaterna. Och så startar cykeln om och om igen, tills den nödvändiga konvergensen uppnås.

Allt! Nu är vi redo att ge oss ut på jakt efter Mariangravens djupaste ravin. Låt oss börja.

Kod för gradientnedstigning

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

Vi dök till botten av Mariangraven och där hittade vi alla samma koefficientvärden и , vilket är precis vad som var att vänta.

Låt oss ta ett dyk till, men den här gången kommer vårt djuphavsfordon att fyllas med andra teknologier, nämligen ett bibliotek numpy.

Kod för gradientnedstigning (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

Koefficientvärden и oföränderlig.

Låt oss titta på hur felet förändrades under gradientnedstigning, det vill säga hur summan av kvadrerade avvikelser förändrades med varje steg.

Kod för att plotta summor av kvadrerade avvikelser

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()

Diagram nr 4 "Summa av kvadrerade avvikelser under gradientnedstigning"

På grafen ser vi att för varje steg minskar felet, och efter ett visst antal iterationer observerar vi en nästan horisontell linje.

Låt oss slutligen uppskatta skillnaden i kodexekveringstiden:

Kod för att bestämma tid för beräkning av gradientnedstigning

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)

Vi kanske gör något fel, men återigen är det en enkel "hemskriven" funktion som inte använder biblioteket numpy överträffar beräkningstiden för en funktion som använder biblioteket numpy.

Men vi står inte stilla, utan går mot att studera ett annat spännande sätt att lösa den enkla linjära regressionsekvationen. Möt oss!

Stokastisk gradientnedstigning

För att snabbt förstå principen om drift av stokastisk gradientnedstigning är det bättre att bestämma dess skillnader från vanlig gradientnedstigning. Vi, i fallet med gradientnedstigning, i ekvationerna av derivator av и använde summan av värdena för alla funktioner och sanna svar som finns tillgängliga i urvalet (det vill säga summan av alla и ). I stokastisk gradientnedstigning kommer vi inte att använda alla värden som finns i provet, utan istället pseudo-slumpmässigt välja det så kallade provindexet och använda dess värden.

Till exempel, om indexet bestäms vara nummer 3 (tre), så tar vi värdena и , sedan ersätter vi värdena i derivatekvationerna och bestämmer nya koordinater. Sedan, efter att ha bestämt koordinaterna, bestämmer vi återigen pseudo-slumpmässigt provindexet, ersätter värdena som motsvarar indexet i de partiella differentialekvationerna och bestämmer koordinaterna på ett nytt sätt и etc. tills konvergensen blir grön. Vid första anblicken kanske det inte verkar som att det här skulle kunna fungera alls, men det gör det. Det är sant att det är värt att notera att felet inte minskar för varje steg, men det finns säkert en tendens.

Vilka är fördelarna med stokastisk gradientnedstigning jämfört med konventionell? Om vår urvalsstorlek är mycket stor och mäts i tiotusentals värden, är det mycket lättare att bearbeta, säg, ett slumpmässigt tusental av dem, snarare än hela urvalet. Det är här stokastisk gradientnedstigning kommer in i bilden. I vårt fall kommer vi naturligtvis inte att märka någon större skillnad.

Låt oss titta på koden.

Kod för stokastisk gradientnedstigning

# определим функцию стох.град.шага
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])

Vi tittar noga på koefficienterna och får oss själva att ställa frågan "Hur kan detta vara?" Vi har andra koefficientvärden и . Kanske har stokastisk gradientnedstigning hittat mer optimala parametrar för ekvationen? Tyvärr inte. Det räcker med att titta på summan av kvadrerade avvikelser och se att med nya värden på koefficienterna är felet större. Vi har ingen brådska att misströsta. Låt oss bygga en graf över felförändringen.

Kod för att plotta summan av kvadratiska avvikelser i stokastisk gradientnedstigning

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()

Diagram nr 5 "Summa av kvadrerade avvikelser under stokastisk gradientnedstigning"

Tittar man på schemat så faller allt på plats och nu ska vi fixa allt.

Så vad hände? Följande hände. När vi slumpmässigt väljer en månad, är det för den valda månaden som vår algoritm försöker minska felet vid beräkning av intäkter. Sedan väljer vi ytterligare en månad och upprepar beräkningen, men vi minskar felet för den andra valda månaden. Kom nu ihåg att de första två månaderna avviker signifikant från linjen i den enkla linjära regressionsekvationen. Detta innebär att när någon av dessa två månader väljs, genom att minska felet för var och en av dem, ökar vår algoritm allvarligt felet för hela urvalet. Så vad ska man göra? Svaret är enkelt: du måste minska nedstigningssteget. När allt kommer omkring, genom att minska nedstigningssteget kommer felet också att sluta "hoppa" upp och ner. Eller snarare, "hoppnings"-felet kommer inte att sluta, men det kommer inte att göra det så snabbt :) Låt oss kolla.

Kod för att köra SGD med mindre steg

# запустим функцию, уменьшив шаг в 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()

Diagram nr 6 "Summa av kvadrerade avvikelser under stokastisk gradientnedstigning (80 tusen steg)"

Koefficienterna har förbättrats, men är fortfarande inte idealiska. Hypotetiskt kan detta korrigeras på detta sätt. Vi väljer till exempel i de senaste 1000 iterationerna värdena på koefficienterna med vilka det minsta felet gjordes. Det är sant, för detta måste vi också skriva ner värdena för själva koefficienterna. Vi kommer inte att göra detta, utan snarare uppmärksamma schemat. Det ser smidigt ut och felet verkar minska jämnt. Detta är faktiskt inte sant. Låt oss titta på de första 1000 iterationerna och jämföra dem med de senaste.

Kod för SGD-diagram (första 1000 stegen)

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()

Diagram nr 7 "Summa av kvadrerade avvikelser SGD (första 1000 stegen)"

Diagram nr 8 "Summa av kvadrerade avvikelser SGD (senaste 1000 stegen)"

Allra i början av nedstigningen observerar vi en ganska enhetlig och brant minskning av felet. I de sista iterationerna ser vi att felet går runt och runt värdet 1,475 och vid vissa tillfällen till och med är lika med detta optimala värde, men sedan går det fortfarande upp... Jag upprepar, du kan skriva ner värdena för koefficienter и , och välj sedan de där felet är minimalt. Men vi hade ett allvarligare problem: vi var tvungna att ta 80 tusen steg (se kod) för att få värden nära optimala. Och detta motsäger redan idén om att spara beräkningstid med stokastisk gradientnedstigning i förhållande till gradientnedstigning. Vad kan korrigeras och förbättras? Det är inte svårt att lägga märke till att vi i de första iterationerna med säkerhet går ner och därför bör vi lämna ett stort steg i de första iterationerna och minska steget när vi går framåt. Vi kommer inte att göra detta i den här artikeln - den är redan för lång. De som vill kan själva fundera på hur man gör detta, det är inte svårt :)

Låt oss nu utföra stokastisk gradientnedstigning med hjälp av biblioteket numpy (och låt oss inte snubbla över stenarna som vi identifierade tidigare)

Kod för Stokastisk Gradient Descent (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

Värdena visade sig vara nästan desamma som när de gick ner utan att använda numpy. Detta är dock logiskt.

Låt oss ta reda på hur lång tid nedgångar med stokastiska gradienter tog oss.

Kod för att bestämma SGD-beräkningstid (80 tusen steg)

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)

Ju längre in i skogen, desto mörkare blir molnen: återigen visar den "självskrivna" formeln det bästa resultatet. Allt detta tyder på att det måste finnas ännu mer subtila sätt att använda biblioteket numpy, vilket verkligen påskyndar beräkningsoperationer. I den här artikeln kommer vi inte att lära oss om dem. Det kommer att finnas något att tänka på på fritiden :)

sammanfatta

Innan jag sammanfattar skulle jag vilja svara på en fråga som med största sannolikhet uppstod från vår kära läsare. Varför, faktiskt, sådan "tortyr" med nedstigningar, varför behöver vi gå upp och ner för berget (för det mesta ner) för att hitta det värdefulla låglandet, om vi har i våra händer en så kraftfull och enkel anordning, i form av en analytisk lösning, som omedelbart teleporterar oss till rätt plats?

Svaret på denna fråga ligger på ytan. Nu har vi tittat på ett mycket enkelt exempel, där det sanna svaret är beror på ett tecken . Du ser inte detta ofta i livet, så låt oss föreställa oss att vi har 2, 30, 50 eller fler tecken. Låt oss lägga till tusentals eller till och med tiotusentals värden för varje attribut. I detta fall kan det hända att analyslösningen inte klarar testet och misslyckas. I sin tur kommer gradientnedstigning och dess variationer sakta men säkert att föra oss närmare målet - funktionens minimum. Och oroa dig inte för hastighet - vi kommer förmodligen att titta på sätt som gör att vi kan ställa in och reglera steglängden (det vill säga hastighet).

Och nu den faktiska korta sammanfattningen.

För det första hoppas jag att materialet som presenteras i artikeln kommer att hjälpa nybörjare "dataforskare" att förstå hur man löser enkla (och inte bara) linjära regressionsekvationer.

För det andra tittade vi på flera sätt att lösa ekvationen. Nu kan vi, beroende på situationen, välja den som är bäst lämpad för att lösa problemet.

För det tredje såg vi kraften i ytterligare inställningar, nämligen steglängden för gradientnedstigning. Denna parameter kan inte försummas. Som noterats ovan, för att minska kostnaden för beräkningar, bör steglängden ändras under nedstigningen.

För det fjärde, i vårt fall, visade "hemskrivna" funktioner de bästa tidsresultaten för beräkningar. Detta beror förmodligen på att bibliotekets kapacitet inte används mest professionellt numpy. Men hur som helst, följande slutsats antyder sig själv. Å ena sidan är det ibland värt att ifrågasätta etablerade åsikter, och å andra sidan är det inte alltid värt att komplicera allt – tvärtom, ibland är ett enklare sätt att lösa ett problem mer effektivt. Och eftersom vårt mål var att analysera tre metoder för att lösa en enkel linjär regressionsekvation, räckte användningen av "självskrivna" funktioner för oss.