Lahendame lihtsa lineaarse regressiooni võrrandi

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

Siin käsitletud kõik lahendusmeetodid põhinevad väikseimate ruutude meetodil. Määratleme meetodid järgmiselt:

• Analüütiline lahendus
• Gradientne laskumine
• Stohhastiline gradientne laskumine

Iga jooksu lahenduseks, artiklis on toodud erinevad funktsioonid, mis peamiselt jagunevad nende järgi, mis on kirjutatud ilma raamatukogude kasutamiseta NumPy ja nende, mis kasutavad arvutuste tegemiseks NumPy. Arvatakse, et oskuslik kasutamine NumPy aitab vähendada arvutuskulusid.

Kogu artiklis esitatud kood on kirjutatud keeles python 2.7 kasutades Jupyter Notebook. Algne kood ja andmefail on üles pandud GitHubis

Artikkel on suunatud rohkem nii algajatele kui ka neile, kes on juba natuke hakanud õppima väga ulatuslikku teemat tehisintellekti valdkonnas — masinõppe.

Kasutame materjali illustreerimiseks väga lihtsat näidet.

Näite tingimused

Meil on viis väärtust, mis iseloomustavad sõltuvust Y alates X (Tabel nr 1):

Tabel nr 1 „Näite tingimused“

Oletame, et väärtused — see on aasta kuu, ja — käive sellel kuul. Teisisõnu, käive sõltub aasta kuust, ja — ainus tunnus, millest käive sõltub.

Näide on kesine, nii tingimusliku sõltuvuse osas käibest aasta kuust kui ka väärtuste arvu osas — neid on väga vähe. Siiski võimaldab see lihtsustamine, nagu öeldakse, kahe sõrme vahelt selgitada, et mitte alati algajatele kergesti omandatav teema. Samuti võimaldavad numbrite lihtsus soovijatel lahendada näidet „paberil“ ilma suuremate tööjõukuludeta.

Eeldame, et näites toodud sõltuvust saab piisavalt hästi lähendada lihtsa (paarilise) regressioonijoone matemaatilise võrrandiga järgmisel kujul:

kus — see on kuu, mille jooksul käive saadud, — kuukäive, ja — regressiooni koefitsientide hinnatud joone kohta.

Tuleb märkida, et koefitsienti tavaliselt nimetatakse hinnatud joone kaldenurgaks või gradientideks; see kujutab endast suurust, mille võrra muutub muutmise korral .

Ilmselt on meie ülesanne näites valida sellised koefitsiendid ja , mille korral meie arvutatud tulude kõrvalekalded kuude jooksul tõelistest vastustest, st. andmetes esitatud väärtustest, on minimaalset.

Vähenenud ruutude meetod

Vastavalt vähenenud ruutude meetodile tuleks kõrvalekalle arvutada, tõstes selle ruutu. Selline lähenemine väldib kõrvalekallete vastastikust üksteise tühistamist, juhul kui need omavad vastupidiseid märke. Näiteks, kui ühes juhul on kõrvalekalle +5 (pluss viis) ja teises -5 (miinus viis), siis kõrvalekallete summa tühistab üksteist ja on 0 (null). Ei ole hädavajalik tõsta kõrvalekallet ruutu, vaid võib kasutada mooduli omadust, siis kõik meie kõrvalekalded on positiivsed ja kogunevad. Me ei peatugu sellel hetkel põhjalikult, vaid lihtsalt märkime, et arvutuste mugavuse huvides on tavaks tõsta kõrvalekalle ruutu.

Nii näeb välja valem, mille abil me määrame minimaalsete ruutude summa kõrvalekaldeid (vigu):

kus — see on funktsioon, mis läheneb tõelistele vastustele (st meie arvutatud tulud),

— need on tõelised vastused (andmetes esitatud tulud),

— see on valimi indeks (kuu number, milles tehakse kõrvalekalde määramine)

Teeme funktsiooni diferentsieerimise, määrame osade tuletiste võrrandid ja oleme valmis minema analüütilisele lahendusele. Kuid alustuseks teeme väikese ekskursiooni, mis on diferentsieerimine ja meenutame tuletise geomeetrilist mõtet.

Diferentsieerimine

