Článok rozoberá niekoľko spôsobov, ako určiť matematickú rovnicu jednoduchej (párovej) regresnej priamky.
Všetky metódy riešenia tu diskutovanej rovnice sú založené na metóde najmenších štvorcov. Označme metódy takto:
- Analytický roztok
- Gradientný zostup
- Stochastický gradientový zostup
Pre každú metódu riešenia rovnice priamky poskytuje článok rôzne funkcie, ktoré sú rozdelené najmä na tie, ktoré sú napísané bez použitia knižnice nemotorný a tie, ktoré sa používajú na výpočty nemotorný. Verí sa, že zručné použitie nemotorný zníži náklady na výpočtovú techniku.
Všetok kód uvedený v článku je napísaný v jazyku python 2.7 s Jupyter Notebook. Zdrojový kód a súbor so vzorovými údajmi sú zverejnené na
Článok je viac zameraný ako na začiatočníkov, tak aj na tých, ktorí už postupne začali zvládať štúdium veľmi širokej sekcie umelej inteligencie – strojového učenia.
Na ilustráciu materiálu použijeme veľmi jednoduchý príklad.
Príklady podmienok
Máme päť hodnôt, ktoré charakterizujú závislosť Y od X (Tabuľka č. 1):
Tabuľka č. 1 „Príkladové podmienky“
Budeme predpokladať, že hodnoty je mesiac v roku a — príjmy za tento mesiac. Inými slovami, príjmy závisia od mesiaca v roku a - jediné znamenie, od ktorého závisia príjmy.
Príklad je taký-tak, a to z hľadiska podmienenej závislosti výnosov od mesiaca v roku, ako aj z hľadiska počtu hodnôt - je ich veľmi málo. Takéto zjednodušenie však umožní, ako sa hovorí, vysvetliť, nie vždy ľahko, materiál, ktorý si začiatočníci osvojujú. A tiež jednoduchosť čísel umožní tým, ktorí chcú vyriešiť príklad na papieri bez značných nákladov na prácu.
Predpokladajme, že závislosť uvedenú v príklade možno celkom dobre aproximovať matematickou rovnicou jednoduchej (párovej) regresnej priamky v tvare:
kde je mesiac, v ktorom bola tržba prijatá, — príjem zodpovedajúci mesiacu, и sú regresné koeficienty odhadovanej čiary.
Všimnite si, že koeficient často nazývaný sklon alebo sklon odhadovanej čiary; predstavuje sumu, o ktorú keď sa zmení .
Je zrejmé, že našou úlohou v príklade je vybrať takéto koeficienty v rovnici и , pri ktorých sú odchýlky nami vypočítaných hodnôt tržieb za mesiac od skutočných odpovedí, t.j. hodnoty uvedené vo vzorke budú minimálne.
Metóda najmenších štvorcov
Podľa metódy najmenších štvorcov by sa odchýlka mala vypočítať jej umocnením. Táto technika vám umožňuje vyhnúť sa vzájomnému zrušeniu odchýlok, ak majú opačné znamienka. Napríklad, ak je v jednom prípade odchýlka +5 (plus päť) a v druhom -5 (mínus päť), potom sa súčet odchýlok navzájom vyruší a bude mať hodnotu 0 (nula). Je možné odchýlku neodmocniť, ale využiť vlastnosť modulu a potom budú všetky odchýlky kladné a budú sa kumulovať. Nebudeme sa týmto bodom podrobne zaoberať, ale jednoducho naznačíme, že pre pohodlie výpočtov je zvyčajné umocniť odchýlku na druhú.
Takto vyzerá vzorec, pomocou ktorého určíme najmenší súčet štvorcových odchýlok (chýb):
kde je funkciou aproximácie pravdivých odpovedí (t. j. príjmov, ktoré sme vypočítali),
sú pravdivé odpovede (príjmy uvedené vo vzorke),
je výberový index (číslo mesiaca, v ktorom sa odchýlka určuje)
Diferencujme funkciu, definujme parciálne diferenciálne rovnice a buďme pripravení prejsť na analytické riešenie. Najprv si však urobme krátky exkurz o tom, čo je to diferenciácia a pripomeňme si geometrický význam derivácie.
Diferenciácia
Diferenciácia je operácia hľadania derivácie funkcie.
Na čo sa derivát používa? Derivácia funkcie charakterizuje rýchlosť zmeny funkcie a udáva jej smer. Ak je derivácia v danom bode kladná, funkcia sa zvyšuje, v opačnom prípade funkcia klesá. A čím väčšia je hodnota absolútnej derivácie, tým vyššia je rýchlosť zmeny funkčných hodnôt, ako aj strmší sklon grafu funkcie.
Napríklad v podmienkach karteziánskeho súradnicového systému sa hodnota derivácie v bode M(0,0) rovná + 25 znamená, že v danom bode, keď je hodnota posunutá vpravo o konvenčnú jednotku, hodnotu zvyšuje o 25 konvenčných jednotiek. Na grafe to vyzerá na pomerne strmý nárast hodnôt z daného bodu.
Ďalší príklad. Hodnota derivátu je rovnaká -0,1 znamená, že pri premiestnení na jednu konvenčnú jednotku, hodnotu klesá len o 0,1 konvenčnej jednotky. Zároveň na grafe funkcie môžeme pozorovať sotva badateľný sklon smerom nadol. Ak nakreslíme analógiu s horou, je to, ako keby sme veľmi pomaly schádzali miernym svahom z hory, na rozdiel od predchádzajúceho príkladu, kde sme museli vyliezť na veľmi strmé štíty :)
Teda po diferenciácii funkcie podľa šance и , definujeme parciálne diferenciálne rovnice 1. rádu. Po určení rovníc dostaneme sústavu dvoch rovníc, ktorých riešením budeme môcť vybrať také hodnoty koeficientov и , pre ktoré sa hodnoty zodpovedajúcich derivácií v daných bodoch menia o veľmi, veľmi malé množstvo a v prípade analytického riešenia sa nemenia vôbec. Inými slovami, chybová funkcia pri nájdených koeficientoch dosiahne minimum, pretože hodnoty parciálnych derivácií v týchto bodoch sa budú rovnať nule.
Takže podľa pravidiel diferenciácie parciálna derivačná rovnica 1. rádu vzhľadom na koeficient bude mať podobu:
Parciálna derivačná rovnica 1. rádu vzhľadom na bude mať podobu:
V dôsledku toho sme dostali systém rovníc, ktorý má pomerne jednoduché analytické riešenie:
začať{rovnica*}
začať{prípady}
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}
koniec{rovnica*}
Pred riešením rovnice predbežne načítajme, skontrolujeme správnosť načítania a naformátujeme údaje.
Načítavanie a formátovanie údajov
Treba poznamenať, že vzhľadom na skutočnosť, že pre analytické riešenie a následne pre gradientový a stochastický gradientový zostup, použijeme kód v dvoch variantoch: pomocou knižnice nemotorný a bez toho, aby sme ho použili, potom budeme potrebovať vhodné formátovanie údajov (pozri kód).
Kód načítania a spracovania údajov
# импортируем все нужные нам библиотеки
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 '********************************************'Vizualizácia
Teraz, keď sme po prvé načítali údaje, po druhé skontrolovali správnosť načítania a nakoniec údaje naformátovali, vykonáme prvú vizualizáciu. Často sa na to používa metóda párový pozemok knižnica morský. V našom príklade kvôli obmedzenému počtu nemá zmysel používať knižnicu morský. Budeme používať bežnú knižnicu matplotlib a stačí sa pozrieť na bodový graf.
Kód bodového grafu
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()Graf č. 1 „Závislosť príjmov od mesiaca v roku“
Analytický roztok
Využime najbežnejšie nástroje v krajta a vyriešiť sústavu rovníc:
začať{rovnica*}
začať{prípady}
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}
koniec{rovnica*}
Podľa Cramerovho pravidla nájdeme všeobecný determinant, ako aj determinanty podľa a , po ktorom sa determinant delí o k všeobecnému determinantu - nájdite koeficient , podobne nájdeme koeficient .
Kód analytického roztoku
# определим функцию для расчета коэффициентов 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)Tu je to, čo máme:
Takže boli nájdené hodnoty koeficientov, bol stanovený súčet štvorcových odchýlok. Narysujme na histograme rozptylu priamku podľa zistených koeficientov.
Kód regresnej čiary
# определим функцию для формирования массива рассчетных значений выручки
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()Tabuľka č. 2 „Správne a vypočítané odpovede“
Môžete sa pozrieť na graf odchýlok pre každý mesiac. V našom prípade z toho nevyvodíme žiadnu významnú praktickú hodnotu, ale uspokojíme našu zvedavosť, ako dobre jednoduchá lineárna regresná rovnica charakterizuje závislosť tržieb od mesiaca v roku.
Kód grafu odchýlok
# определим функцию для формирования массива отклонений в процентах
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()Graf č. 3 „Odchýlky, %“
Nie je to dokonalé, ale svoju úlohu sme splnili.
Napíšme funkciu, ktorá určí koeficienty и používa knižnicu nemotorný, presnejšie napíšeme dve funkcie: jednu pomocou pseudoinverznej matice (v praxi sa neodporúča, keďže proces je výpočtovo zložitý a nestabilný), druhú pomocou maticovej rovnice.
Kód analytického riešenia (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_npPorovnajme čas strávený určovaním koeficientov и , v súlade s 3 prezentovanými metódami.
Kód pre výpočet času výpočtu
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)
S malým množstvom údajov sa pred nami objaví „vlastne napísaná“ funkcia, ktorá nájde koeficienty pomocou Cramerovej metódy.
Teraz môžete prejsť na iné spôsoby hľadania koeficientov и .
Gradientný zostup
Najprv si definujme, čo je gradient. Jednoducho povedané, gradient je segment, ktorý udáva smer maximálneho rastu funkcie. Analogicky s výstupom na horu, kde stúpanie je tam, kde je najstrmšie stúpanie na vrchol hory. Pri rozvíjaní príkladu s horou si pamätáme, že v skutočnosti potrebujeme najstrmší zostup, aby sme čo najrýchlejšie dosiahli nížinu, to znamená minimum - miesto, kde sa funkcia nezvýši ani nezníži. V tomto bode bude derivácia rovná nule. Preto nepotrebujeme gradient, ale antigradient. Ak chcete nájsť antigradient, stačí vynásobiť gradient -1 (mínus jedna).
Venujme pozornosť skutočnosti, že funkcia môže mať niekoľko miním a po zostúpení do jedného z nich pomocou nižšie navrhnutého algoritmu nenájdeme ďalšie minimum, ktoré môže byť nižšie ako nájdené. Uvoľnime sa, toto nám nehrozí! V našom prípade máme do činenia s jediným minimom, keďže naša funkcia na grafe je pravidelná parabola. A ako všetci dobre vieme z nášho školského kurzu matematiky, parabola má len jedno minimum.
Potom, čo sme zistili, prečo potrebujeme gradient, a tiež, že gradient je segment, teda vektor s danými súradnicami, čo sú presne tie isté koeficienty и môžeme realizovať gradientný zostup.
Pred začatím odporúčam prečítať si len pár viet o zostupovom algoritme:
- Súradnice koeficientov určujeme pseudonáhodným spôsobom и . V našom príklade určíme koeficienty blízke nule. Toto je bežná prax, ale každý prípad môže mať svoju vlastnú prax.
- Zo súradnice odpočítajte hodnotu parciálnej derivácie 1. rádu v bode . Ak je teda derivácia kladná, funkcia sa zvyšuje. Odčítaním hodnoty derivátu sa teda posunieme v opačnom smere rastu, teda v smere zostupu. Ak je derivácia záporná, funkcia v tomto bode klesá a odčítaním hodnoty derivácie sa pohybujeme v smere klesania.
- Podobnú operáciu vykonávame so súradnicou : odpočítajte hodnotu parciálnej derivácie v bode .
- Aby ste nepreskočili minimum a neleteli do hlbokého vesmíru, je potrebné nastaviť veľkosť kroku v smere zostupu. Vo všeobecnosti by sa dal napísať celý článok o tom, ako správne nastaviť krok a ako ho zmeniť počas procesu zostupu, aby sa znížili výpočtové náklady. Teraz však máme pred sebou trochu inú úlohu a veľkosť kroku určíme pomocou vedeckej metódy „poke“ alebo, ako sa hovorí v bežnej reči, empiricky.
- Keď už sme z daných súradníc и odčítaním hodnôt derivátov dostaneme nové súradnice и . Urobíme ďalší krok (odčítanie), už z vypočítaných súradníc. A tak sa cyklus začína znova a znova, až kým sa nedosiahne požadovaná konvergencia.
Všetky! Teraz sme pripravení ísť hľadať najhlbšiu roklinu priekopy Mariana. Začnime.
Kód pre gradientový zostup
# напишем функцию градиентного спуска без использования библиотеки 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
Ponorili sme sa až na samé dno Mariánskej priekopy a tam sme našli všetky rovnaké hodnoty koeficientov и , čo sa presne dalo očakávať.
Poďme ešte raz, len tentoraz bude naše hlbokomorské vozidlo naplnené inými technológiami, a to knižnicou nemotorný.
Kód pre gradientový zostup (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
Hodnoty koeficientov и nezmeniteľné.
Pozrime sa, ako sa menila chyba počas klesania gradientu, teda ako sa menil súčet štvorcových odchýlok s každým krokom.
Kód na vykreslenie súčtu kvadrátov odchýlok
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()Graf č. 4 „Súčet štvorcových odchýlok počas klesania gradientom“
Na grafe vidíme, že každým krokom sa chyba zmenšuje a po určitom počte iterácií pozorujeme takmer vodorovnú čiaru.
Nakoniec odhadnime rozdiel v čase vykonania kódu:
Kód na určenie času výpočtu zostupu gradientu
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)
Možno robíme niečo zle, ale opäť ide o jednoduchú „domácu“ funkciu, ktorá nepoužíva knižnicu nemotorný prekonáva čas výpočtu funkcie pomocou knižnice nemotorný.
Ale nestojíme na mieste, ale smerujeme k štúdiu ďalšieho vzrušujúceho spôsobu riešenia jednoduchej lineárnej regresnej rovnice. Zoznámte sa!
Stochastický gradientový zostup
Aby sme rýchlo pochopili princíp fungovania stochastického gradientového zostupu, je lepšie určiť jeho rozdiely od bežného gradientového zostupu. My, v prípade gradientu klesania, v rovniciach derivácií o и použil súčty hodnôt všetkých funkcií a pravdivých odpovedí dostupných vo vzorke (to znamená súčty všetkých и ). Pri stochastickom zostupe gradientu nepoužijeme všetky hodnoty prítomné vo vzorke, ale namiesto toho pseudonáhodne vyberieme takzvaný index vzorky a použijeme jeho hodnoty.
Napríklad, ak je index určený ako číslo 3 (tri), potom vezmeme hodnoty и , potom dosadíme hodnoty do derivačných rovníc a určíme nové súradnice. Potom, po určení súradníc, opäť pseudonáhodne určíme index vzorky, nahradíme hodnoty zodpovedajúce indexu do parciálnych diferenciálnych rovníc a určíme súradnice novým spôsobom. и atď. kým sa konvergencia nezmení na zelenú. Na prvý pohľad sa môže zdať, že to vôbec nemôže fungovať, ale je to tak. Je pravda, že stojí za zmienku, že chyba sa neznižuje každým krokom, ale určite existuje tendencia.
Aké sú výhody stochastického gradientového zostupu oproti konvenčnému? Ak je veľkosť našej vzorky veľmi veľká a meria sa v desiatkach tisíc hodnôt, potom je oveľa jednoduchšie spracovať, povedzme, náhodných tisíc z nich, a nie celú vzorku. Tu vstupuje do hry stochastický gradientový zostup. V našom prípade si, samozrejme, nevšimneme veľký rozdiel.
Pozrime sa na kód.
Kód pre stochastický gradientový zostup
# определим функцию стох.град.шага
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])
Pozorne sa pozrieme na koeficienty a pristihneme sa, ako si kladieme otázku "Ako je to možné?" Dostali sme iné hodnoty koeficientov и . Možno stochastický gradientový zostup našiel optimálnejšie parametre pre rovnicu? Bohužiaľ nie. Stačí sa pozrieť na súčet štvorcových odchýlok a vidieť, že s novými hodnotami koeficientov je chyba väčšia. Nikam sa neponáhľame do zúfalstva. Zostavme graf zmeny chyby.
Kód na vykreslenie súčtu kvadrátov odchýlok v zostupe stochastického gradientu
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()Graf č. 5 „Súčet štvorcových odchýlok počas stochastického gradientového zostupu“
Pri pohľade na harmonogram všetko zapadá a teraz všetko napravíme.
Takže, čo sa stalo? Stalo sa nasledovné. Keď náhodne vyberieme mesiac, potom sa náš algoritmus snaží znížiť chyby vo výpočte výnosov pre vybraný mesiac. Potom vyberieme ďalší mesiac a výpočet zopakujeme, ale znížime chybu pre druhý vybraný mesiac. Teraz si pamätajte, že prvé dva mesiace sa výrazne odchyľujú od línie jednoduchej lineárnej regresnej rovnice. To znamená, že keď sa vyberie ktorýkoľvek z týchto dvoch mesiacov, znížením chyby každého z nich náš algoritmus vážne zvýši chybu pre celú vzorku. Čo teda robiť? Odpoveď je jednoduchá: musíte znížiť krok zostupu. Koniec koncov, znížením kroku zostupu, chyba tiež prestane „skákať“ hore a dole. Alebo skôr, chyba „skákania“ sa nezastaví, ale neurobí to tak rýchlo :) Poďme skontrolovať.
Kód na spustenie SGD s menšími prírastkami
# запустим функцию, уменьшив шаг в 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()
Graf č. 6 „Súčet štvorcových odchýlok pri stochastickom gradientovom zostupe (80 tisíc krokov)“
Koeficienty sa zlepšili, ale stále nie sú ideálne. Hypoteticky sa to dá takto napraviť. Vyberáme napríklad v posledných 1000 iteráciách hodnoty koeficientov, s ktorými bola urobená minimálna chyba. Je pravda, že na to budeme musieť zapísať aj hodnoty samotných koeficientov. Nebudeme to robiť, ale radšej dbajte na harmonogram. Vyzerá to hladko a zdá sa, že chyba sa znižuje rovnomerne. V skutočnosti to nie je pravda. Pozrime sa na prvých 1000 iterácií a porovnajme ich s poslednými.
Kód pre graf SGD (prvých 1000 krokov)
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()Graf č. 7 „Súčet štvorcových odchýlok SGD (prvých 1000 krokov)“
Graf č. 8 „Súčet štvorcových odchýlok SGD (posledných 1000 krokov)“
Hneď na začiatku klesania pozorujeme pomerne rovnomerný a strmý pokles chyby. V posledných iteráciách vidíme, že chyba sa pohybuje okolo a okolo hodnoty 1,475 a v niektorých momentoch sa dokonca rovná tejto optimálnej hodnote, ale potom stále stúpa... Opakujem, môžete si zapísať hodnoty koeficienty и a potom vyberte tie, pri ktorých je chyba minimálna. Mali sme však vážnejší problém: museli sme urobiť 80 XNUMX krokov (pozri kód), aby sme dosiahli hodnoty blízke optimálnym. A to je už v rozpore s myšlienkou šetrenia výpočtového času pomocou stochastického klesania gradientu v porovnaní s klesaním gradientu. Čo sa dá opraviť a zlepšiť? Nie je ťažké si všimnúť, že v prvých iteráciách s istotou klesáme nadol, a preto by sme v prvých iteráciách mali nechať veľký krok a pri postupe vpred krok znižovať. V tomto článku to neurobíme - je už príliš dlhý. Kto chce, môže si sám premyslieť, ako na to, nie je to ťažké :)
Teraz urobme stochastický gradientový zostup pomocou knižnice nemotorný (a nezakopnime o kamene, ktoré sme predtým identifikovali)
Kód pre zostup stochastického gradientu (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
Hodnoty sa ukázali byť takmer rovnaké ako pri zostupe bez použitia nemotorný. Je to však logické.
Poďme zistiť, ako dlho nám trvali stochastické gradientové zjazdy.
Kód na určenie času výpočtu SGD (80 tisíc krokov)
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)
Čím ďalej do lesa, tým sú oblaky tmavšie: najlepší výsledok opäť ukazuje „samomopísaný“ vzorec. To všetko naznačuje, že musia existovať ešte jemnejšie spôsoby využitia knižnice nemotorný, ktoré skutočne urýchľujú výpočtové operácie. V tomto článku sa o nich nedozvieme. Vo voľnom čase bude o čom premýšľať :)
Zhrnieme
Pred zhrnutím by som rád odpovedal na otázku, ktorá s najväčšou pravdepodobnosťou vyvstala od nášho drahého čitateľa. Prečo vlastne také „mučenie“ zostupmi, prečo musíme chodiť hore a dolu (väčšinou dole), aby sme našli vzácnu nížinu, ak máme v rukách taký výkonný a jednoduchý prístroj, v formou analytického riešenia, ktoré nás okamžite teleportuje na správne miesto?
Odpoveď na túto otázku leží na povrchu. Teraz sme sa pozreli na veľmi jednoduchý príklad, v ktorom je pravdivá odpoveď závisí od jedného znamenia . Toto sa v živote často nevidí, takže si predstavme, že máme 2, 30, 50 alebo viac znamení. Pridajme k tomu tisíce alebo dokonca desaťtisíce hodnôt pre každý atribút. V tomto prípade analytický roztok nemusí vydržať test a zlyhať. Gradientový zostup a jeho variácie nás zase pomaly, ale isto približujú k cieľu – minimu funkcie. A nebojte sa rýchlosti – pravdepodobne sa pozrieme na spôsoby, ktoré nám umožnia nastaviť a regulovať dĺžku kroku (teda rýchlosť).
A teraz skutočné krátke zhrnutie.
Po prvé, dúfam, že materiál prezentovaný v článku pomôže začínajúcim „dátovým vedcom“ pochopiť, ako riešiť jednoduché (nielen) lineárne regresné rovnice.
Po druhé, pozreli sme sa na niekoľko spôsobov riešenia rovnice. Teraz, v závislosti od situácie, môžeme vybrať ten, ktorý je najvhodnejší na vyriešenie problému.
Po tretie, videli sme silu ďalších nastavení, konkrétne dĺžky kroku zostupu. Tento parameter nemožno zanedbať. Ako je uvedené vyššie, aby sa znížili náklady na výpočty, dĺžka kroku by sa mala počas zostupu meniť.
Po štvrté, v našom prípade „domáce písané“ funkcie vykazovali najlepšie časové výsledky pre výpočty. Je to zrejme spôsobené nie práve najprofesionálnejším využitím možností knižnice nemotorný. Ale nech je to akokoľvek, nasledujúci záver naznačuje sám seba. Na jednej strane sa niekedy oplatí spochybniť ustálené názory a na druhej strane nie vždy sa oplatí všetko komplikovať – naopak, niekedy je jednoduchší spôsob riešenia problému efektívnejší. A keďže naším cieľom bolo analyzovať tri prístupy k riešeniu jednoduchej lineárnej regresnej rovnice, úplne nám stačilo použitie „vlastne napísaných“ funkcií.
literatúra (alebo niečo podobné)
1. Lineárna regresia
2. Metóda najmenších štvorcov
3. Derivát
4. Gradient
5. Gradientný zostup
6. Knižnica NumPy
Zdroj: hab.com