Zgjidhja e ekuacionit të regresionit linear të thjeshtë

Në këtë artikull shqyrtohen disa mënyra për të përcaktuar ekuacionin matematikor të vijës së regresionit të thjeshtë (binjak).

Të gjitha metodat e shqyrtuara këtu për zgjidhjen e ekuacionit janë të bazuara në metodën e katrorëve të vogla. Le të shënjojmë metodat si më poshtë:

• Zgjidhja analitike
• Rënia gradiente
• Rënia gradiente stohastike

Për çdo metodë zgjidhjeje të ekuacionit të drejtpërdrejtë, artikulli paraqet funksione të ndryshme, të cilat kryesisht ndahen në ato që janë shkruar pa përdorimin e bibliotekës NumPy dhe ato që përdorin NumPy. Mendohet se përdorimi i aftë i NumPy do të lejojë të zvogëlohen kostot e llogaritjeve.

I gjithë kodi i paraqitur në artikull është shkruar në gjuhën python 2.7 duke përdorur Jupyter Notebook. Kodi burimor dhe skedari me të dhënat e mostrës janë publikuar në GitHub.

Artikulli është kryesisht i orientuar si për fillestarët ashtu edhe për ata që kanë filluar të njihen me këtë seksion të gjerë në inteligjencën artificiale — të mësuarit e makinave.

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

Kushtet e shembullit

Ne kemi pesë vlera, të cilat karakterizojnë varësinë Y nga X (Tabela nr. 1):

Tabela nr. 1 "Kushtet e shembullit"

Le të supozojmë se vlerat janë muaji i vitit, ndërsa është të ardhurat në atë muaj. Në fjalë të tjera, të ardhurat varen nga muaji i vitit, ndërsa është shenja e vetme nga e cila varen të ardhurat.

Shembulli është i tillë, si nga pikëpamja e varësisë hipotetike të të ardhurave nga muaji i vitit, ashtu edhe nga numri i vlerave — ato janë shumë pak. Megjithatë, kjo thjeshtësi do të lejojë që, siç thuhet, të shpjegohet material që nuk është gjithmonë i lehtë për t'u kuptuar nga fillestarët. Gjithashtu, thjeshtësia e numrave do të lejojë ata që dëshirojnë të zgjidhin shembullin në "papir" pa shumë përpjekje.

Supozoni se varësia e dhënë në shembull mund të aproksimohet mjaft mirë me një ekuacion matematikor të vijës së regresionit të thjeshtë (binjak) të tipit:

ku është muaji kur janë arritur të ardhurat, është të ardhurat e lidhura me muajin, dhe është koeficienti i regresionit të vijës së vlerësuar.

Vlen të theksohet se koeficienti shpesh njihet si koeficienti këndor ose gradienti i vijës së vlerësuar; përfaqëson sasinë me të cilën do të ndryshojë me ndryshimin e .

E qartë është se detyra jonë në shembull është të përcaktojmë në ekuacion koeficientët e tillë dhe , të cilat do të minimizojnë devijimet e vlerave tona të llogaritura të të ardhurave sipas muajve nga përgjigjet e vërteta, pra vlerat e paraqitura në mostër.

Metoda e katrorëve të vogla

Sipas metodës së katrorëve të vogla, devijimi duhet të llogaritet duke e ngritur këtë në katror. Ky hile lejon që të shmangen shlyerjet e ndërsjella të devijimeve, në rastin kur ato kanë shenja të kundërta. Për shembull, nëse në një rast devijimi është +5 (plus pesë), dhe në tjetrin -5 (minus pesë), atëherë shuma e devijimeve do të shlyhet dhe do të jetë 0 (zero). Mund të mos e ngremë devijimin në katror, por të përdorim pronësinë e modulusit dhe atëherë të gjitha devijimet do të jenë pozitive dhe do të akumulohet. Ne nuk do të ndalemi në këtë pikë në detaje, por thjesht do të shënojmë se për lehtësinë e llogaritjeve, është pranuar që devijimi të ngrihet në katror.

Ja si duket formula me ndihmën e së cilës ne do të përcaktojmë shumën më të vogël të katrorëve të devijimeve (gabimeve):

