🥇Lihtsa lineaarse regressiooni võrrandi lahendamine

Artiklis käsitletakse mitmeid viise lihtsa (paaris) regressioonisirge matemaatilise võrrandi määramiseks.

Kõik siin käsitletud võrrandi lahendamise meetodid põhinevad vähimruutude meetodil. Tähistame meetodeid järgmiselt:

Analüütiline lahendus
Gradiendi laskumine
Stohhastilise gradiendi laskumine

Iga sirgjoone võrrandi lahendamise meetodi jaoks pakub artikkel erinevaid funktsioone, mis jagunevad peamiselt nendeks, mis on kirjutatud ilma teeki kasutamata tuim ja need, mida kasutatakse arvutuste tegemiseks tuim. Arvatakse, et oskuslik kasutamine tuim vähendab arvutuskulusid.

Kogu artiklis antud kood on sisse kirjutatud püüton 2.7 koos Jupyteri sülearvuti. Lähtekood ja fail näidisandmetega on postitatud Github

Artikkel on rohkem suunatud nii algajatele kui ka neile, kes on juba tasapisi hakanud õppima väga laia tehisintellekti osa - masinõpet.

Materjali illustreerimiseks kasutame väga lihtsat näidet.

Näidistingimused

Meil on viis väärtust, mis iseloomustavad sõltuvust Y pärit X (tabel nr 1):

Tabel nr 1 “Näidistingimused”

Eeldame, et väärtused on aasta kuu ja — selle kuu tulud. Teisisõnu sõltub tulu aasta kuust ja - ainus märk, millest tulu sõltub.

Näide on nii ja naa, nii tulude tingimusliku sõltuvuse seisukohalt aasta kuust kui ka väärtuste arvu seisukohast - neid on väga vähe. Selline lihtsustamine võimaldab aga, nagu öeldakse, mitte alati hõlpsalt seletada materjali, mida algajad omastavad. Ja ka numbrite lihtsus võimaldab neil, kes soovivad, lahendada näide paberil ilma märkimisväärsete tööjõukuludeta.

Oletame, et näites toodud sõltuvust saab üsna hästi lähendada vormi lihtsa (paaritud) regressioonisirge matemaatilise võrrandiga:

kus on tulu laekumise kuu, — kuul vastav tulu, и on hinnangulise sirge regressioonikoefitsiendid.

Pange tähele, et koefitsient sageli nimetatakse seda hinnangulise joone kaldeks või gradiendiks; tähistab summat, mille võrra kui see muutub .

Ilmselgelt on meie ülesanne näites valida võrrandis sellised koefitsiendid и , mille puhul meie arvutatud tuluväärtuste kõrvalekalded kuude lõikes tõesetest vastustest, s.o. proovis esitatud väärtused on minimaalsed.

Vähima ruudu meetod

Vähimruutude meetodi kohaselt tuleks kõrvalekalle arvutada selle ruudustamisel. See tehnika võimaldab vältida kõrvalekallete vastastikust tühistamist, kui neil on vastupidised märgid. Näiteks kui ühel juhul on kõrvalekalle +5 (pluss viis) ja teises -5 (miinus viis), siis hälvete summa tühistab üksteist ja on 0 (null). Hälvet saab mitte ruutu panna, vaid kasutada mooduli omadust ja siis on kõik hälbed positiivsed ja akumuleeruvad. Me ei peatu sellel punktil üksikasjalikult, vaid lihtsalt näitame, et arvutuste hõlbustamiseks on tavaks kõrvalekalde ruudustamiseks.

Nii näeb välja valem, mille abil määrame hälbete (vigade) väikseima ruudusumma:

kus on tõeste vastuste (st meie arvutatud tulu) ligikaudse väärtuse funktsioon,

kas vastused on õiged (näidis toodud tulud),

on valimi indeks (hälbe kindlaksmääramise kuu number)

Diferentseerime funktsiooni, defineerime osadiferentsiaalvõrrandid ja oleme valmis liikuma edasi analüütilise lahenduse juurde. Kuid kõigepealt teeme lühikese ekskursiooni selle kohta, mis on diferentseerimine, ja tuletagem meelde tuletise geomeetrilist tähendust.