Diferentsieerimine on operatsioon, mis määrab funktsiooni tuletise.

Miks on tuletis vajalik? Funktsiooni tuletis iseloomustab funktsiooni muutumise kiirus ja näitab meile selle suunda. Kui tuletis on antud punktis positiivne, siis funktsioon kasvab, muidu - funktsioon langeb. Ja mida suurem on tuletise väärtus moodulis, seda suurem on funktsiooni väärtuste muutumise kiirus ning järsem on funktsiooni graafiku kaldenurk.

Näiteks, dekartaalses koordinaatsüsteemis, tuletise väärtus punktis M(0,0) on +25 tähendab, et antud punktis, väärtuse nihutamise korral paremale tingimuslikku ühikut, väärtus kasvab 25 tingimuslikku ühikut. Graham nad näevad välja nagu üsna järsk väärtuste tõusunurk antud punktist.

Teine näide. Tuletise väärtus on võrdne -0,1 tähendab, et nihke korral ühe tingimusliku ühiku võrra, väärtus väheneb vaid 0,1 tingimuslikku ühikut. Samal ajal näeme funktsiooni graafikul vaevu märkamatut allapoole kallakut. Tehes analoogia mäega, on see nagu me aeglaselt laskumegi mäe lohku, erinevalt eelmisest näitest, kus pidime ronima väga järskude tippude poole :)

Seega, diferentseerides funktsiooni koefitsientide järgi ja , määrame esmakordse osatuletise võrrandid. Pärast võrrandite määratlemist saame süsteemi kahest võrrandist, mille lahendamiseks suudame leida sellised koefitsiendid ja , mille korral vastavate osatuletiste väärtused antud punktides muutuvad väga ja väga vähe, ja analüütilise lahenduse korral ei muutu nad üldse. Teisisõnu, vigafunktsioon leidmiseks saavutame minimaalse väärtuse, kuna nende punktide osatuletiste väärtused on null.

Nii, diferentseerimise reeglite järgi, esmakordse osatuletise võrrand koefitsiendist võtab järgmist vormi:

esmakordse osatuletise võrrand võtab järgmist vormi:

Kokkuvõttes saime võrrandite süsteemi, millel on üsna lihtne analüütiline lahendus:

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*}

Enne võrrandi lahendamist laadime eelnevalt, kontrollime laadimise õigsust ja vormindame andmed.

Andmete laadimine ja vormindamine

Tuleb märkida, et seoses sellega, et analüütiliseks lahenduseks ja hiljem gradientse ning stohhastilise gradientse langemise jaoks kasutame koodi kahes variandis: koos NumPy ja ilma selle kasutamiseta, vajame vastavat andmete vormindamist (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, pärast seda, kui oleme kõigepealt andmed laadinud, teiseks kontrollinud laadimise õigsust ja lõpuks andmed vormindanud, teeme esimese visualiseerimise. Selleks kasutatakse tihti meetodit pairplot teeki Seaborn. Meie näites, arvestades numbrite piiranguid, pole mõtet raamatukogu kasutada Seaborn. Kasutame tavalist raamatukogu Matplotlib ja vaatame ainult hajusdiagrammi.

Hajusdiagrammi kood

print 'Graafik №1 "Aasta kuu müügi sõltuvus"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Kuu$', size=16)
plt.ylabel('$Müük$', size=16)
plt.show()

Graafik №1 «Aasta kuu müügi sõltuvus»

Analüütiline lahendus

Kasutame kõige tavalisemaid tööriistu python ja lahendame võrrandite süsteemi:

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*}

Krämer'i reegli kohaselt leiame üldise määratleja ning ka määratlejad ja , seejärel, jagades määratleja üldise määratlejaga — leiame koefitsiendi , sarnaselt 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)

Nii, et koefitsientide väärtused on leitud, ruutsete kõrvalekalde summad on kindlaks määratud. Joonistame hajutusdiagrammile sirge vastavalt leitud koefitsientidele.

Regressiooni joone kood

Graafik №2 «Õiged ja arvutatud vastused»

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

Saame vaadata iga kuu kõrvalekalde graafikut. Meie puhul ei saa me sealt olulist praktilist väärtust, kuid rahuldame uudishimu, kui hästi lihtne lineaarne regressioonivõrrand iseloomustab müügi sõltuvust aasta kuust.

