🥇Zgjidhja e ekuacionit të regresionit të thjeshtë linear

Artikulli diskuton disa mënyra për të përcaktuar ekuacionin matematikor të një linje regresioni të thjeshtë (të çiftuar).

Të gjitha metodat e zgjidhjes së ekuacionit të diskutuar këtu bazohen në metodën e katrorëve më të vegjël. Le t'i shënojmë metodat si më poshtë:

Zgjidhje analitike
Zbritja me gradient
Zbritja e gradientit stokastik

Për secilën metodë të zgjidhjes së ekuacionit të një vije të drejtë, artikulli ofron funksione të ndryshme, të cilat kryesisht ndahen në ato që janë shkruar pa përdorur bibliotekën i mprehtë dhe ato që përdorin për llogaritje i mprehtë. Besohet se përdorimi i aftë i mprehtë do të zvogëlojë kostot e llogaritjes.

I gjithë kodi i dhënë në artikull është shkruar në gjuhën python 2.7 me Fletore Jupyter. Kodi burimor dhe skedari me të dhënat e mostrës janë postuar në Github

Artikulli synon më shumë si për fillestarët ashtu edhe për ata që tashmë kanë filluar gradualisht të zotërojnë studimin e një seksioni shumë të gjerë në inteligjencën artificiale - mësimin e makinerisë.

Për të ilustruar materialin, ne përdorim një shembull shumë të thjeshtë.

Shembull i kushteve

Kemi pesë vlera që karakterizojnë varësinë Y nga X (Tabela nr. 1):

Tabela nr. 1 “Shembull kushtet”

Do të supozojmë se vlerat është muaji i vitit, dhe — të ardhurat këtë muaj. Me fjalë të tjera, të ardhurat varen nga muaji i vitit, dhe - e vetmja shenjë nga e cila varen të ardhurat.

Shembulli është kështu, si nga pikëpamja e varësisë së kushtëzuar të të ardhurave nga muaji i vitit, ashtu edhe nga pikëpamja e numrit të vlerave - ka shumë pak prej tyre. Sidoqoftë, një thjeshtim i tillë do të bëjë të mundur, siç thonë ata, të shpjegojë, jo gjithmonë me lehtësi, materialin që përvetësojnë fillestarët. Dhe gjithashtu thjeshtësia e numrave do t'u lejojë atyre që dëshirojnë të zgjidhin shembullin në letër pa kosto të konsiderueshme të punës.

Le të supozojmë se varësia e dhënë në shembull mund të përafrohet mjaft mirë nga ekuacioni matematik i një linje regresioni të thjeshtë (të çiftuar) të formës:

ku është muaji në të cilin janë marrë të ardhurat, - të ardhurat që korrespondojnë me muajin, и janë koeficientët e regresionit të vijës së vlerësuar.

Vini re se koeficienti shpesh quhet pjerrësia ose gradienti i vijës së vlerësuar; paraqet shumën me të cilën kur ndryshon .

Natyrisht, detyra jonë në shembull është të zgjedhim koeficientë të tillë në ekuacion и , në të cilat devijimet e vlerave tona të llogaritura të të ardhurave sipas muajve nga përgjigjet e vërteta, d.m.th. vlerat e paraqitura në mostër do të jenë minimale.

Metoda me katrorin më të vogël

Sipas metodës së katrorëve më të vegjël, devijimi duhet të llogaritet duke e katroruar atë. Kjo teknikë ju lejon të shmangni anulimin e ndërsjellë të devijimeve nëse ato kanë shenja të kundërta. Për shembull, nëse në një rast, devijimi është +5 (plus pesë), dhe në tjetrën -5 (minus pesë), atëherë shuma e devijimeve do të anulojë njëra-tjetrën dhe do të arrijë në 0 (zero). Është e mundur të mos vendoset në katror devijimi, por të përdoret vetia e modulit dhe atëherë të gjitha devijimet do të jenë pozitive dhe do të grumbullohen. Ne nuk do të ndalemi në këtë pikë në detaje, por thjesht tregojmë se për lehtësinë e llogaritjeve, është zakon që devijimi të vendoset në katror.