ku është funksioni që aproksimon përgjigjet e vërteta (dmth të ardhurat e llogaritura nga ne),

është përgjigjet e vërteta (të ardhurat e ofruara në mostër),

është indeksi i mostrës (numri i muajit, në të cilin ndodh përcaktimi i devijimit)

Të diferencojmë funksionin, ta përcaktojmë ekuacionin e pjesshëm dhe të jemi gati për të kaluar në zgjidhjen analitike. Por fillimisht do të bëjmë një ekskurs të vogël mbi atë se çfarë është diferencimi dhe do të kujtojmë kuptimin gjeometrik të derivatës.

Diferencimi

Diferencimi është operacioni i gjetjes së derivatës së funksionit.

Për çfarë na nevojitet derivata? Derivata e funksionit karakterizon shpejtësinë e ndryshimit të funksionit dhe na tregon drejtimin e tij. Nëse derivata në një pikë të caktuar është pozitive, atëherë funksioni rritet; përndryshe — funksioni zvogëlohet. Sa më e madhe të jetë vlera e derivatës sipas modulusit, aq më e lartë është shpejtësia e ndryshimit të vlerave të funksionit, si dhe më ngjitës është këndi i përplasjeve të grafikut të funksionit.

Për shembull, në kushte të sistemit të koordinatave Dekart, vlera e derivatës në pikën M(0,0) e barabartë me +25 tregon se në pikën e caktuar, me zhvendosjen e vlerës dhe më tutje në të djathtë me një njësi të dhënë, vlera rritet me 25 njësi të dhëna. Në grafik, kjo duket si një kënd mjaft i pjerrët i rritjes së vlerave nga pika e caktuar.

Një tjetër shembull. Vlera e derivatit të barabartë me -0,1 do të thotë se në një zhvendosje me njësinë e kushteve, vlera zvogëlohet vetëm me 0,1 njësi të kushteve. Në të njëjtën kohë, në grafik të funksionit, mund të vëzhgojmë një kënd të lehtë në rënie. Duke bërë një analogji me malin, duket se po zbresim shumë ngadalë në një pjerrësi të butë të malit, për dallim nga shembulli i mëparshëm, ku na duhej të kapnim maja shumë të pjerrëta :)

Kështu, duke kryer diferencimin e funksionit sipër coefficientëve dhe , do të përcaktojmë ekuacionet e derivatëve pjesorë të rendit të parë. Pas përcaktimit të ekuacioneve, ne do të marrim një sistem me dy ekuacione, duke zgjidhur të cilat do të jemi në gjendje të përshtatim vlera të tilla të coefficientëve dhe , për të cilat vlerat e derivatëve përkatës në pikët e caktuara ndryshojnë në një vlerë shumë të vogël, dhe në rastin e zgjidhjes analitike nuk ndryshojnë fare. Në fjalë të tjera, funksioni i gabimit me coefficientët e gjetur do të arrijë minimumin, pasi vlerat e derivatëve pjesorë në këto pika do të jenë të barabarta me zero.

Pra, sipas rregullave të diferencimit, ekuacioni i derivatit pjesor të rendit të parë sipas coeficientit do të marrë formën:

ekuacioni i derivatit pjesor të rendit të parë sipas do të marrë formën:

Në fund, ne morëm një sistem ekuacionesh, i cili ka një zgjidhje analitike mjaft të thjeshtë:

begin{equation*}
begin{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
end{cases}
end{equation*}

Para se të zgjidhim ekuacionin, le të ngarkojmë paraprakisht, të kontrollojmë saktësinë e ngarkesës dhe të formatizojmë të dhënat.

Ngarkimi dhe formatizimi i të dhënave

Duhet të theksohet se, për shkak se për zgjidhjen analitike, dhe më pas për gradientin dhe gradientin stokastik, do të aplikojmë kodin në dy variante: me përdorimin e bibliotekës NumPy dhe pa përdorim të saj, ne do të kërkojmë formatizimin e duhur të të dhënave (shih kodin).

Kodi për ngarkimin dhe përpunimin e 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 që, në radhë të parë, kemi ngarkuar të dhënat, në të dytë, kemi verifikuar saktësinë e ngarkesës dhe më në fund kemi formatizuar të dhënat, do të realizojmë vizualizimin e parë. Shpesh për këtë përdoret metoda pairplot bibliotekës Seaborn. Në shembullin tonë, për shkak të kufizimit të numrave nuk ka kuptim të përdorim bibliotekën Seaborn. Ne do të përdorim bibliotekën e zakonshme Matplotlib dhe do të shikojmë vetëm diagramin e shpërndarjes.

Kodi i diagramit të shpërndarjes

print 'Grafiku №1 "Varësia e të ardhurave nga muaji i vitit"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()

Grafiku №1 «Varësia e të ardhurave nga muaji i vitit»

Zgjidhja analitike

Do të përdorim veglat më të zakonshme në python dhe do të zgjidhim sistemin e ekuacioneve:

begin{equation*}
begin{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
end{cases}
end{equation*}

Sipas rregullit të Cramer do të gjejmë determinantin e përgjithshëm, si dhe determinantët sipas dhe sipas , pas të cilës, duke ndarë determinantin sipas në determinantin e përgjithshëm — do të gjejmë coefficientin , po ashtu do të gjejmë coefficientin .

Kodi për zgjidhjen 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 arritur:

Pra, vlerat e coefficientëve janë gjetur, shuma e katrorëve të devijimeve është e përcaktuar. Do të vizatojmë një vijë të drejtuar në histogramin e shpërndarjes sipas coefficientëve të gjetur.

Kodi për vijën e 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 №2 «Përgjigjet e sakta dhe të llogaritura»

Mund të shikojmë grafikun e devijimeve për çdo muaj. Në rastin tonë, ndonjë vlerë praktike të rëndësishme nga ai nuk do të nxjerrim, por do të plotësojmë kuriozitetin për sa mirë, ekuacioni i regresionit të thjeshtë karakterizon varësinë e të ardhurave nga muaji i vitit.

Kodi për grafikun e devijimeve

# определим функцию для формирования массива отклонений в процентах
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 №3 «Devijimet, %»

Nuk është perfekt, por ne realizuam detyrën tonë.

Do të shkruajmë një funksion që për përcaktimin e coefficientëve dhe përdor bibliotekën NumPy, saktësisht — do të shkruajmë dy funksione: një me përdorimin e matricës pseudo-inverse (nuk rekomandohet në praktikë, pasi procesi është llogarimor dhe i paqëndrueshëm), një tjetër me përdorimin e ekuacionit matriçor.

Kodi për zgjidhjen 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

Do të krahasojmë kohën e kaluar për të përcaktuar coefficientët dhe , në përputhje me 3 mënyrat e paraqitura.

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

print ' 33[1m' + ' 33[4m' + "Koha e ekzekutimit të llogaritjeve të coefficientëve pa përdorimin e bibliotekës NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Koha e ekzekutimit të llogaritjeve të coefficientëve me përdorimin e matricës pseudo-inverse:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Koha e ekzekutimit të llogaritjeve të coefficientëve me përdorimin e ekuacionit matriçor:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Në një numër të vogël të të dhënave, përpara del funksioni "i shkruar vetë", i cili gjen coefficientët duke përdorur metodën e Cramer.

Tani tashmë mund të kalojmë në mënyra të tjera për të gjetur koeficientët dhe .

Rënia gradiente

Le të fillojmë duke përcaktuar se çfarë është gradienti. Në mënyrë të thjeshtë, gradienti është një segment që tregon drejtimin e rritjes maksimale të funksionit. Në analogji me ngjitjen e një mali, aty ku shikon gradienti, është edhe ngjitja më e kujdesshme drejt majës së malit. Duke zhvilluar shembullin me malin, kujtojmë se në të vërtetë na nevojitet zbritja më e pjerrët, për të arritur sa më shpejt në nivelin më të ulët, pra minimumin — vendi ku funksioni nuk rritet dhe nuk zvogëlohet. Në këtë vend, derivati do të jetë zero. Si rrjedhojë, na nevojitet jo gradienti, por antigradienti. Për të gjetur antigradientin duhet thjesht të shumzojmë gradientin me -1 (minus një).

Vlen të theksohet se funksioni mund të ketë disa minimale dhe, duke zbritur në një nga ato përmes algoritmit të propozuar më poshtë, nuk do të mund të gjejmë një minimum tjetër që mund të jetë më i ulët se ai që kemi gjetur. Le të qetësohemi, kjo nuk na kërcënon! Në rastin tonë, kemi të bëjmë me një minimum të vetëm, pasi funksioni ynë në grafik është një parabolë e zakonshme. Siç e dimë të gjithë nga kursi i matematikës në shkollë — parabolës i takon vetëm një minimum.

Pasi zbuluam se përse na nevojitet gradienti, dhe gjithashtu që gradienti është një segment, pra një vektor me koordinata të caktuara, ato janë pikërisht koeficientët dhe mund të realizojmë zbretjen përmes gradientit.

Para se të fillojmë, propozoj të lexojmë disa fjali mbi algoritmin e zbritjes:

• Përcaktojmë në mënyrë pseudo-rastësore koordinatat e koeficientëve dhe . Në shembullin tonë do të përcaktojmë koeficientët afër zeros. Kjo është një praktikë e zakonshme, megjithatë për çdo rast mund të jetë parashikuar një praktikë e tij.
• Nga koordinata heqim vlerën e derivatit të parë në pikën . Pra, nëse derivati është pozitiv, funksioni rritet. Si rrjedhim, duke hequr vlerën e derivatit, ne do të lëvizim në drejtimin e 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 hequr vlerën e derivatit ne lëvizim në drejtim të zbritjes.
• Këtu realizojmë një operacion të ngjashëm me koordinatën : heqim vlerën e derivatit në pikën .
• Për të mos kaluar mbi minimumin dhe për të mos u larguar në kozmos, është e nevojshme të vendosim madhësinë e hapit në drejtim të zbritjes. Në tërësi, mund të shkruajmë një artikull të tërë sesi përcaktohet më mirë hapat dhe si të ndryshojnë ata gjatë zbritjes për të ulur kostot e llogaritjeve. Por tani para nesh kemi një detyrë të ndryshme, dhe përmes metodës shkencore "të goditjes" ose siç thonë ndryshe, nëpërmjet një procesi empirike, do të përcaktojmë madhësinë e hapit.
• Pas që kemi hequr nga koordinatat e caktuara dhe vlerat e derivatëve, fitojmë koordinatat e reja dhe . Kemi një hap tjetër (heqje), tashmë nga koordinatat e llogaritura. Dhe kështu, cikli fillon përsëri e përsëri, derisa të arrijmë konvergjencën e kërkuar.

Tani jemi gati të nisim kërkimin për kanionin më të thellë të Shkretëtirës Mariana. Le të fillojmë.

Kodi për zbretjen përmes 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

Ne jemi zhytur në fundin e Shkretëtirës Mariana dhe atje gjetëm të njëjtat vlera të koeficientëve dhe , që përfaqësonte atë që pritej.

Do të realizojmë një tjetër zhytje, vetëm këtë herë, përmbajtja e aparatit tonë të thellë do të jetë teknologji të tjera, domethënë biblioteka NumPy.

Kodi për zbretjen përmes gradientit (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 koeficientëve dhe nuk ndryshojnë.

Le të shohim se si ndryshoi gabimi gjatë zbretjes përmes gradientit, në kuptimin se si ndryshoi shuma e katrorëve të devijimeve me çdo hap.

Kodi për grafikun e shumës së katrorëve të devijimeve

print 'Grafiku№4 "Shuma e katrorëve të devijimeve hap pas hapi"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Hapat (Iterimi)', size=16)
plt.ylabel('Shuma e katrorëve të devijimeve', size=16)
plt.show()

Grafiku №4 «Shuma e katrorëve të devijimeve gjatë zbritjes përmes gradientit»

Në grafik shohim se me çdo hap gabimi zvogëlohet, dhe pas disa iteracioneve vërejmë një linjë pothuajse horizontale.

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

Kodi për përcaktimin e kohës së llogaritjes së gradientit

print ' 33[1m' + ' 33[4m' + "Koha e ekzekutimit të gradientit pa përdorimin e bibliotekës 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' + "Koha e ekzekutimit të gradientit me përdorimin e bibliotekës 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, funksioni i thjeshtë 'i shkruar vetë' që nuk përdor librari NumPy kalon në kohën e ekzekutimit një funksion që përdor bibliotekën NumPy.

Por ne nuk qëndrojmë vend, por lëvizim drejt studimit të një mënyre tjetër tërheqëse për zgjidhjen e ekuacionit të regresionit linear të thjeshtë. Mirëpresim!

Rënia gradiente stohastike

Për të kuptuar më shpejt parimin e funksionimit të rënies stohastike të gradientit, është më mirë të përcaktojmë dallimet e saj nga rënia normale e gradientit. Ne, në rastin e rënies së gradientit, në ekuacionet e derivatave dhe kemi përdorur shumat e vlerave të gjitha karakteristikave dhe përgjigjeve të vërteta që ekzistojnë në mostrën (dmth shumat e të gjitha dhe ). Në rënien stohastike të gradientit ne nuk do të përdorim të gjitha vlerat që kemi në mostrë, por në vend të kësaj, do të zgjedhim në mënyrë pseudofundamentale një indeks të njohur si indeks i mostrës dhe do të përdorim vlerat e tij.

Për shembull, nëse indeksi u përcaktua si numri 3 (tre), atëherë ne marrim vlerat dhe , pastaj i vendosim vlerat në ekuacionet e derivatave dhe përcaktojmë koordinatat e reja. Më pas, pasi të kemi përcaktuar koordinatat, përsëri përcaktojmë në mënyrë pseudofundamentale indeksin e mostrës, vendosim vlerat që i përgjigjen indeksit në ekuacionet e përgjithshme të derivatave, përcaktojmë përsëri koordinatat dhe etj. deri në arrijmë konvergjencën. Me sa duket, mund të duket si mund të funksionojë kjo, por ajo funksionon. E vërteta është, se nuk e ulim gabimin me çdo hap, por tendenca pa dyshim ekziston.

Cilat janë përfitimet e rënies stohastike të gradientit në krahasim me të zakonshmen? Në rast se kemi një madhësi të mostrës shumë të madhe, që matet me dhjetëra mijëra vlera, atëherë është shumë më e lehtë të trajtojmë, për shembull, një mijë nga ato vlera, sesa të trajtojmë të gjithë mostrën. Kështu që këtu ndodh rënia stohastike e gradientit. Në rastin tonë, sigurisht, nuk do të vërejmë ndonjë ndryshim të madh.

Shikojmë kodin.

Kodi për rënien stohastike të gradientit

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

Shikojmë me vëmendje koeficientët dhe na vjen në mendje pyetja 'Si ndodhi kjo?'. Na ndodhi të kemi vlera të tjera koeficientësh dhe . Ndoshta rënia stohastike e gradientit gjeti parametra më optimalë të ekuacionit? Fatkeqësisht, jo. Mjafton të shohim shumën e katrorëve të devijimeve dhe të shohim se me vlerat e reja të koeficientëve, gabimi është më i madh. Mos u ngutni të dëshpëroheni. Do të ndërtojmë një grafik të ndryshimeve të gabimeve.

Kodi për grafikun e shumës së katrorëve të devijimeve gjatë rënies stohastike të gradientit

print 'Grafiku №5 "Shuma e katrorëve të devijimeve hap pas hapi"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Hapat (Iterimi)', size=16)
plt.ylabel('Shuma e katrorëve të devijimeve', size=16)
plt.show()

Grafiku №5 'Shuma e katrorëve të devijimeve gjatë rënies stohastike të gradientit'

Pasi e shikojmë grafikun, gjithçka vendoset në vend dhe tani do t’i bëjmë të gjitha korrigjimet.

Pra, çfarë ndodhi? Ndodhi kjo. Kur ne zgjedhim në mënyrë të rastësishme një muaj, për këtë muaj të zgjedhur algoritmi ynë përpiqet të reduktojë gabimin në llogaritjen e të ardhurave. Më pas zgjedhim një muaj tjetër dhe përsërisim llogaritjen, por tani e ulim gabimin për muajin e dytë të zgjedhur. Tani kujtojmë se dy muajt e parë shpesh devijojnë nga linja e ekuacionit të regresionit linear të thjeshtë. Kjo do të thotë se çdo herë që zgjidhet një nga këta dy muaj, ndërsa reduktojmë gabimin për secilin, algoritmi ynë rrit ndjeshëm gabimin për të gjithë mostrën. Po çfarë duhet të bëjmë? Pështu është një përgjigje e thjeshtë: duhet ta ulim hapin e rënies. Sepse nëse e ulim hapin e rënies, gabimi gjithashtu do të ndalojë ‘skakada’ lart e poshtë. Në fakt, gabimi nuk do të ndalojë ‘skakada’, por do ta bëjë këtë më ngadalë :) Do të bëjmë një provë.

