🥇 Paprastosios tiesinės regresijos lygties sprendimas

Straipsnyje aptariami keli paprastos (porinės) regresijos tiesės matematinės lygties nustatymo būdai.

Visi čia aptariami lygties sprendimo būdai yra pagrįsti mažiausių kvadratų metodu. Pažymime metodus taip:

Analitinis sprendimas
Gradiento nusileidimas
Stochastinis gradiento nusileidimas

Kiekvienam tiesės lygties sprendimo būdui straipsnyje pateikiamos įvairios funkcijos, kurios daugiausia skirstomos į tas, kurios parašytos nenaudojant bibliotekos „NumPy“ ir tie, kurie naudojami skaičiavimams „NumPy“. Manoma, kad sumaniai naudojamas „NumPy“ sumažins skaičiavimo išlaidas.

Visas straipsnyje pateiktas kodas parašytas ta kalba pitonas 2.7 naudojant „Jupyter“ užrašų knygelė. Šaltinio kodas ir failas su duomenų pavyzdžiais yra paskelbti Github

Straipsnis labiau skirtas ir pradedantiesiems, ir tiems, kurie jau pamažu pradėjo mokytis labai plačios dirbtinio intelekto dalies – mašininio mokymosi.

Norėdami iliustruoti medžiagą, naudojame labai paprastą pavyzdį.

Sąlygų pavyzdžiai

Turime penkias vertes, apibūdinančias priklausomybę Y nuo X (Lentelė Nr. 1):

Lentelė Nr. 1 „Sąlygų pavyzdžiai“

Darysime prielaidą, kad vertybės yra metų mėnuo ir – šio mėnesio pajamos. Kitaip tariant, pajamos priklauso nuo metų mėnesio ir – vienintelis ženklas, nuo kurio priklauso pajamos.

Pavyzdys yra toks, tiek sąlyginės pajamų priklausomybės nuo metų mėnesio, tiek verčių skaičiaus požiūriu - jų yra labai mažai. Tačiau toks supaprastinimas leis, kaip sakoma, ne visada lengvai paaiškinti pradedančiųjų įsisavinamą medžiagą. Be to, skaičių paprastumas leis norintiems išspręsti pavyzdį popieriuje be didelių darbo sąnaudų.

Tarkime, kad pavyzdyje pateiktą priklausomybę gana gerai galima aproksimuoti paprastos (suporuotos) formos regresijos tiesės matematine lygtimi:

kur yra mėnuo, kurį buvo gautos pajamos, — mėnesio pajamos, и yra įvertintos tiesės regresijos koeficientai.

Atkreipkite dėmesį, kad koeficientas dažnai vadinamas įvertintos linijos nuolydžiu arba nuolydžiu; reiškia sumą, kuria kai pasikeičia .

Akivaizdu, kad mūsų užduotis pavyzdyje yra pasirinkti tokius koeficientus lygtyje и , kuriame mūsų apskaičiuotų pajamų reikšmių nukrypimai pagal mėnesius nuo tikrųjų atsakymų, t.y. pavyzdyje pateiktos vertės bus minimalios.

Mažiausio kvadrato metodas

Pagal mažiausių kvadratų metodą nuokrypis turėtų būti apskaičiuojamas jį padalijus kvadratu. Ši technika leidžia išvengti abipusio nukrypimų panaikinimo, jei jie turi priešingų ženklų. Pavyzdžiui, jei vienu atveju, nuokrypis yra +5 (plius penki), o kitame -5 (minus penkis), tada nukrypimų suma panaikins vienas kitą ir sudarys 0 (nulis). Galima nuokrypio ne kvadratuoti, o panaudoti modulio savybę ir tada visi nuokrypiai bus teigiami ir kaupsis. Mes nenagrinėsime šio punkto išsamiai, o tiesiog nurodysime, kad skaičiavimų patogumui įprasta nuokrypį skaičiuoti kvadratu.

Taip atrodo formulė, su kuria nustatysime mažiausią nuokrypių (klaidų) kvadratu sumą:

kur yra tikrų atsakymų (ty mūsų apskaičiuotų pajamų) apytikslės funkcija,

ar atsakymai yra teisingi (pajamos pateiktos pavyzdyje),

yra imties indeksas (mėnesio, kurį buvo nustatytas nuokrypis, skaičius)