Kõrvalekalde graafiku kood

Graafik №3 «Kõrvalekalded, %»

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

Ei ole ideaalne, kuid meie ülesanne on täidetud.

Kirjutame funktsiooni, mis koefitsientide määramiseks

kasutab raamatukogu ja , täpsemalt — kirjutame kaks funktsiooni: üks pseudo-ristkülikulise maatriksi kasutamisega (mida ei soovitata praktiliselt, kuna see on arvutuslikult keeruline ja ebatavaline), teine maatrikuvõrrandi abil. NumPyAnalüütilise lahenduse kood (NumPy)

Võrdleme aega, mis kulus koefitsientide määramiseks

# для начала добавим столбец с не изменяющимся значением в 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

, vastavalt kolmele esitatud viisile. ja Arvutuste aja kood

Arvutuste aja määramise kood

print ' 33[1m' + ' 33[4m' + "Aegud koefitsentide arvutamise aeg ilma NumPy raamatukoguta:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Aegud koefitsentide arvutamiseks, kasutades pseudoinverse maatriksi:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Aegud koefitsentide arvutamiseks, kasutades maatrikvite võrrandit:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Väikese andmemahu puhul on ees teatud „enda kirjutatud“ funktsioon, mis leiab koefitsendid Krameri meetodi abil.

Nüüd saame liikuda teistele koefitsentide leidmise viisidele. ja .

Gradientne laskumine

Esialgu määratleme, mis on gradient. Lihtsalt öeldes, gradient on joon, mis näitab funktsiooni maksimaalse kasvu suunda. Sarnaselt mäkke tõusule, kuhu gradient osutab, seal on kõige järsem tõus mäe tipule. Jätkates mägi näidet, tuletame meelde, et tegelikult vajame kõige järsemat langust, et võimalikult kiiresti madalale jõuda, st miinimumini — kohani, kus funktsioon ei tõuse ega lange. Selles kohas on tuletis null. Seega vajame me mitte gradienti, vaid antigradiendi. Antigradiendi leidmiseks peab lihtsalt gradienti korrutama -1 (miinus üks).

Tuleb märkida, et funktsioonil võib olla mitu miinimumi ning kui me langeme ühte neist allapoole pakutud algoritmi järgi, ei suuda me leida teist miinimumi, mis võib olla madalam leitud miinimumist. Rahunege, see ei ohusta meid! Meie juhul on meil tegemist ainulaadse miinimumiga, kuna meie funktsioon graafikul esindab tavalist parabooli. Kuidas me kõik peaksime olema hästi teadlikud kooli matemaatika kursusest — paraboolil on ainult üks miinimum.

Pärast seda, kui oleme välja selgitanud, miks me vajame gradienti, ja samuti, et gradient on joon, st vektor, millel on kindlad koordinaadid, mis on just need koefitsendid ja saame me rakendada gradientide langust.

Enne käivitamist soovitan lugeda vaid mõne lause allakäigu algoritmist:

• Määrame koefitsentide koordinaadid pseudojuhuslikult. ja . Meie näites määrame koefitsiente nulli läheduses. See on tavaline praktika, kuid igaks juhuks võib olla ette nähtud oma praktika.
• Koordinaadist lahutame esimese järgu osatuletise väärtuse punktis . Nii et kui tuletis on positiivne, siis funktsioon kasvab. Seega, lahutades tuletise väärtuse, liigume kasvu vastassuunda, st languse suunda. Kui tuletis on negatiivne, siis funktsioon väheneb ja lahutades tuletise väärtuse liigume languse suunda.
• Teeme sarnast toimingut koordinaadiga : lahutame osatuletise väärtuse punktis .
• Kuna ei ületaks miinimumpunkti ja ei lendaks kaugesse kosmosesse, on vaja seadistada sammu suurus languse suunas. Üldiselt võiks kirjutada terve artikli selle kohta, kuidas oskuslikult sammu seada ja kuidas seda langemise protsessis muuta, et vähendada arvutuskulusid. Kuid praegu seisame silmitsi veidi teistsuguse ülesandega ja teadusliku meetodi "tunne" ehk nagu rahvasuu ütleb, empiirilise meetodi abil määrame sammu suuruse.
• Pärast seda, kui oleme antud koordinaatidest ja lahutanud tuletiste väärtused, saame uued koordinaadid ja . Teeme järgmise sammu (lahutamine), juba arvutatud koordinaatidest. Ja nii tsükkel käivitatakse uuesti ja uuesti, kuni saavutatakse vajalik konvergentsus.