Kodi për nisjen e SGD me një hap më të vogël

# запустим функцию, уменьшив шаг в 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 №6 'Shuma e katrorëve të devijimeve gjatë rënies stohastike të gradientit (80 mijë hapa)'

Vlerat e koeficientëve u përmirësuan, por ende nuk janë të shkëlqyera. Hipotetikisht, mund të rregullohen kështu. Zgjedhim, për shembull, gjatë 1000 iteracioneve të fundit vlerat e koeficientëve, për të cilat është bërë gabimi minimal. Në fakt, për këtë do na duhet të regjistrojmë edhe vetë vlerat e koeficientëve. Nuk do ta bëjmë këtë, por është më mirë të kthejmë vëmendjen në grafik. Ai duket i qetë, dhe gabimi duket se po zvogëlohet njëlloj. Në të vërtetë, nuk është kështu. Do të shikojmë 1000 të parat dhe do t'i krahasojmë me të fundit.

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

print 'Grafiku №7 "Shuma e katrorëve të devijimeve hap pas hapi. 1000 iteracionet e para"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Hapat (Iteracione)', size=16)
plt.ylabel('Shuma e katrorëve të devijimeve', size=16)
plt.show()

print 'Grafiku №7 "Shuma e katrorëve të devijimeve hap pas hapi. 1000 iteracionet e fundit"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Hapat (Iteracione)', size=16)
plt.ylabel('Shuma e katrorëve të devijimeve', size=16)
plt.show()