Išskirkime funkciją, apibrėžkime dalines diferencialines lygtis ir būkime pasiruošę pereiti prie analitinio sprendimo. Tačiau pirmiausia trumpai apžvelkime, kas yra diferenciacija, ir prisiminkime geometrinę išvestinės reikšmę.

Diferencijavimas

Diferenciacija yra funkcijos išvestinės radimo operacija.

Kam naudojamas darinys? Funkcijos išvestinė apibūdina funkcijos kitimo greitį ir nurodo jos kryptį. Jei išvestinė tam tikrame taške yra teigiama, tada funkcija didėja, priešingu atveju funkcija mažėja. Ir kuo didesnė absoliučios išvestinės reikšmė, tuo didesnis funkcijos reikšmių kitimo greitis, o kartu ir statesnis funkcijos grafiko nuolydis.

Pavyzdžiui, Dekarto koordinačių sistemos sąlygomis išvestinės reikšmė taške M(0,0) yra lygi +25 reiškia, kad tam tikrame taške, kai vertė pasislenka į dešinę sutartiniu vienetu, verte padidėja 25 įprastiniais vienetais. Grafike atrodo gana staigus verčių kilimas nuo tam tikro taško.

Kitas pavyzdys. Išvestinė vertė yra lygi -0,1 reiškia, kad kai yra perkeltas vienam sutartiniam vienetui, vertė sumažėja tik 0,1 sutartinio vieneto. Tuo pačiu metu funkcijos grafike galime stebėti vos pastebimą nuolydį žemyn. Pabrėžiant analogiją su kalnu, atrodo, kad labai lėtai leidžiamės nuo kalno švelniu šlaitu, skirtingai nei ankstesniame pavyzdyje, kur reikėjo kopti į labai stačias viršukalnes :)

Taigi, diferencijavus funkciją pagal šansus и , apibrėžiame 1 eilės dalines diferencialines lygtis. Nustačius lygtis gausime dviejų lygčių sistemą, kurią išsprendę galėsime parinkti tokias koeficientų reikšmes и , kuriai atitinkamų išvestinių vertės tam tikruose taškuose keičiasi labai, labai mažai, o analitinio sprendimo atveju visai nesikeičia. Kitaip tariant, rastų koeficientų paklaidos funkcija pasieks minimumą, nes dalinių išvestinių vertės šiuose taškuose bus lygios nuliui.

Taigi pagal diferenciacijos taisykles I eilės dalinė išvestinė lygtis koeficiento atžvilgiu bus tokia forma:

1 eilės dalinės išvestinės lygtis atžvilgiu bus tokia forma:

Dėl to gavome lygčių sistemą, kuri turi gana paprastą analitinį sprendimą:

pradėti{lygtis*}
pradėti{atvejai}
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
pabaiga{atvejai}
pabaiga{lygtis*}

Prieš spręsdami lygtį, iš anksto įkelkime, patikrinkime, ar įkėlimas teisingas, ir suformatuosime duomenis.

Duomenų įkėlimas ir formatavimas

Reikėtų pažymėti, kad dėl to, kad analitiniam sprendimui, o vėliau ir gradiento bei stochastinio gradiento nusileidimui, kodą naudosime dviem variantais: naudodami biblioteką „NumPy“ ir jo nenaudojant, tuomet mums reikės atitinkamo duomenų formatavimo (žr. kodą).

Duomenų įkėlimo ir apdorojimo kodas

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

Vizualizacija

Dabar, kai, pirma, įkelsime duomenis, antra, patikrinsime įkėlimo teisingumą ir galiausiai suformatuosime duomenis, atliksime pirmąją vizualizaciją. Tam dažnai naudojamas metodas porinis sklypas bibliotekos Seabornas. Mūsų pavyzdyje dėl riboto skaičiaus biblioteka naudotis nėra prasmės Seabornas. Naudosimės įprasta biblioteka matplotlib ir tik pažiūrėk į sklaidos diagramą.

Sklaidos kodas

print 'График №1 "Зависимость выручки от месяца года"'

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

1 diagrama „Pajamų priklausomybė nuo metų mėnesio“

Analitinis sprendimas

Naudokime dažniausiai naudojamus įrankius pitonas ir išspręskite lygčių sistemą:

pradėti{lygtis*}
pradėti{atvejai}
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
pabaiga{atvejai}
pabaiga{lygtis*}

Pagal Cramerio taisyklę rasime bendrąjį determinantą, taip pat determinantus pagal ir , po kurio determinantą dalijant iš prie bendro determinanto – raskite koeficientą , panašiai randame koeficientą .

Analitinio sprendimo kodas

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

Štai ką gavome:

Taigi, buvo rastos koeficientų reikšmės, nustatyta kvadratinių nuokrypių suma. Pagal rastus koeficientus nubrėžkime sklaidos histogramoje tiesią liniją.

Regresijos linijos kodas

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

Diagrama Nr. 2 „Teisingi ir apskaičiuoti atsakymai“

Galite peržiūrėti kiekvieno mėnesio nuokrypių grafiką. Mūsų atveju iš to negausime reikšmingos praktinės vertės, tačiau patenkinsime savo smalsumą, kaip paprasta tiesinės regresijos lygtis apibūdina pajamų priklausomybę nuo metų mėnesio.

Nukrypimų diagramos kodas

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

Grafikas Nr. 3 „Nukrypimai, %“

Ne tobula, bet savo užduotį atlikome.

Parašykime funkciją, kuri koeficientams nustatyti и naudojasi biblioteka „NumPy“, tiksliau, parašysime dvi funkcijas: vieną naudojant pseudoinversinę matricą (praktikoje nerekomenduojama, nes procesas skaičiavimo požiūriu sudėtingas ir nestabilus), kitą naudojant matricos lygtį.

Analitinio sprendimo kodas (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

Palyginkime laiką, sugaištą nustatant koeficientus и , pagal 3 pateiktus metodus.

Skaičiavimo laiko skaičiavimo kodas

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)

Turint nedidelį duomenų kiekį, išeina „savarankiškai parašyta“ funkcija, kuri randa koeficientus naudojant Cramerio metodą.

Dabar galite pereiti prie kitų būdų, kaip rasti koeficientus и .

Gradiento nusileidimas

Pirma, apibrėžkime, kas yra gradientas. Paprasčiau tariant, gradientas yra segmentas, nurodantis maksimalaus funkcijos augimo kryptį. Pagal analogiją su kopimu į kalną, kur nuolydis yra ten, kur stačiausias kopimas į kalno viršūnę. Plėtodami pavyzdį su kalnu, prisimename, kad iš tikrųjų mums reikia stačiausio nusileidimo, kad kuo greičiau pasiektume žemumą, tai yra minimumą - vietą, kur funkcija nedidėja ir nemažėja. Šiuo metu išvestinė bus lygi nuliui. Todėl mums reikia ne gradiento, o antigradiento. Norėdami rasti antigradientą, tereikia gradientą padauginti iš -1 (minus vienas).

Atkreipkime dėmesį į tai, kad funkcija gali turėti kelis minimumus, o nusileidę į vieną iš jų pagal žemiau pateiktą algoritmą, nerasime kito minimumo, kuris gali būti mažesnis už rastą. Atsipalaiduokime, tai mums negresia! Mūsų atveju mes turime reikalą su vienu minimumu, nes mūsų funkcija grafike yra taisyklinga parabolė. Ir kaip mes visi turėtume labai gerai žinoti iš mokyklinio matematikos kurso, parabolė turi tik vieną minimumą.

Po to, kai išsiaiškinome, kodėl mums reikia gradiento, taip pat kad gradientas yra atkarpa, tai yra vektorius su nurodytomis koordinatėmis, kurios yra lygiai tokie patys koeficientai и galime įgyvendinti gradiento nusileidimą.

Prieš pradedant siūlau perskaityti tik kelis sakinius apie nusileidimo algoritmą:

Pseudoatsitiktiniu būdu nustatome koeficientų koordinates и . Mūsų pavyzdyje apibrėžsime koeficientus, artimus nuliui. Tai įprasta praktika, tačiau kiekvienas atvejis gali turėti savo praktiką.
Iš koordinatės taške atimti 1 eilės dalinės išvestinės vertės vertę . Taigi, jei išvestinė yra teigiama, funkcija didėja. Todėl atėmę išvestinės vertę, judėsime priešinga augimo, tai yra, nusileidimo kryptimi. Jei išvestinė neigiama, tai funkcija šiame taške mažėja ir atėmus išvestinės reikšmę judame nusileidimo kryptimi.
Panašią operaciją atliekame su koordinate : atimkite taško dalinės išvestinės vertės vertę .
Norint neperšokti minimumo ir skristi į gilųjį kosmosą, būtina nustatyti žingsnio dydį nusileidimo kryptimi. Apskritai, galite parašyti visą straipsnį apie tai, kaip teisingai nustatyti žingsnį ir kaip jį pakeisti nusileidimo proceso metu, kad sumažintumėte skaičiavimo išlaidas. Tačiau dabar mūsų laukia kiek kitokia užduotis ir žingsnio dydį nustatysime moksliniu „poke“ metodu arba, kaip sakoma bendrine kalba, empiriškai.
Kai esame nuo duotų koordinačių и atimdami išvestinių reikšmes, gauname naujas koordinates и . Mes atliekame kitą žingsnį (atimtį), jau iš apskaičiuotų koordinačių. Ir taip ciklas prasideda vėl ir vėl, kol pasiekiama reikiama konvergencija.

Viskas! Dabar esame pasiruošę leistis ieškoti giliausio Marianų įdubos tarpeklio. Pradėkime.

Nusileidimo nuo gradiento kodas

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

Nardėme į patį Marianų įdubos dugną ir ten radome visas tas pačias koeficientų reikšmes и , būtent to ir reikėjo tikėtis.

Pasinerkime dar kartą, tik šį kartą mūsų giliavandenė transporto priemonė bus užpildyta kitomis technologijomis, būtent biblioteka „NumPy“.

Gradiento nusileidimo kodas (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

Koeficientų reikšmės и nekeičiamas.

Pažiūrėkime, kaip pakito klaida gradiento nusileidimo metu, tai yra, kaip pasikeitė kvadratinių nuokrypių suma su kiekvienu žingsniu.

Kvadratų nuokrypių sumų braižymo kodas

print 'График№4 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

4 grafikas „Nukrypimų kvadratu suma nusileidžiant nuolydžiu“

Grafike matome, kad su kiekvienu žingsniu paklaida mažėja, o po tam tikro iteracijų skaičiaus stebime beveik horizontalią liniją.

Galiausiai įvertinkime kodo vykdymo laiko skirtumą:

Kodas, skirtas nustatyti gradiento nusileidimo skaičiavimo laiką

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)

Galbūt mes darome kažką ne taip, bet tai vėlgi yra paprasta „parašyta namuose“ funkcija, kuri nenaudoja bibliotekos „NumPy“ pranoksta funkcijos, naudojančios biblioteką, skaičiavimo laiką „NumPy“.

Tačiau mes nestovime vietoje, o judame link kito įdomaus būdo, kaip išspręsti paprastą tiesinės regresijos lygtį. Susitikti!

Stochastinis gradiento nusileidimas

Norint greitai suprasti stochastinio gradiento nusileidimo veikimo principą, geriau nustatyti jo skirtumus nuo įprasto gradiento nusileidimo. Mes, gradiento nusileidimo atveju, išvestinių lygtyse и naudojo visų pavyzdyje esančių funkcijų ir teisingų atsakymų verčių sumas (tai yra visų и ). Stochastinio gradiento nusileidime nenaudosime visų imtyje esančių reikšmių, o pseudoatsitiktiniu būdu pasirinksime vadinamąjį imties indeksą ir naudosime jo reikšmes.

Pavyzdžiui, jei indeksas yra nustatytas kaip skaičius 3 (trys), tada imame reikšmes и , tada reikšmes pakeičiame į išvestines lygtis ir nustatome naujas koordinates. Tada, nustatę koordinates, vėl pseudo-atsitiktinai nustatome imties indeksą, pakeičiame indeksą atitinkančias reikšmes į dalines diferencialines lygtis ir nustatome koordinates nauju būdu. и ir tt kol konvergencija taps žalia. Iš pirmo žvilgsnio gali atrodyti, kad tai visai neveikia, bet taip. Tiesa, verta pastebėti, kad paklaida nemažėja su kiekvienu žingsniu, tačiau tendencija tikrai yra.