Eristumine

Diferentseerimine on funktsiooni tuletise leidmise operatsioon.

Milleks tuletist kasutatakse? Funktsiooni tuletis iseloomustab funktsiooni muutumise kiirust ja ütleb meile selle suuna. Kui tuletis antud punktis on positiivne, siis funktsioon suureneb, vastasel juhul funktsioon väheneb. Ja mida suurem on absoluuttuletise väärtus, seda suurem on funktsiooni väärtuste muutumise kiirus ja seda järsem on funktsiooni graafiku kalle.

Näiteks Descartes'i koordinaatsüsteemi tingimustes on tuletise väärtus punktis M(0,0) võrdne + 25 tähendab, et antud punktis, kui väärtust nihutatakse paremale kokkuleppelise ühiku võrra, väärtus suureneb 25 tavaühiku võrra. Graafikul näib see väärtuste üsna järsu tõusuna antud punktist.

Veel üks näide. Tuletisväärtus on võrdne -0,1 tähendab, et kui nihutatakse ühe kokkuleppelise ühiku kohta, väärtus väheneb vaid 0,1 tavaühiku võrra. Samal ajal võime funktsiooni graafikul jälgida vaevumärgatavat langust. Mäega analoogiat tuues tundub, et erinevalt eelmisest näitest, kus pidime väga järskudele tippudele ronima, laskume väga aeglaselt mööda lauget nõlva mäest alla :)

Seega pärast funktsiooni eristamist koefitsientide järgi и , defineerime esimest järku osadiferentsiaalvõrrandid. Pärast võrrandite määramist saame kahe võrrandi süsteemi, mille lahendamisel saame valida sellised koefitsientide väärtused и , mille puhul vastavate tuletiste väärtused antud punktides muutuvad väga-väga vähe ja analüütilise lahenduse puhul ei muutu üldse. Teisisõnu, leitud koefitsientide veafunktsioon jõuab miinimumini, kuna osatuletisi väärtused nendes punktides on võrdsed nulliga.

Nii et vastavalt diferentseerimisreeglitele on esimest järku osatuletisvõrrand koefitsiendi suhtes toimub järgmisel kujul:

1. järku osatuletise võrrand suhtes toimub järgmisel kujul:

Selle tulemusena saime võrrandisüsteemi, millel on üsna lihtne analüütiline lahendus:

algus{võrrand*}
alustada{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
lõpp{juhtumid}
lõpp{võrrand*}

Enne võrrandi lahendamist teeme eellaadimise, kontrollime laadimise õigsust ja vormindame andmed.

Andmete laadimine ja vormindamine

Tuleb märkida, et kuna analüütilise lahenduse ning seejärel gradiendi ja stohhastilise gradiendi laskumise jaoks kasutame koodi kahes variandis: teeki kasutades. tuim ja ilma seda kasutamata, siis vajame sobivat andmevormingut (vt koodi).

Andmete laadimise ja töötlemise kood

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

Visualiseerimine

Nüüd, kui oleme esiteks andmed laadinud, teiseks laadimise õigsuse kontrollinud ja lõpuks andmed vormindanud, viime läbi esimese visualiseerimise. Selleks sageli kasutatav meetod on paarisplaan raamatukogud Merine. Meie näites pole teeki piiratud arvu tõttu mõtet kasutada Merine. Kasutame tavaraamatukogu matplotlib ja lihtsalt vaadake hajuvusgraafikut.

Scatterplot kood

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

Tabel nr 1 “Tulu sõltuvus aasta kuust”

Analüütiline lahendus

Kasutame kõige tavalisemaid tööriistu püüton ja lahendage võrrandisüsteem:

algus{võrrand*}
alustada{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
lõpp{juhtumid}
lõpp{võrrand*}

Crameri reegli järgi leiame ülddeterminandi, samuti determinandid poolt ja , mille järel jagades determinandi arvuga ülddeterminandile - leidke koefitsient , samamoodi leiame koefitsiendi .

Analüütilise lahenduse kood

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

Saime järgmist:

Niisiis, koefitsientide väärtused on leitud, hälvete ruudu summa on kindlaks tehtud. Joonistame hajumise histogrammile sirgjoone vastavalt leitud koefitsientidele.

Regressioonijoone kood

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

Tabel nr 2 “Õiged ja arvutatud vastused”

Saate vaadata iga kuu hälbegraafikut. Meie puhul ei saa me sellest olulist praktilist väärtust, kuid rahuldame oma uudishimu selle vastu, kui hästi iseloomustab lihtne lineaarse regressiooni võrrand tulude sõltuvust aasta kuust.

Kõrvalekallete diagrammi kood

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

Diagramm nr 3 “hälbed, %”

Mitte täiuslik, kuid täitsime oma ülesande.

Kirjutame koefitsientide määramiseks funktsiooni и kasutab raamatukogu tuim, täpsemalt kirjutame kaks funktsiooni: üks kasutab pseudoinversset maatriksit (praktikas pole soovitatav, kuna protsess on arvutuslikult keeruline ja ebastabiilne), teine kasutab maatriksvõrrandit.

Analüütilise lahenduse kood (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

Võrdleme koefitsientide määramisele kuluvat aega и , vastavalt 3 esitatud meetodile.

Kood arvutusaja arvutamiseks

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)

Väikese andmehulga korral tuleb ette “isekirjutatud” funktsioon, mis leiab Crameri meetodil koefitsiendid.

Nüüd saate koefitsientide leidmiseks liikuda muude viiside juurde и .

Gradiendi laskumine

Esiteks määratleme, mis on gradient. Lihtsamalt öeldes on gradient segment, mis näitab funktsiooni maksimaalse kasvu suunda. Analoogiliselt mäe otsa ronimisega, kus kaldpinnad on seal, kus on kõige järsem tõus mäe tippu. Mäega näidet arendades peame meeles, et tegelikult vajame võimalikult kiiresti madalikule jõudmiseks kõige järsemat laskumist ehk minimaalset - kohta, kus funktsioon ei suurene ega vähene. Sel hetkel on tuletis võrdne nulliga. Seetõttu ei vaja me gradienti, vaid antigradienti. Antigradiendi leidmiseks peate gradiendi lihtsalt korrutama -1 (miinus üks).

Pöörakem tähelepanu asjaolule, et funktsioonil võib olla mitu miinimumi ja kui oleme allpool pakutud algoritmi abil laskunud ühte neist, ei suuda me leida teist miinimumi, mis võib olla leitust madalam. Lõdvestume, see ei ohusta meid! Meie puhul on tegemist ühe miinimumiga, kuna meie funktsioon graafikul on tavaline parabool. Ja nagu me kõik peaksime oma kooli matemaatikakursusest väga hästi teadma, on paraboolil ainult üks miinimum.

Pärast seda, kui saime teada, miks me vajame gradienti ja ka seda, et gradient on segment, st vektor antud koordinaatidega, mis on täpselt samad koefitsiendid и saame rakendada gradiendi laskumist.

Enne alustamist soovitan lugeda vaid paar lauset laskumisalgoritmi kohta:

Määrame pseudojuhuslikult koefitsientide koordinaadid и . Meie näites määrame nullilähedased koefitsiendid. See on tavaline praktika, kuid igal juhtumil võib olla oma praktika.
Koordinaadist lahutada punktis 1. järku osatuletise väärtus . Seega, kui tuletis on positiivne, siis funktsioon suureneb. Seega, lahutades tuletise väärtuse, liigume kasvu vastassuunas ehk laskumise suunas. Kui tuletis on negatiivne, siis funktsioon selles punktis väheneb ja tuletise väärtuse lahutamisel liigume laskumise suunas.
Teeme sarnase toimingu koordinaadiga : lahutab punktis osatuletise väärtuse .
Et mitte üle miinimumi hüpata ja süvakosmosesse lennata, on vaja seada sammu suurus laskumise suunas. Üldiselt võiks kirjutada terve artikli sellest, kuidas sammu õigesti määrata ja kuidas seda laskumisprotsessi käigus muuta, et arvutuskulusid vähendada. Kuid nüüd on meil pisut teistsugune ülesanne ja me määrame sammu suuruse teaduslikul meetodil "torkamine" või, nagu tavakeeles öeldakse, empiiriliselt.
Kui oleme antud koordinaatidest и Lahutage tuletiste väärtused, saame uued koordinaadid и . Järgmise sammu (lahutamine) teeme juba arvutatud koordinaatidest. Ja nii algab tsükkel ikka ja jälle, kuni saavutatakse vajalik konvergents.

Kõik! Nüüd oleme valmis Mariaani süviku sügavaimat kuru otsima. Alustame.

Gradiendi laskumise kood

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

Sukeldusime Mariaani süviku põhja ja leidsime sealt kõik samad koefitsientide väärtused и , mida oli täpselt oodata.

Teeme veel ühe sukeldumise, ainult seekord täidetakse meie süvameresõiduk muude tehnoloogiatega, nimelt raamatukoguga tuim.

Gradiendi laskumise kood (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

Koefitsiendi väärtused и muutumatuks.

Vaatame, kuidas muutus gradiendi laskumisel viga ehk kuidas muutus iga sammuga hälvete ruudu summa.

Ruuthälvete summade joonistamise kood

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

Graafik nr 4 “Halbade ruudu summa gradiendi laskumisel”

Graafikul näeme, et iga sammuga viga väheneb ja pärast teatud arvu iteratsioone jälgime peaaegu horisontaalset joont.

Lõpuks hindame koodi täitmisaja erinevust:

Kood gradiendi laskumise arvutamise aja määramiseks

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)

Võib-olla teeme midagi valesti, kuid see on jällegi lihtne "kodus kirjutatud" funktsioon, mis ei kasuta raamatukogu tuim ületab teeki kasutava funktsiooni arvutusaega tuim.

Kuid me ei seisa paigal, vaid liigume edasi ühe põneva lihtsa lineaarse regressioonivõrrandi lahendamise viisi uurimise poole. Saage tuttavaks!

Stohhastilise gradiendi laskumine

Stohhastilise gradiendi laskumise tööpõhimõtte kiireks mõistmiseks on parem määrata selle erinevused tavalisest gradiendi laskumisest. Meie, gradiendi laskumise korral, tuletiste võrrandites и kasutas kõigi proovis saadaolevate tunnuste ja tõeste vastuste väärtuste summasid (st kõigi и ). Stohhastilise gradiendi laskumisel ei kasuta me kõiki valimis olevaid väärtusi, vaid valime pseudojuhuslikult nn valimiindeksi ja kasutame selle väärtusi.

Näiteks kui indeks on määratud numbriga 3 (kolm), siis võtame väärtused и , siis asendame väärtused tuletisvõrranditesse ja määrame uued koordinaadid. Seejärel, pärast koordinaatide määramist, määrame pseudojuhuslikult uuesti valimiindeksi, asendame indeksile vastavad väärtused osadiferentsiaalvõrranditega ja määrame koordinaadid uuel viisil и jne. kuni konvergents muutub roheliseks. Esmapilgul võib tunduda, et see ei tööta üldse, kuid see toimib. Tõsi, tasub tähele panna, et viga ei vähene iga sammuga, kuid tendents on kindlasti olemas.

Millised on stohhastilise gradiendi laskumise eelised tavapärase laskumise ees? Kui meie valimi suurus on väga suur ja seda mõõdetakse kümnetes tuhandetes väärtustes, siis on palju lihtsam töödelda näiteks juhuslikku tuhat neist, mitte kogu valimit. Siin tuleb mängu stohhastiline gradiendi laskumine. Meie puhul me muidugi suurt erinevust ei märka.

Vaatame koodi.

Stohhastilise gradiendi laskumise kood

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

Vaatame hoolikalt koefitsiente ja tabame end küsimast: "Kuidas see nii saab?" Saime teised koefitsiendi väärtused и . Võib-olla on stohhastiline gradiendi laskumine leidnud võrrandi jaoks optimaalsemad parameetrid? Kahjuks ei. Piisab, kui vaadata hälbete ruudu summat ja näha, et koefitsientide uute väärtuste korral on viga suurem. Me ei kiirusta meeleheitesse. Koostame veamuutuse graafiku.

Stohhastilise gradiendi laskumise ruudu hälvete summa joonistamise kood

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

Graafik nr 5 “Stohhastilise gradiendi laskumise ajal hälvete ruudu summa”

Ajakava vaadates loksub kõik paika ja nüüd teeme kõik korda.

Mis juhtus? Juhtus järgmine. Kui valime kuu juhuslikult, püüab meie algoritm tulude arvutamisel viga vähendada just valitud kuu kohta. Seejärel valime teise kuu ja kordame arvutust, kuid vähendame teise valitud kuu viga. Nüüd pidage meeles, et esimesed kaks kuud erinevad oluliselt lihtsa lineaarse regressiooni võrrandi joonest. See tähendab, et kui valitakse mõni neist kahest kuust, suurendab meie algoritm vigade arvu kõigis nendes kuudes oluliselt kogu valimi puhul. Mida siis teha? Vastus on lihtne: peate vähendama laskumise sammu. Lõppude lõpuks, vähendades laskumisastet, lõpetab viga ka üles-alla "hüppamise". Õigemini, "hüppamise" viga ei peatu, kuid see ei tee seda nii kiiresti :) Vaatame üle.

Kood SGD käitamiseks väiksemate sammudega

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

Graafik nr 6 “Stohhastilise gradiendi laskumise (80 tuhat sammu) hälvete ruudu summa”

Koefitsiendid on paranenud, kuid pole endiselt ideaalsed. Hüpoteetiliselt saab seda niimoodi parandada. Valime näiteks viimase 1000 iteratsiooni koefitsientide väärtused, millega tehti minimaalne viga. Tõsi, selleks peame üles kirjutama ka koefitsientide endi väärtused. Me seda ei tee, vaid pigem pöörame tähelepanu ajakavale. See näeb välja sile ja viga väheneb ühtlaselt. Tegelikult pole see tõsi. Vaatame esimest 1000 iteratsiooni ja võrdleme neid viimasega.

SGD diagrammi kood (esimesed 1000 sammu)

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

Graafik nr 7 „SGD hälvete ruudu summa (esimesed 1000 sammu)”

Graafik nr 8 „SGD hälvete ruudu summa (viimased 1000 sammu)”

Kohe laskumise alguses täheldame vea üsna ühtlast ja järsku vähenemist. Viimastes iteratsioonides näeme, et viga käib ümber ja ümber väärtuse 1,475 ja mõnel hetkel on isegi võrdne selle optimaalse väärtusega, kuid siis see ikkagi tõuseb... Ma kordan, võite üles kirjutada koefitsiendid и ja seejärel valige need, mille puhul viga on minimaalne. Meil oli aga tõsisem probleem: pidime tegema 80 tuhat sammu (vt koodi), et saada optimaalsele lähedased väärtused. Ja see on juba vastuolus ideega säästa arvutusaega stohhastilise gradiendi laskumise ja gradiendi laskumise suhtes. Mida saab parandada ja parandada? Pole raske märgata, et esimestes iteratsioonides liigume enesekindlalt allapoole ja seetõttu peaksime esimestes iteratsioonides jätma suure sammu ja edasi liikudes sammu vähendama. Selles artiklis me seda ei tee – see on juba liiga pikk. Soovijad võivad ise mõelda, kuidas seda teha, see pole keeruline :)