Ja si duket formula me të cilën do të përcaktojmë shumën më të vogël të devijimeve në katror (gabimet):

ku është një funksion i përafrimit të përgjigjeve të vërteta (domethënë të ardhurat që kemi llogaritur),

janë përgjigjet e vërteta (të ardhurat e dhëna në mostër),

është indeksi i mostrës (numri i muajit në të cilin përcaktohet devijimi)

Le të diferencojmë funksionin, të përcaktojmë ekuacionet diferenciale të pjesshme dhe të jemi gati të kalojmë në zgjidhjen analitike. Por së pari, le të bëjmë një ekskursion të shkurtër se çfarë është diferencimi dhe të kujtojmë kuptimin gjeometrik të derivatit.

Diferencimi

Diferencimi është operacioni i gjetjes së derivatit të një funksioni.

Për çfarë përdoret derivati? Derivati i një funksioni karakterizon shpejtësinë e ndryshimit të funksionit dhe na tregon drejtimin e tij. Nëse derivati në një pikë të caktuar është pozitiv, atëherë funksioni rritet, përndryshe funksioni zvogëlohet. Dhe sa më e madhe të jetë vlera e derivatit absolut, aq më e lartë është shkalla e ndryshimit të vlerave të funksionit, si dhe sa më e madhe të jetë pjerrësia e grafikut të funksionit.

Për shembull, në kushtet e një sistemi koordinativ kartezian, vlera e derivatit në pikën M(0,0) është e barabartë me +25 do të thotë se në një pikë të caktuar, kur vlera është zhvendosur në të djathtë nga një njësi konvencionale, vlera rritet me 25 njësi konvencionale. Në grafik duket si një rritje mjaft e madhe e vlerave nga një pikë e caktuar.

Një shembull tjetër. Vlera e derivatit është e barabartë -0,1 do të thotë se kur zhvendoset për një njësi konvencionale, vlera zvogëlohet me vetëm 0,1 njësi konvencionale. Në të njëjtën kohë, në grafikun e funksionit, ne mund të vëzhgojmë një pjerrësi zbritëse mezi të dukshme. Duke nxjerrë një analogji me një mal, është sikur po zbresim shumë ngadalë një shpat të butë nga një mal, ndryshe nga shembulli i mëparshëm, ku na u desh të ngjiteshim në maja shumë të pjerrëta :)

Kështu, pas diferencimit të funksionit sipas gjasave и , ne përcaktojmë ekuacione diferenciale të pjesshme të rendit të parë. Pas përcaktimit të ekuacioneve, do të marrim një sistem prej dy ekuacionesh, duke zgjidhur të cilat do të mund të zgjedhim vlera të tilla të koeficientëve и , për të cilat vlerat e derivateve përkatëse në pikat e dhëna ndryshojnë me një sasi shumë, shumë të vogël, dhe në rastin e një zgjidhjeje analitike nuk ndryshojnë fare. Me fjalë të tjera, funksioni i gabimit në koeficientët e gjetur do të arrijë një minimum, pasi vlerat e derivateve të pjesshme në këto pika do të jenë të barabarta me zero.

Pra, sipas rregullave të diferencimit, ekuacioni derivat i pjesshëm i rendit të parë në lidhje me koeficientin do të marrë formën:

Ekuacioni derivat i pjesshëm i rendit të parë në lidhje me do të marrë formën:

Si rezultat, ne morëm një sistem ekuacionesh që ka një zgjidhje mjaft të thjeshtë analitike:

fillimi{ekuacioni*}
fillimi{rastet}
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
fund{rastet}
fundi{ekuacioni*}

Përpara se të zgjidhim ekuacionin, le të ngarkojmë paraprakisht, të kontrollojmë që ngarkimi është i saktë dhe të formatojmë të dhënat.

Ngarkimi dhe formatimi i të dhënave