Kokie yra stochastinio gradiento nusileidimo pranašumai prieš įprastą? Jei mūsų imties dydis yra labai didelis ir matuojamas dešimtimis tūkstančių reikšmių, tada daug lengviau apdoroti, tarkime, atsitiktinį jų tūkstantį, o ne visą imtį. Čia atsiranda stochastinis gradiento nusileidimas. Mūsų atveju, žinoma, didelio skirtumo nepastebėsime.

Pažiūrėkime į kodą.

Stochastinio gradiento nusileidimo kodas

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

Atidžiai žiūrime į koeficientus ir užduodame sau klausimą „Kaip tai gali būti? Gavome kitas koeficientų reikšmes и . Gal stochastinio gradiento nusileidimas rado daugiau optimalių lygties parametrų? Deja, ne. Pakanka pažvelgti į kvadratinių nuokrypių sumą ir pamatyti, kad su naujomis koeficientų reikšmėmis paklaida yra didesnė. Į neviltį neskubame. Sukurkime klaidos pasikeitimo grafiką.

Stochastinio gradiento nusileidimo kvadratinių nuokrypių sumos braižymo kodas

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

Grafikas Nr. 5 „Nukrypimų kvadratu suma stochastinio gradiento nusileidimo metu“

Žiūrint į grafiką viskas stoja į savo vietas ir dabar viską sutvarkysime.

Taigi, kas atsitiko? Atsitiko taip. Kai atsitiktinai pasirenkame mėnesį, tada būtent pasirinktam mėnesiui mūsų algoritmas siekia sumažinti pajamų skaičiavimo klaidą. Tada pasirenkame kitą mėnesį ir kartojame skaičiavimą, bet sumažiname antram pasirinktam mėnesiui paklaidą. Dabar atminkite, kad pirmieji du mėnesiai labai nukrypsta nuo paprastos tiesinės regresijos lygties linijos. Tai reiškia, kad pasirinkus bet kurį iš šių dviejų mėnesių, sumažinus kiekvieno iš jų paklaidą, mūsų algoritmas labai padidina visos imties paklaidą. Taigi ką daryti? Atsakymas paprastas: reikia sumažinti nusileidimo žingsnį. Galų gale, sumažinus nusileidimo žingsnį, klaida taip pat nustos "šokinėti" aukštyn ir žemyn. Tiksliau, "šokinėja" klaida nesustos, bet tai padarys ne taip greitai :) Patikrinkim.

Kodas, skirtas SGD paleisti mažesniais žingsniais

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

Grafikas Nr. 6 „Nukrypimų kvadratu suma stochastinio gradiento nusileidimo metu (80 tūkst. žingsnių)“

Koeficientai pagerėjo, bet vis dar nėra idealūs. Hipotetiškai tai galima ištaisyti tokiu būdu. Pavyzdžiui, per pastaruosius 1000 iteracijų pasirenkame koeficientų, su kuriais buvo padaryta mažiausia klaida, reikšmes. Tiesa, tam turėsime užrašyti ir pačių koeficientų reikšmes. Mes to nedarysime, o kreipsime dėmesį į tvarkaraštį. Atrodo sklandžiai, o klaida mažėja tolygiai. Tiesą sakant, tai netiesa. Pažvelkime į pirmuosius 1000 pakartojimų ir palyginkime juos su paskutiniaisiais.

SGD diagramos kodas (pirmieji 1000 žingsnių)

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

Grafikas Nr. 7 „Nukrypimų kvadratu suma SGD (pirmieji 1000 žingsnių)“

Grafikas Nr. 8 „Nukrypimų kvadratu suma SGD (pastarieji 1000 žingsnių)“

Pačioje nusileidimo pradžioje pastebime gana vienodą ir staigų paklaidos sumažėjimą. Paskutinėse iteracijose matome, kad klaida svyruoja apie 1,475 ir kai kuriais momentais netgi prilygsta šiai optimaliai vertei, bet tada ji vis tiek didėja... Pasikartosiu, galite užsirašyti koeficientai и , tada pasirinkite tuos, kurių klaida yra minimali. Tačiau turėjome rimtesnę problemą: turėjome nueiti 80 tūkstančių žingsnių (žr. kodą), kad vertės būtų artimos optimalioms. Ir tai jau prieštarauja idėjai sutaupyti skaičiavimo laiką taikant stochastinį gradiento nusileidimą, palyginti su gradiento nusileidimu. Ką galima pataisyti ir patobulinti? Nesunku pastebėti, kad pirmose iteracijose užtikrintai leidžiamės žemyn, todėl pirmosiose iteracijose turėtume palikti didelį žingsnį, o eidami į priekį žingsnį sumažinti. Šiame straipsnyje to nedarysime – jis jau per ilgas. Norintys gali patys pagalvoti, kaip tai padaryti, nesunku :)