Nüüd teeme teeki kasutades stohhastilise gradiendi laskumise tuim (ja ärgem komistagem kivide otsa, mille me varem tuvastasime)

Stohhastilise gradiendi laskumise kood (NumPy)

# для начала напишем функцию градиентного шага
def stoch_grad_step_numpy(vector_init, X, ind, y, l):
    x = X[ind]
    y_pred = np.dot(x,vector_init)
    err = y_pred - y[ind]
    grad_a = err
    grad_b = x[1]*err
    return vector_init - l*np.array([grad_a, grad_b])

# определим функцию стохастического градиентного спуска
def stoch_grad_descent_numpy(X, y, l=0.1, steps = 800):
    vector_init = np.array([[np.random.randint(X.shape[0])], [np.random.randint(X.shape[0])]])
    errors = []
    for i in range(steps):
        ind = np.random.randint(X.shape[0])
        new_vector = stoch_grad_step_numpy(vector_init, X, ind, y, l)
        vector_init = new_vector
        errors.append(error_square_numpy(vector_init,X,y))
    return (vector_init), (errors)

# запишем массив значений 
list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print



print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print

Väärtused osutusid peaaegu samadeks kui ilma kasutamata laskumisel tuim. See on aga loogiline.

Uurime, kui kaua stohhastilise gradiendi laskumine meil aega võttis.