Duhet të theksohet se për shkak të faktit se për zgjidhjen analitike, dhe më pas për zbritjen e gradientit dhe stokastik të gradientit, ne do të përdorim kodin në dy variacione: duke përdorur bibliotekën i mprehtë dhe pa e përdorur atë, atëherë do të na duhet formatimi i duhur i të dhënave (shih kodin).

Kodi i ngarkimit dhe përpunimit të të dhënave

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

Vizualizimi

Tani, pasi së pari të kemi ngarkuar të dhënat, së dyti, të kemi kontrolluar korrektësinë e ngarkimit dhe në fund të formatojmë të dhënat, do të kryejmë vizualizimin e parë. Metoda e përdorur shpesh për këtë është parcela çift bibliotekë I lindur në det. Në shembullin tonë, për shkak të numrit të kufizuar, nuk ka kuptim përdorimi i bibliotekës I lindur në det. Ne do të përdorim bibliotekën e zakonshme matplotlib dhe thjesht shikoni skemën e shpërndarjes.

Kodi Scatterplot

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

Grafiku nr. 1 “Varësia e të ardhurave nga muaji i vitit”

Zgjidhje analitike

Le të përdorim mjetet më të zakonshme në piton dhe zgjidhni sistemin e ekuacioneve:

fillimi{ekuacioni*}
fillimi{rastet}
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
fund{rastet}
fundi{ekuacioni*}

Sipas rregullit të Cramer-it do të gjejmë përcaktorin e përgjithshëm, si dhe përcaktorët nga dhe , pas së cilës, pjesëtimi i përcaktorit me tek përcaktorja e përgjithshme - gjeni koeficientin , në mënyrë të ngjashme gjejmë koeficientin .

Kodi i zgjidhjes analitike

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

Ja çfarë kemi:

Pra, janë gjetur vlerat e koeficientëve, është vendosur shuma e devijimeve në katror. Le të vizatojmë një vijë të drejtë në histogramin e shpërndarjes në përputhje me koeficientët e gjetur.

Kodi i linjës së regresionit

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

Grafiku nr. 2 “Përgjigje të sakta dhe të llogaritura”

Ju mund të shikoni grafikun e devijimit për çdo muaj. Në rastin tonë, nuk do të nxjerrim ndonjë vlerë praktike të rëndësishme prej tij, por do të kënaqim kureshtjen tonë se sa mirë ekuacioni i thjeshtë i regresionit linear karakterizon varësinë e të ardhurave nga muaji i vitit.

Kodi i grafikut të devijimit

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

Grafiku nr. 3 “Devijimet, %”

Jo perfekte, por ne e përfunduam detyrën tonë.

Le të shkruajmë një funksion që, për të përcaktuar koeficientët и përdor bibliotekën i mprehtë, më saktë, do të shkruajmë dy funksione: njëri duke përdorur një matricë pseudoinverse (nuk rekomandohet në praktikë, pasi procesi është llogaritarisht kompleks dhe i paqëndrueshëm), tjetri duke përdorur një ekuacion matricë.

Kodi i zgjidhjes analitike (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

Le të krahasojmë kohën e shpenzuar për përcaktimin e koeficientëve и , në përputhje me 3 metodat e paraqitura.

Kodi për llogaritjen e kohës së llogaritjes

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)

Me një sasi të vogël të dhënash, del përpara një funksion "i vetë-shkruar", i cili gjen koeficientët duke përdorur metodën e Cramer.

Tani mund të kaloni në mënyra të tjera për të gjetur koeficientët и .

Zbritja me gradient

Së pari, le të përcaktojmë se çfarë është një gradient. E thënë thjesht, gradienti është një segment që tregon drejtimin e rritjes maksimale të një funksioni. Për analogji me ngjitjen e një mali, ku fytyrat e gradientit janë aty ku është ngjitja më e pjerrët në majën e malit. Duke zhvilluar shembullin me malin, kujtojmë se në fakt na duhet zbritja më e pjerrët për të arritur sa më shpejt në ultësirë, domethënë minimumin - vendin ku funksioni nuk rritet ose zvogëlohet. Në këtë pikë derivati do të jetë i barabartë me zero. Prandaj, nuk kemi nevojë për një gradient, por një antigradient. Për të gjetur antigradientin ju duhet vetëm të shumëzoni gradientin me -1 (minus një).