Kõik! Nüüd oleme valmis minema otsima kõige sügavamate süvikute Mariannide kraavi. Alustame.

Gradientlanguse 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 Mariannide kraavi põhja ja seal avastasime kõik need sama koefitsendi väärtused ja , nagu tegelikult oodati.

Teeme veel ühe sukeldumise, kuid seekord on meie süvemeretankuri sisuks teised tehnoloogiad, nimelt raamatukogu NumPy.

Gradientlanguse 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 ja on muutumatud.

Vaadakem, kuidas viga muutus gradientlanguses, st kuidas kõrvalekalde ruutude summa muutus iga sammu jooksul.

Graafiku kood kõrvalekalde ruutude summast

print 'Graafik№4 "K kõrvalekalde ruutude summa samm-sammult"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Sammud (Itaeratsioon)', size=16)
plt.ylabel('Kvideeritud kõrvalekalded', size=16)
plt.show()

Graafik nr 4 „Nüanside ruutude summa gradientide langemisel“

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

Lõpuks hindame koodi täitmise aja erinevust:

Kood gradientide langemise arvutamise aja määramiseks

print ' 33[1m' + ' 33[4m' + "Gradientide langemise täitmise aeg ilma NumPy teeki kasutamata:" + ' 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' + "Gradientide langemise täitmise aeg NumPy teegi kasutamisega:" + ' 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 taas see lihtne „oma” funktsioon, mis ei kasuta teeki, NumPy ületab arvutuste täitmise ajal teeki kasutava funktsiooni. NumPy.

Kuid me ei seisa paigal, vaid liigume edasi, et uurida veel ühte põnevat viisi lihtsa lineaarse regressioonivõrrandi lahendamiseks. Tere tulemast!

Stohhastiline gradientne laskumine

Et paremini mõista stohhastilise gradientide langemise toimimise põhimõtet, on parem määrata selle erinevused tavalisest gradientide langemisest. Gradientide langemise puhul kasutasime me tuletiste valemites ja kõigi tunnuste ja tõeliste vastuste summasid, mis on valimis olemas (ehk kõikide summasid ja ). Stohhastilises gradientide langemises ei kasuta me kõikide väärtusi, mis on valimis olemas, vaid valime nii-öelda pseudojuhuslikult valimi indeksi ja kasutame selle väärtusi.

Näiteks, kui indeks määratakse numbriga 3 (kolm), siis võtame väärtused ja , seejärel asendame väärtused tuletiste valemites ja määrame uued koordinaadid. Seejärel, määrates koordinaadid, määrame jälle pseudojuhuslikult valimi indeksi, asendame väärtused, mis vastavad indeksile, tuletiste valemites, ja määrame uued koordinaadid ja jne. kuni konvergentsi hõbeduseni. Esmapilgul võib tunduda, et kuidas see kõik toimib, kuid see töötab. Tõsi, tuleb mainida, et mitte iga sammuga ei vähene viga, kuid suundumust on kindlasti märgata.

Millised on stohhastilise gradientide langemise eelised tavalise ees? Kui meie valimi suurus on väga suur ja mõõdetud kümnetes tuhandetes väärtustes, on oluliselt lihtsam töödelda näiteks juhuslikku tuhat neist, kui tervet valimit. Just siis käivitub stohhastiline gradientide langemine. Meie juhul me erilist erinevust ei märka.

Vaatame koodi.

Kood stohhastilise gradientide langemise jaoks

# определим функцию стох.град.шага
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 küsime endalt: "Kuidas nii?" Meil on saanud teised koefitsiendi väärtused. ja . Võib-olla leidis stohhastiline gradientide langemine optimaalsemad võrrandi parameetrid? Kahjuks mitte. Piisab, kui vaatame ruutude summat ja näeme, et uute koefitsientide väärtustega on viga suurem. Ärge kiirustage meeleheitele. Joonistame vea muutuse graafiku.