Kood SGD arvutusaja määramiseks (80 tuhat sammu)

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)

Mida kaugemale metsa, seda tumedamaks muutuvad pilved: jällegi näitab “isekirjutatud” valem parimat tulemust. Kõik see viitab sellele, et raamatukogu kasutamiseks peab olema veelgi peenemaid viise tuim, mis tõesti kiirendavad arvutustoiminguid. Selles artiklis me nende kohta ei õpi. Vabal ajal on mille üle mõelda :)

Kokkuvõtteks

Enne kokkuvõtte tegemist tahaksin vastata küsimusele, mis suure tõenäosusega tekkis meie armsalt lugejalt. Milleks õigupoolest selline laskumistega “piinamine”, miks on vaja mäest üles ja alla (enamasti alla) kõndida, et leida väärtuslik madalik, kui meie käes on nii võimas ja lihtne seade, analüütilise lahenduse vormis, mis teleporteerib meid koheselt õigesse kohta?

Vastus sellele küsimusele peitub pinnal. Nüüd vaatasime väga lihtsat näidet, milles on õige vastus oleneb ühest märgist . Te ei näe seda elus sageli, nii et kujutame ette, et meil on 2, 30, 50 või enam märki. Lisame sellele iga atribuudi jaoks tuhandeid või isegi kümneid tuhandeid väärtusi. Sel juhul ei pruugi analüütiline lahendus testile vastu pidada ja ebaõnnestuda. Gradient laskumine ja selle variatsioonid viivad meid omakorda aeglaselt, kuid kindlalt lähemale eesmärgile - funktsiooni miinimumile. Ja ärge muretsege kiiruse pärast – tõenäoliselt uurime võimalusi, mis võimaldavad meil sammu pikkust (st kiirust) määrata ja reguleerida.