Le t'i kushtojmë vëmendje faktit që një funksion mund të ketë disa minima, dhe pasi të kemi zbritur në njërën prej tyre duke përdorur algoritmin e propozuar më poshtë, nuk do të mund të gjejmë një minimum tjetër, i cili mund të jetë më i ulët se ai i gjetur. Le të pushojmë, ky nuk është një kërcënim për ne! Në rastin tonë kemi të bëjmë me një minimum të vetëm, që nga funksioni ynë në grafik është një parabolë e rregullt. Dhe siç duhet të dimë shumë mirë të gjithë nga kursi ynë i matematikës në shkollë, një parabolë ka vetëm një minimum.

Pasi zbuluam pse na duhej një gradient, dhe gjithashtu se gradienti është një segment, domethënë një vektor me koordinata të dhëna, të cilat janë saktësisht të njëjtët koeficientë и ne mund të zbatojmë zbritjen gradient.

Para fillimit, unë sugjeroj të lexoni vetëm disa fjali në lidhje me algoritmin e zbritjes:

Ne përcaktojmë në mënyrë pseudo rastësore koordinatat e koeficientëve и . Në shembullin tonë, ne do të përcaktojmë koeficientët afër zeros. Kjo është një praktikë e zakonshme, por secili rast mund të ketë praktikën e vet.
Nga koordinata zbres vlerën e derivatit të pjesshëm të rendit të parë në pikë . Pra, nëse derivati është pozitiv, atëherë funksioni rritet. Prandaj, duke zbritur vlerën e derivatit, do të lëvizim në drejtim të kundërt të rritjes, pra në drejtim të zbritjes. Nëse derivati është negativ, atëherë funksioni në këtë pikë zvogëlohet dhe duke zbritur vlerën e derivatit lëvizim në drejtim të zbritjes.
Ne kryejmë një operacion të ngjashëm me koordinatën : zbres vlerën e derivatit të pjesshëm në pikë .
Në mënyrë që të mos hidheni mbi minimumin dhe të fluturoni në hapësirën e thellë, është e nevojshme të vendosni madhësinë e hapit në drejtim të zbritjes. Në përgjithësi, mund të shkruani një artikull të tërë se si të vendosni hapin në mënyrë korrekte dhe si ta ndryshoni atë gjatë procesit të zbritjes në mënyrë që të reduktoni kostot llogaritëse. Por tani kemi një detyrë paksa të ndryshme përpara nesh, dhe do të vendosim madhësinë e hapit duke përdorur metodën shkencore të "poke" ose, siç thonë në gjuhën e zakonshme, në mënyrë empirike.
Pasi jemi nga koordinatat e dhëna и zbresim vlerat e derivateve, marrim koordinata të reja и . Ne bëjmë hapin tjetër (zbritje), tashmë nga koordinatat e llogaritura. Dhe kështu cikli fillon përsëri dhe përsëri, derisa të arrihet konvergjenca e kërkuar.

Të gjitha! Tani jemi gati të shkojmë në kërkim të grykës më të thellë të Hendekut Mariana. Le të fillojmë.

Kodi për zbritjen e gradientit

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

U zhytëm në fund të Hendekut Mariana dhe aty gjetëm të njëjtat vlera të koeficientit и , që është pikërisht ajo që pritej.

Le të bëjmë një zhytje tjetër, vetëm se këtë herë, automjeti ynë në det të thellë do të mbushet me teknologji të tjera, përkatësisht një bibliotekë i mprehtë.