Kood ruutude summa graafiku jaoks stohhastilise gradientide langemise puhul

print 'Graafik nr 5 "Ruutude summa samm-sammult"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Sammud (Iteratsioon)', size=16)
plt.ylabel('Ruutude summa', size=16)
plt.show()

Graafik nr 5 «Ruutude summa stohhastilise gradientide langemise puhul»

Vaadates graafikut, hakkab kõik oma kohale tagasi tulema ja nüüd me parandame kõik.

Nii, mis juhtus? Juhtus see, et kui me valime juhuslikult kuu, siis just selle kuu jaoks püüab meie algoritm vähendada tulude arvutamise viga. Siis valime teise kuu ja kordame arvutust, kuid nüüd vähendame viga teise valitud kuu jaoks. Ja nüüd meenutame, et meie esimesed kaks kuud kaldusid oluliselt kõrvalekalduma lihtsast lineaarregressiooni võrrandi realt. See tähendab, et kui valitakse ükskõik milline neist kahest kuust, siis vähendades igaühe viga, suurendab meie algoritm tõsiselt vea suurust kogu valimi puhul. Mida siis teha? Vastus on lihtne: tuleb vähendada langemise sammu. Vähendades langemise sammu, lõpetab viga ka „hüppamise“ ühte või teise suunda. Tõsi, viga ei lõpeta hüppamist, kuid see ei juhtu enam nii kiiresti :) Kontrollime.

Kood SGD käivitamiseks väiksema sammuga

# запустим функцию, уменьшив шаг в 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 «Ruutude summa stohhastilise gradientide langemise puhul (80 tuhat sammu)»

Kordite väärtused on paranenud, kuid siiski mitte ideaalsed. Hüpotetoolselt saaksime seda parandada järgmiselt. Valime näiteks viimase 1000 iteratsiooni kordite väärtused, millega on tehtud minimaalne viga. Tõsi, selleks peame me salvestama ka kordite väärtused ise. Me ei tee seda, vaid pöörame tähelepanu graafikule. See näeb välja sujuv ja viga näib vähenevat ühtlaselt. Tegelikult ei ole see nii. Vaatame esimesed 1000 iteratsiooni ja võrdleme neid viimastega.

SGD graafiku kood (esimesed 1000 sammu)

print 'Graafik №7 "Kvardide summade kõrvalekalded samm-sammult. Esimesed 1000 iteratsiooni"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Sammud (Iteratsioon)', size=16)
plt.ylabel('Kvardide summade kõrvalekalded', size=16)
plt.show()

print 'Graafik №7 "Kvardide summade kõrvalekalded samm-sammult. Viimased 1000 iteratsiooni"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Sammud (Iteratsioon)', size=16)
plt.ylabel('Kvardide summade kõrvalekalded', size=16)
plt.show()

Graafik №7 «Kvardide summade kõrvalekalded SGD (esimesed 1000 sammu)»

Graafik №8 «Kvardide summade kõrvalekalded SGD (viimased 1000 sammu)»

Alguses jälgime me piisavalt ühtlast ja järsku vea vähenemist. Viimastel iteratsioonidel näeme, et viga liigub umbes 1,475 väärtuse ümber ja mõnel hetkel isegi võrdub selle optimaalse väärtusega, kuid see tõuseb ikkagi üles... Kordan, et kordite väärtuseid võiks salvestada ja , ja seejärel valida need, millega viga on minimaalne. Siiski tekkis meil tõsisem probleem: me pidime tegema 80,000 sammu (vt koodi), et saada väärtusi, mis on optimaalsele lähedased. See aga on juba vastuolus ideeaga, et arvestuste aja säästmine on stohhastilisest gradientide langemisest võrreldes gradientide langemisega. Mida saaks parandada ja täiustada? Ei ole raske märgata, et esimestel iteratsioonidel liigume kindlalt allapoole ja seega peaksime esimestel iteratsioonidel säilitama suure sammu ning järk-järgult sammu vähendama. Me ei tee seda käesolevas artiklis — see on juba piisavalt veninud. Soovijad saavad ise mõelda, kuidas seda teha, see ei ole keeruline 🙂