Grafiku №7 «Shuma e katrorëve të devijimeve SGD (1000 hapat e para)»

Grafiku №8 «Shuma e katrorëve të devijimeve SGD (1000 hapat e fundit)»

Në fillim të zbritjes ne vërejmë një zvogëlim të dukshëm dhe të shpejtë të gabimit. Në iteracionet e fundit shohim se gabimi luhatet rreth vlerës së 1,475 dhe në disa momente madje arrin këtë vlerë optimale, por pastaj përsëri rritet… Mund të theksoj se mund të regjistrojmë vlerat e koeficientëve dhe , dhe më pas të zgjedhim ato për të cilat gabimi është minimal. Megjithatë, kemi një problem më serioz: na duhet të bëjmë 80,000 hapa (shq. kod) për të marrë vlera që janë afër optimumit. Kjo, tashmë, kundërshton idenë e kursimit të kohës së llogaritjeve në krahasim me zbritjen gradient. Çfarë mund të korrigjojmë dhe përmirësojmë? Nuk është e vështirë të vëresh se në iteracionet e para po shkojmë me siguri poshtë dhe, prandaj, duhet të mbajmë një hap të madh në iteracionet e para dhe ndërsa shkojmë përpara, hapi të zvogëlohet. Nuk do ta bëjmë këtë në këtë artikull - ajo tashmë është zgjatur. Të interesuarit mund ta mendojnë vetë se si të bëjnë këtë, nuk është e vështirë 🙂