Kodi për zbritjen gradient (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

Vlerat e koeficientit и e pandryshueshme.

Le të shohim se si ndryshoi gabimi gjatë zbritjes së gradientit, domethënë se si ndryshoi shuma e devijimeve në katror me çdo hap.

Kodi për vizatimin e shumave të devijimeve në katror

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

Grafiku nr. 4 “Shuma e devijimeve në katror gjatë zbritjes me gradient”

Në grafik shohim se me çdo hap gabimi zvogëlohet, dhe pas një numri të caktuar përsëritjesh vërejmë një vijë pothuajse horizontale.

Më në fund, le të vlerësojmë ndryshimin në kohën e ekzekutimit të kodit:

Kodi për të përcaktuar kohën e llogaritjes së zbritjes së gradientit

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)

Ndoshta po bëjmë diçka të gabuar, por përsëri është një funksion i thjeshtë "i shkruar në shtëpi" që nuk përdor bibliotekën i mprehtë e tejkalon kohën e llogaritjes së një funksioni duke përdorur bibliotekën i mprehtë.

Por ne nuk po qëndrojmë ende, por po shkojmë drejt studimit të një mënyre tjetër emocionuese për të zgjidhur ekuacionin e thjeshtë të regresionit linear. Takohuni!

Zbritja e gradientit stokastik

Për të kuptuar shpejt parimin e funksionimit të zbritjes stokastike të gradientit, është më mirë të përcaktohen dallimet e tij nga zbritja e zakonshme e gradientit. Ne, në rastin e zbritjes së gradientit, në ekuacionet e derivateve të и përdori shumat e vlerave të të gjitha veçorive dhe përgjigjet e vërteta të disponueshme në mostër (d.m.th., shumat e të gjitha и ). Në zbritjen e gradientit stokastik, ne nuk do të përdorim të gjitha vlerat e pranishme në mostër, por në vend të kësaj, zgjedhim në mënyrë pseudo rastësisht të ashtuquajturin indeks të mostrës dhe përdorim vlerat e tij.

Për shembull, nëse indeksi përcaktohet të jetë numri 3 (tre), atëherë marrim vlerat и , pastaj i zëvendësojmë vlerat në ekuacionet e derivateve dhe përcaktojmë koordinatat e reja. Pastaj, pasi kemi përcaktuar koordinatat, ne përsëri përcaktojmë në mënyrë pseudo rastësisht indeksin e mostrës, zëvendësojmë vlerat që korrespondojnë me indeksin në ekuacionet diferenciale të pjesshme dhe përcaktojmë koordinatat në një mënyrë të re и etj. derisa konvergjenca të marrë ngjyrë të gjelbër. Në pamje të parë, mund të duket se kjo nuk mund të funksionojë fare, por funksionon. Është e vërtetë që vlen të theksohet se gabimi nuk ulet me çdo hap, por sigurisht që ka një tendencë.

Cilat janë avantazhet e zbritjes së gradientit stokastik ndaj atij konvencional? Nëse madhësia e kampionit tonë është shumë e madhe dhe matet në dhjetëra mijëra vlera, atëherë është shumë më e lehtë të përpunosh, të themi, një mijë prej tyre të rastësishëm, sesa të gjithë kampionin. Këtu hyn në lojë zbritja e gradientit stokastik. Në rastin tonë, natyrisht, ne nuk do të vërejmë shumë ndryshim.

Le të shohim kodin.

Kodi për zbritjen e gradientit stokastik

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

Ne shikojmë me kujdes koeficientët dhe kapemi duke bërë pyetjen "Si mund të jetë kjo?" Ne morëm vlera të tjera koeficienti и . Ndoshta zbritja e gradientit stokastik ka gjetur parametra më optimalë për ekuacionin? Fatkeqësisht jo. Mjafton të shikoni shumën e devijimeve në katror dhe të shihni se me vlerat e reja të koeficientëve, gabimi është më i madh. Nuk po nxitojmë të dëshpërohemi. Le të ndërtojmë një grafik të ndryshimit të gabimit.

Kodi për vizatimin e shumës së devijimeve në katror në zbritjen e gradientit stokastik

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

Grafiku nr. 5 “Shuma e devijimeve në katror gjatë zbritjes stokastike të gradientit”

Duke parë orarin, gjithçka bie në vend dhe tani do të rregullojmë gjithçka.