Nüüd teeme stohhastilise gradientide langemise, kasutades teeki NumPy (ja ei komistame kividele, mille me varem avastasime)

Stohhastilise gradientlanguse 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

Saadud väärtused olid peaaegu samad kui languse korral ilma kasutamata NumPy. Kuid see on loogiline.

Uurime, kui kaua võttis meil aega stohhastiline gradientlangus.

Stohhastilise gradientlanguse arvutamise aja määramise kood (80 tuhat sammu)

print ' 33[1m' + ' 33[4m' +
"Stohhastilise gradientlanguse teostamise aeg ilma NumPy raamatukogu kasutamata:"
+ ' 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' +
"Stohhastilise gradientlanguse teostamise aeg NumPy raamatukogu kasutades:"
+ ' 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 pilved muutuvad: taas näitab „oma valmistatud” valem parimat tulemust. Kõik see paneb mõtlema, et peaks olema olemas veelgi peenemaid viise raamatukogu kasutamiseks NumPy, mis tõeliselt kiirendavad arvutusoperatsioone. Selles artiklis me neist juba ei kuule. Mõtle seda oma vabalt valitud ajal :)

Kokkuvõtteks

Enne kokkuvõtte tegemist sooviksin vastata küsimusele, mis ilmselt on tekkinud meie kallitel lugejatel. Miks on vajalik selline „vaev” languste osas, miks me peame ronima mäge üles ja alla (peamiselt alla), et leida ihaldatud madalikku, kui meie käes on nii võimas ja lihtne seade, nagu analüütiline lahendus, mis transportib meid koheselt vajalikku kohta?

Sellele küsimusele on vastus pinnapealne. Oleme käsitlenud väga lihtsat näidet, kus tõeline vastus sõltub ühest tunnusest . Elus ei juhtu selliseid olukordi tihti, seega kujutame ette, et meil on 2, 30, 50 või rohkem tunnust. Lisame sellele tuhandeid, kui mitte kümneid tuhandeid väärtusi iga tunnuse jaoks. Sellisel juhul ei pruugi analüütiline lahendus eksamile vastu pidada ja võib läbikukkuda. Samal ajal läheneb gradientlangus ja selle variatsioonid aeglaselt, kuid kindlalt meie eesmärgile — funktsiooni miinimumile. Ja kiirusest ei tasu muretseda — me võtame kindlasti veelkord läbi meetodeid, mis võimaldavad meil seada ja reguleerida sammu pikkust (st kiirus).

Nüüd on lühike kokkuvõte.

Esiteks, loodan, et artiklis esitatud materjal aitab algaval "andmete teadlastel" mõista, kuidas lahendada lihtsa (ja mitte ainult) lineaarse regressiooni võrrandeid.

Teiseks, me vaatasime üle mitmeid viise võrrandi lahendamiseks. Nüüd, sõltuvalt olukorrast, saame valida selle, mis sobib kõige paremini püstitatud ülesande lahendamiseks.

Kolmandaks, nägime täiendavate seadistuste jõudu, nimelt gradientide pikkust. Seda parameetrit ei tohi alahinnata. Nagu eespool märgitud, peaksite kulude vähendamiseks muutma sammupikkust langemise käigus.

Neljandaks, meie puhul näitasid "oma kirjutatud" funktsioonid parima ajakulu arvutustes. Tõenäoliselt on see seotud mitte kõige professionaalsema raamatukogu võimaluste rakendamisega. NumPy. Kuid igal juhul järeldus on järgmine. Ühelt poolt väärib mõnikord kahtlustama väljakujunenud arvamusi ja teiselt poolt ei ole alati mõistlik kõike keeruliseks muuta - vastupidi, mõnikord osutub efektiivsemaks lihtsam lahendus. Ja kuna meie eesmärk oli uurida kolme lähenemist lihtsa lineaarse regressiooni võrrandi lahendamiseks, siis piisab meile "oma kirjutatud" funktsioonide kasutamisest.

Literatuur (või midagi sarnast)

1. Lineaarne regressioon

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

2. Vähemate ruutude meetod

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Tuletis

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Gradientide tõus

habr.com/ru/post/471458

habr.com/ru/post/307312

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

6. NumPy teek

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

Allikas: habr.com