Ja nüüd tegelik lühikokkuvõte.

Esiteks loodan, et artiklis esitatud materjal aitab alustavatel "andmeteadlastel" mõista, kuidas lahendada lihtsaid (ja mitte ainult) lineaarse regressiooni võrrandeid.

Teiseks vaatlesime võrrandi lahendamiseks mitmeid viise. Nüüd saame sõltuvalt olukorrast valida selle, mis probleemi lahendamiseks kõige paremini sobib.

Kolmandaks nägime lisaseadete võimsust, nimelt gradiendi laskumise sammu pikkust. Seda parameetrit ei saa tähelepanuta jätta. Nagu eespool märgitud, tuleks arvutuste maksumuse vähendamiseks laskumise ajal sammu pikkust muuta.

Neljandaks näitasid meie puhul parimad ajatulemused arvutuste tegemiseks “kodukirjutatud” funktsioonid. Selle põhjuseks on ilmselt raamatukogu võimaluste mitte kõige professionaalsem kasutamine tuim. Kuid olgu nii, järgmine järeldus viitab iseenesest. Ühelt poolt tasub vahel kahtluse alla seada väljakujunenud arvamusi, teisalt ei tasu alati kõike keeruliseks ajada – vastupidi, vahel on lihtsam viis probleemi lahendamiseks tõhusam. Ja kuna meie eesmärk oli analüüsida kolme lähenemist lihtsa lineaarse regressioonivõrrandi lahendamiseks, siis piisas meie jaoks "isekirjutatud" funktsioonide kasutamisest.