Pra, çfarë ndodhi? Ndodhi si më poshtë. Kur zgjedhim rastësisht një muaj, atëherë është për muajin e zgjedhur që algoritmi ynë kërkon të reduktojë gabimin në llogaritjen e të ardhurave. Pastaj zgjedhim një muaj tjetër dhe përsërisim llogaritjen, por zvogëlojmë gabimin për muajin e dytë të zgjedhur. Tani mbani mend se dy muajt e parë devijojnë ndjeshëm nga vija e ekuacionit të thjeshtë të regresionit linear. Kjo do të thotë që kur zgjidhet ndonjë nga këta dy muaj, duke reduktuar gabimin e secilit prej tyre, algoritmi ynë rrit seriozisht gabimin për të gjithë kampionin. Pra, çfarë të bëni? Përgjigja është e thjeshtë: ju duhet të zvogëloni hapin e zbritjes. Në fund të fundit, duke zvogëluar hapin e zbritjes, gabimi gjithashtu do të ndalojë "kërcimin" lart e poshtë. Ose më mirë, gabimi "kërcim" nuk do të ndalet, por nuk do ta bëjë kaq shpejt :) Le të kontrollojmë.

Kodi për të ekzekutuar SGD me rritje më të vogla

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

Grafiku nr. 6 “Shuma e devijimeve në katror gjatë zbritjes stokastike të gradientit (80 mijë hapa)”

Koeficientët janë përmirësuar, por ende nuk janë idealë. Në mënyrë hipotetike, kjo mund të korrigjohet në këtë mënyrë. Ne zgjedhim, për shembull, në 1000 përsëritjet e fundit vlerat e koeficientëve me të cilët është bërë gabimi minimal. Vërtetë, për këtë do të duhet të shkruajmë edhe vetë vlerat e koeficientëve. Ne nuk do ta bëjmë këtë, por përkundrazi do t'i kushtojmë vëmendje orarit. Duket e qetë dhe gabimi duket se zvogëlohet në mënyrë të barabartë. Në fakt kjo nuk është e vërtetë. Le të shohim 1000 përsëritjet e para dhe t'i krahasojmë ato me të fundit.

Kodi për grafikun SGD (1000 hapat e parë)

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

Grafiku nr. 7 “Shuma e devijimeve në katror SGD (1000 hapat e parë)”

Grafiku nr. 8 “Shuma e devijimeve në katror SGD (1000 hapat e fundit)”

Në fillim të zbritjes, vërejmë një rënie mjaft uniforme dhe të pjerrët të gabimit. Në përsëritjet e fundit, shohim se gabimi shkon rreth e rrotull vlerës 1,475 dhe në disa momente madje barazohet me këtë vlerë optimale, por më pas ai ende rritet... E përsëris, mund të shkruani vlerat e koeficientët и , dhe më pas zgjidhni ato për të cilat gabimi është minimal. Sidoqoftë, ne patëm një problem më serioz: na u desh të bënim 80 mijë hapa (shih kodin) për t'i afruar vlerat optimale. Dhe kjo tashmë bie në kundërshtim me idenë e kursimit të kohës së llogaritjes me zbritjen stokastike të gradientit në krahasim me zbritjen e gradientit. Çfarë mund të korrigjohet dhe përmirësohet? Nuk është e vështirë të vërehet se në përsëritjet e para ne po zbresim me besim dhe, për rrjedhojë, duhet të lëmë një hap të madh në përsëritjet e para dhe të zvogëlojmë hapin ndërsa ecim përpara. Ne nuk do ta bëjmë këtë në këtë artikull - tashmë është shumë e gjatë. Ata që dëshirojnë mund të mendojnë vetë se si ta bëjnë këtë, nuk është e vështirë :)

Tani le të kryejmë zbritjen stokastike të gradientit duke përdorur bibliotekën i mprehtë (dhe le të mos pengohemi mbi gurët që kemi identifikuar më parë)