Tani do të kryejmë zbritjen gradient me rastësi, duke përdorur bibliotekën NumPy (dhe nuk do të pengohemi nga pengesat që identifikuam më parë)

Kodi për zbritjen gradient me rastësi (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 e marra ishin pothuajse të njëjta si gjatë zbritjes pa përdorimin e NumPy. Megjithatë, kjo është logjike.

Të shohim sa kohë i morën zbritjet gradient me rastësi.

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

print ' 33[1m' + ' 33[4m' +
"Koha e ekzekutimit të zbritjes gradient me rastësi pa përdorimin e bibliotekës 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' +
"Koha e ekzekutimit të zbritjes gradient me rastësi me përdorimin e bibliotekës 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 që shkojmë, aq më të errëta bëhen re: përsëri formula 'e vetë-shkruar' tregon rezultat më të mirë. Të gjitha këto na bëjnë të mendojmë se duhet të ekzistojnë mënyra edhe më të hollë për të përdorur bibliotekën NumPy, të cilat në të vërtetë e përshpejtojnë operacionet llogaritëse. Në këtë artikull ne nuk do të mësojmë për to. Do të ketë diçka për të menduar në kohën e lirë:)

Përmbledhje

Para se të përmbledhim, do të doja të përgjigjem për pyetjen që gjashtë është bërë e qartë për lexuesin tonë të dashur. Pse, në të vërtetë, këto 'mundime' me zbritjet, pse duhet të ngjitemi dhe zbresim (përgjithësisht poshtë) për të gjetur uljen e dëshiruar, nëse kemi në duar një mjet të tillë të fuqishëm dhe të thjeshtë, në formën e zgjidhjes analitike, e cila na teleportohet menjëherë në vendin e duhur?

Përgjigja në këtë pyetje është e qartë. Ne tani po shqyrtojmë një shembull shumë të thjeshtë, në të cilin përgjigjja e vërtetë varet nga një faktor . Në jetën reale një gjë e tillë nuk e takon shpesh, prandaj le të supozojmë se kemi 2, 30, 50 ose më shumë faktorë. Shto këtë mijëra, madje edhe dhjetëra mijëra vlera për çdo faktor. Në këtë rast, zgjidhja analitike mund të mos përballojë provimin dhe të dështojë. Në anën tjetër, zbritja gradient dhe variacionet e saj ngadalë, por sigurt, na afrojnë drejt qëllimit - minimumit të funksionit. Dhe sa i përket shpejtësisë, mos u shqetësoni - me siguri do të shqyrtojmë mënyra që do të na lejojnë të përcaktojmë dhe rregullojmë gjatësi të hapi (dmth shpejtësinë).

Tani pra, përmbledhja e shkurtër.

Së pari, shpresoj se materiali i paraqitur në artikull do të ndihmojë fillestarët 'data scientists' në kuptimin se si të zgjidhin ekuacionet e regresionit të thjeshtë (dhe jo vetëm).

Në të dytë, shqyrtuam disa mënyra për të zgjidhur ekuacionin. Tani, në varësi të situatës, mund të zgjedhim atë që përshtatet më së miri për zgjidhjen e detyrës së caktuar.

Në të tretë, pamë fuqinë e përmirësimeve të tjera, konkretisht gjatësi e hapi të zbritjes gradient. Ky parameter nuk mund të injorohet. Siç u theksua më parë, për qëllim të reduktimit të kostove të llogaritjeve, gjatësi e hapi duhet të ndryshohet gjatë zbritjes.

Në rastin tonë, funksionet "të shkruara vetë" treguan rezultatet më të mira të llogaritjeve. Kjo ndoshta lidhet me përdorimin jo të profesional të mundësive të bibliotekës. NumPyMegjithatë, përfundimi është ky: nga njëra anë, ndonjëherë është e nevojshme të vihen në pikëpyetje mendimet e zakonshme, ndërsa nga ana tjetër, nuk është gjithmonë e nevojshme të gjithë ta komplikohet — ndonjëherë një mënyrë më e thjeshtë për të zgjidhur një problem është më efektive. Duke pasur parasysh se qëllimi ynë ishte të shqyrtonim tre qasje për zgjidhjen e ekuacionit të regresionit linear të thjeshtë, përdorimi i funksioneve "të shkruara vetë" ishte mjaft i mjaftueshëm.

Literatura (ose diçka e tillë)

1. Regresioni linear

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

2. Metoda e katrorëve të vogla

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Derivata

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradienti

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Rënia e gradientes

habr.com/ru/post/471458

habr.com/ru/post/307312

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

6. Biblioteka 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

Burimi: habr.com