Dabar atlikkime stochastinį gradiento nusileidimą naudodami biblioteką „NumPy“ (ir nesuklupkime už akmenų, kuriuos nustatėme anksčiau)

Stochastinio gradiento nusileidimo kodas (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

Vertės pasirodė beveik tokios pačios, kaip nusileidžiant nenaudojant „NumPy“. Tačiau tai logiška.

Sužinokime, kiek laiko užtruko stochastinio gradiento nusileidimai.

SGD skaičiavimo laiko nustatymo kodas (80 tūkst. žingsnių)

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)

Kuo toliau į mišką, tuo debesys tamsėja: vėlgi „savarankiškai parašyta“ formulė rodo geriausią rezultatą. Visa tai rodo, kad turi būti dar subtilesni bibliotekos naudojimo būdai „NumPy“, kurios tikrai pagreitina skaičiavimo operacijas. Šiame straipsnyje mes apie juos nesužinosime. Laisvalaikiu bus apie ką pamąstyti :)

Mes apibendrinsime

Prieš apibendrindamas norėčiau atsakyti į klausimą, kuris greičiausiai kilo iš mūsų brangaus skaitytojo. Kodėl iš tikrųjų tokie „kankinimai“ nusileidimais, kodėl mums reikia vaikščioti aukštyn ir žemyn nuo kalno (dažniausiai žemyn), kad surastume brangią žemumą, jei rankose turime tokį galingą ir paprastą prietaisą. analitinio sprendimo forma, kuri akimirksniu nuteleportuoja mus į tinkamą vietą?

Atsakymas į šį klausimą slypi paviršiuje. Dabar pažvelgėme į labai paprastą pavyzdį, kuriame yra tikrasis atsakymas priklauso nuo vieno ženklo . Gyvenime tai nedažnai matote, todėl įsivaizduokime, kad turime 2, 30, 50 ar daugiau ženklų. Pridėkime prie to tūkstančius ar net dešimtis tūkstančių kiekvieno atributo verčių. Tokiu atveju analitinis tirpalas gali neatlaikyti bandymo ir nepavykti. Savo ruožtu gradiento nusileidimas ir jo variacijos lėtai, bet užtikrintai priartins mus prie tikslo – funkcijos minimumo. Ir nesijaudinkite dėl greičio – tikriausiai ieškosime būdų, kurie leis mums nustatyti ir reguliuoti žingsnio ilgį (ty greitį).

O dabar tikra trumpa santrauka.

Pirma, tikiuosi, kad straipsnyje pateikta medžiaga padės pradedantiesiems „duomenų mokslininkams“ suprasti, kaip išspręsti paprastas (ir ne tik) tiesinės regresijos lygtis.

Antra, pažvelgėme į kelis lygties sprendimo būdus. Dabar, atsižvelgdami į situaciją, galime pasirinkti tą, kuris geriausiai tinka problemai išspręsti.

Trečia, pamatėme papildomų nustatymų galią, būtent gradiento nusileidimo žingsnio ilgį. Šio parametro negalima nepaisyti. Kaip minėta aukščiau, norint sumažinti skaičiavimų išlaidas, nusileidimo metu žingsnio ilgį reikia keisti.

Ketvirta, mūsų atveju geriausius skaičiavimo laiko rezultatus parodė „namuose parašytos“ funkcijos. Greičiausiai taip yra dėl ne profesionaliausio bibliotekos galimybių panaudojimo „NumPy“. Bet kaip ten bebūtų, tokia išvada leidžia suprasti save. Viena vertus, kartais verta suabejoti nusistovėjusia nuomone, kita vertus, ne visada verta viską komplikuoti – atvirkščiai, kartais paprastesnis problemos sprendimo būdas yra efektyvesnis. Ir kadangi mūsų tikslas buvo išanalizuoti tris paprastos tiesinės regresijos lygties sprendimo būdus, mums visiškai pakako „savarankiškai parašytų“ funkcijų.