Kodi për zbritjen e gradientit stokastik (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

Vlerat rezultuan të jenë pothuajse të njëjta si kur zbrisnin pa përdorur i mprehtë. Megjithatë, kjo është logjike.

Le të zbulojmë se sa kohë na mori zbritjet e gradientit stokastik.

Kodi për përcaktimin e kohës së llogaritjes së SGD (80 mijë hapa)

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)

Sa më larg në pyll, aq më të errëta janë retë: përsëri, formula "vetë-shkruar" tregon rezultatin më të mirë. E gjithë kjo sugjeron se duhet të ketë mënyra edhe më delikate për të përdorur bibliotekën i mprehtë, të cilat me të vërtetë shpejtojnë operacionet llogaritëse. Në këtë artikull ne nuk do të mësojmë rreth tyre. Do të keni diçka për të menduar në kohën tuaj të lirë :)

Ne përmbledhim

Para se të përmbledh, do të doja t'i përgjigjem një pyetjeje që ka shumë të ngjarë të lindi nga lexuesi ynë i dashur. Pse, në fakt, një "torturë" e tillë me zbritje, pse duhet të ecim lart e poshtë malit (kryesisht poshtë) për të gjetur ultësirën e çmuar, nëse kemi në duart tona një pajisje kaq të fuqishme dhe të thjeshtë, në formë e një zgjidhjeje analitike, e cila na teleporton menjëherë në vendin e duhur?

Përgjigja për këtë pyetje qëndron në sipërfaqe. Tani kemi parë një shembull shumë të thjeshtë, në të cilin është përgjigja e vërtetë varet nga një shenjë . Ju nuk e shihni këtë shpesh në jetë, kështu që le të imagjinojmë se kemi 2, 30, 50 ose më shumë shenja. Le t'i shtojmë kësaj mijëra, apo edhe dhjetëra mijëra vlera për çdo atribut. Në këtë rast, zgjidhja analitike mund të mos i rezistojë provës dhe të dështojë. Nga ana tjetër, zbritja e gradientit dhe variacionet e saj ngadalë por me siguri do të na afrojnë me qëllimin - minimumin e funksionit. Dhe mos u shqetësoni për shpejtësinë - ne ndoshta do të shikojmë mënyra që do të na lejojnë të vendosim dhe rregullojmë gjatësinë e hapit (d.m.th., shpejtësinë).

Dhe tani përmbledhja aktuale e shkurtër.

Së pari, shpresoj që materiali i paraqitur në artikull do të ndihmojë në fillimin e "shkencëtarëve të të dhënave" për të kuptuar se si të zgjidhin ekuacionet e thjeshta (dhe jo vetëm) të regresionit linear.

Së dyti, ne shikuam disa mënyra për të zgjidhur ekuacionin. Tani, në varësi të situatës, ne mund të zgjedhim atë që është më e përshtatshme për të zgjidhur problemin.

Së treti, ne pamë fuqinë e cilësimeve shtesë, përkatësisht gjatësinë e hapit të zbritjes gradient. Ky parametër nuk mund të neglizhohet. Siç u përmend më lart, për të ulur koston e llogaritjeve, gjatësia e hapit duhet të ndryshohet gjatë zbritjes.

Së katërti, në rastin tonë, funksionet "të shkruara në shtëpi" treguan rezultatet më të mira kohore për llogaritjet. Kjo është ndoshta për shkak të përdorimit jo më profesional të aftësive të bibliotekës i mprehtë. Por sido që të jetë, përfundimi i mëposhtëm sugjeron vetë. Nga njëra anë, ndonjëherë ia vlen të vihet në dyshim opinionet e vendosura, dhe nga ana tjetër, nuk ia vlen gjithmonë të ndërlikosh gjithçka - përkundrazi, ndonjëherë një mënyrë më e thjeshtë për të zgjidhur një problem është më efektive. Dhe meqenëse qëllimi ynë ishte të analizonim tre qasje për zgjidhjen e një ekuacioni të thjeshtë të regresionit linear, përdorimi i funksioneve "të vetë-shkruara" ishte mjaft i mjaftueshëm për ne.