Článek pojednává o několika způsobech, jak určit matematickou rovnici jednoduché (párové) regresní přímky.
Všechny metody řešení rovnice zde diskutované jsou založeny na metodě nejmenších čtverců. Označme metody takto:
- Analytické řešení
- Gradientní sestup
- Stochastický gradient klesání
Pro každou metodu řešení rovnice přímky poskytuje článek různé funkce, které se dělí především na ty, které jsou psány bez použití knihovny nemotorný a ty, které se používají pro výpočty nemotorný. To je věřil, že dovedné použití nemotorný sníží náklady na výpočetní techniku.
Veškerý kód uvedený v článku je napsán v jazyce python 2.7 použití Jupyter Notebook. Zdrojový kód a soubor s ukázkovými daty jsou zveřejněny na
Článek je zaměřen spíše jak na začátečníky, tak na ty, kteří již postupně začali zvládat studium velmi široké sekce umělé inteligence – strojového učení.
Pro ilustraci materiálu použijeme velmi jednoduchý příklad.
Příklad podmínek
Máme pět hodnot, které charakterizují závislost Y z X (Tabulka č. 1):
Tabulka č. 1 „Příkladové podmínky“
Budeme předpokládat, že hodnoty je měsíc v roce a — tržby za tento měsíc. Jinými slovy, tržby závisí na měsíci v roce a - jediné znamení, na kterém závisí výnos.
Příklad je takový, a to jak z hlediska podmíněné závislosti výnosů na měsíci v roce, tak z hlediska počtu hodnot - je jich velmi málo. Takové zjednodušení však umožní, jak se říká, vysvětlit, ne vždy snadno, látku, kterou začátečníci asimilují. A také jednoduchost čísel umožní těm, kteří chtějí vyřešit příklad na papíře, bez významných nákladů na pracovní sílu.
Předpokládejme, že závislost uvedenou v příkladu lze docela dobře aproximovat matematickou rovnicí jednoduché (párové) regresní přímky tvaru:
kde je měsíc, ve kterém byla tržba přijata, — příjmy odpovídající měsíci, и jsou regresní koeficienty odhadované přímky.
Všimněte si, že koeficient často nazývaný sklon nebo sklon odhadované čáry; představuje částku, o kterou když se to změní .
Je zřejmé, že naším úkolem v příkladu je vybrat takové koeficienty v rovnici и , při které jsou odchylky námi vypočtených hodnot tržeb po měsících od skutečných odpovědí, tzn. hodnoty uvedené ve vzorku budou minimální.
Metoda nejmenších čtverců
Podle metody nejmenších čtverců by se odchylka měla vypočítat jejím umocněním. Tato technika vám umožňuje vyhnout se vzájemnému zrušení odchylek, pokud mají opačná znaménka. Pokud je například v jednom případě odchylka +5 (plus pět) a ve druhém -5 (mínus pět), pak se součet odchylek vzájemně vyruší a bude činit 0 (nula). Je možné odchylku neodmocnit, ale využít vlastnosti modulu a pak budou všechny odchylky kladné a budou se kumulovat. Nebudeme se tímto bodem podrobně zabývat, ale jednoduše naznačíme, že pro usnadnění výpočtů je obvyklé odchylku odmocnit.
Takto vypadá vzorec, pomocí kterého určíme nejmenší součet čtverců odchylek (chyb):
kde je funkcí aproximace pravdivých odpovědí (tj. příjmů, které jsme vypočítali),
jsou pravdivé odpovědi (příjmy uvedené ve vzorku),
je výběrový index (číslo měsíce, ve kterém se odchylka určuje)
Derivujme funkci, definujme parciální diferenciální rovnice a buďme připraveni přejít k analytickému řešení. Nejprve si ale udělejme krátký exkurz o tom, co je to diferenciace, a připomeňme si geometrický význam derivace.
Diferenciace
Diferenciace je operace hledání derivace funkce.
K čemu se derivát používá? Derivace funkce charakterizuje rychlost změny funkce a říká nám její směr. Pokud je derivace v daném bodě kladná, pak funkce roste, v opačném případě funkce klesá. A čím větší je hodnota absolutní derivace, tím vyšší je rychlost změny funkčních hodnot a také strmější sklon grafu funkce.
Například za podmínek kartézského souřadnicového systému je hodnota derivace v bodě M(0,0) rovna +25 znamená, že v daném bodě, když je hodnota posunuta vpravo o konvenční jednotku, hodnotu zvyšuje o 25 konvenčních jednotek. Na grafu to vypadá na poměrně strmý nárůst hodnot z daného bodu.
Další příklad. Hodnota derivace se rovná -0,1 znamená, že při přemístění na jednu konvenční jednotku, hodnotu klesá pouze o 0,1 konvenční jednotky. Na grafu funkce přitom můžeme pozorovat sotva znatelný sklon dolů. Když nakreslíme analogii s horou, je to, jako bychom velmi pomalu sestupovali z mírného svahu z hory, na rozdíl od předchozího příkladu, kde jsme museli šplhat na velmi strmé vrcholy :)
Tedy po diferenciaci funkce podle šance и , definujeme parciální diferenciální rovnice 1. řádu. Po určení rovnic obdržíme soustavu dvou rovnic, jejichž řešením budeme moci vybrat takové hodnoty koeficientů и , u kterého se hodnoty odpovídajících derivací v daných bodech mění velmi, velmi málo a v případě analytického řešení se nemění vůbec. Jinými slovy, chybová funkce u nalezených koeficientů dosáhne minima, protože hodnoty parciálních derivací v těchto bodech se budou rovnat nule.
Takže podle pravidel derivace parciální derivační rovnice 1. řádu vzhledem ke koeficientu bude mít podobu:
Parciální derivační rovnice 1. řádu vzhledem k bude mít podobu:
V důsledku toho jsme dostali systém rovnic, který má poměrně jednoduché analytické řešení:
začít{rovnice*}
začít{pří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
konec{případy}
konec{rovnice*}
Před řešením rovnice předběžně načteme, zkontrolujeme správnost načítání a naformátujeme data.
Načítání a formátování dat
Je třeba poznamenat, že vzhledem k tomu, že pro analytické řešení a následně pro gradientový a stochastický sestup gradientu, použijeme kód ve dvou variantách: pomocí knihovny nemotorný a bez použití, pak budeme potřebovat vhodné formátování dat (viz kód).
Kód načítání a zpracování dat
# импортируем все нужные нам библиотеки
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 '********************************************'Vizualizace
Nyní, když jsme za prvé načetli data, za druhé zkontrolovali správnost načítání a nakonec data naformátovali, provedeme první vizualizaci. K tomu se často používá metoda párový graf knihovna mořský. V našem příkladu kvůli omezenému počtu nemá smysl používat knihovnu mořský. Budeme používat běžnou knihovnu matplotlib a stačí se podívat na bodový graf.
Kód rozptylové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ávislost příjmů na měsíci v roce“
Analytické řešení
Použijme nejběžnější nástroje v krajta a vyřešit soustavu rovnic:
začít{rovnice*}
začít{pří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
konec{případy}
konec{rovnice*}
Podle Cramerova pravidla najdeme obecný determinant, stejně jako determinanty podle a , načež vydělením determinantu k obecnému determinantu - najděte koeficient , podobně zjistíme 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)Zde je to, co máme:
Takže byly nalezeny hodnoty koeficientů, byl stanoven součet čtverců odchylek. Nakreslime na histogramu rozptylu přímku v souladu s nalezenými koeficienty.
Kód regresní čáry
# определим функцию для формирования массива рассчетных значений выручки
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()Tabulka č. 2 „Správné a vypočítané odpovědi“
Můžete se podívat na graf odchylek pro každý měsíc. V našem případě z toho nevyvodíme žádnou významnou praktickou hodnotu, ale uspokojíme svou zvědavost, jak dobře jednoduchá lineární regresní rovnice charakterizuje závislost tržeb na měsíci v roce.
Kód diagramu odchylek
# определим функцию для формирования массива отклонений в процентах
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 „Odchylky, %“
Ne dokonalé, ale svůj úkol jsme splnili.
Napišme funkci, která určí koeficienty и používá knihovnu nemotorný, přesněji řečeno, napíšeme dvě funkce: jednu pomocí pseudoinverzní matice (v praxi se nedoporučuje, protože proces je výpočetně složitý a nestabilní), druhou pomocí maticové rovnice.
Kód analytického řešení (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_npPorovnejme čas strávený stanovením koeficientů и , v souladu se 3 prezentovanými metodami.
Kód pro výpočet doby 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žstvím dat se dopředu objeví „vlastně napsaná“ funkce, která najde koeficienty pomocí Cramerovy metody.
Nyní můžete přejít k dalším způsobům hledání koeficientů и .
Gradientní sestup
Nejprve si definujme, co je gradient. Jednoduše řečeno, gradient je segment, který udává směr maximálního růstu funkce. Analogicky k lezení na horu, kde stoupání je tam, kde je nejstrmější stoupání na vrchol hory. Při rozvíjení příkladu s horou si pamatujeme, že ve skutečnosti potřebujeme nejprudší sestup, abychom co nejrychleji dosáhli nížiny, tedy minima - místa, kde se funkce nezvyšuje ani nesnižuje. V tomto bodě bude derivace rovna nule. Proto nepotřebujeme gradient, ale antigradient. K nalezení antigradientu stačí gradient vynásobit -1 (mínus jedna).
Věnujme pozornost tomu, že funkce může mít několik minim a po sestupu do jednoho z nich pomocí níže navrženého algoritmu nenajdeme další minimum, které může být nižší než to nalezené. Uvolníme se, to nám nehrozí! V našem případě se jedná o jediné minimum, protože naše funkce na grafu je pravidelná parabola. A jak všichni dobře víme z našeho školního kurzu matematiky, parabola má jen jedno minimum.
Poté, co jsme zjistili, proč potřebujeme gradient, a také, že gradient je segment, tedy vektor s danými souřadnicemi, což jsou přesně stejné koeficienty и můžeme implementovat gradientní sestup.
Než začnete, doporučuji přečíst si jen pár vět o sestupovém algoritmu:
- Pseudonáhodným způsobem určujeme souřadnice koeficientů и . V našem příkladu budeme definovat koeficienty blízké nule. Toto je běžná praxe, ale každý případ může mít svou vlastní praxi.
- Z koordinátu odečtěte hodnotu parciální derivace 1. řádu v bodě . Pokud je tedy derivace kladná, funkce se zvyšuje. Odečtením hodnoty derivace se tedy budeme pohybovat v opačném směru růstu, tedy ve směru sestupu. Pokud je derivace záporná, pak funkce v tomto bodě klesá a odečtením hodnoty derivace se pohybujeme ve směru sestupu.
- Podobnou operaci provádíme se souřadnicí : odečtěte hodnotu parciální derivace v bodě .
- Abyste nepřeskočili minimum a neletěli do hlubokého vesmíru, je nutné nastavit velikost kroku ve směru sestupu. Obecně by se dal napsat celý článek o tom, jak správně nastavit krok a jak jej změnit během procesu sestupu, aby se snížily výpočetní náklady. Nyní nás ale čeká trochu jiný úkol a velikost kroku určíme vědeckou metodou „šťouchání“ nebo, jak se v běžné řeči říká, empiricky.
- Jakmile jsme z uvedených souřadnic и odečteme hodnoty derivací, získáme nové souřadnice и . Provedeme další krok (odčítání), již z vypočtených souřadnic. A tak cyklus začíná znovu a znovu, dokud není dosaženo požadované konvergence.
Všechno! Nyní jsme připraveni vydat se hledat nejhlubší soutěsku Marianského příkopu. Začněme.
Kód pro gradientní klesání
# напишем функцию градиентного спуска без использования библиотеки 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
Ponořili jsme se až na samé dno Mariánského příkopu a tam jsme našli všechny stejné hodnoty koeficientů и , což je přesně to, co se dalo čekat.
Pojďme se ještě jednou ponořit, ale tentokrát bude naše hlubokomořské plavidlo naplněno dalšími technologiemi, konkrétně knihovnou nemotorný.
Kód pro gradientní klesání (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 koeficientů и neměnný.
Podívejme se, jak se měnila chyba při klesání gradientu, tedy jak se s každým krokem měnil součet kvadrátů odchylek.
Kód pro vykreslení součtů kvadrátů odchylek
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 „Součet kvadrátů odchylek při klesání gradientem“
Na grafu vidíme, že s každým krokem se chyba zmenšuje a po určitém počtu iterací pozorujeme téměř vodorovnou čáru.
Nakonec odhadneme rozdíl v době provádění kódu:
Kód pro určení doby výpočtu klesání 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žná děláme něco špatně, ale opět je to jednoduchá „doma psaná“ funkce, která nepoužívá knihovnu nemotorný překonává dobu výpočtu funkce pomocí knihovny nemotorný.
Ale nestojíme na místě, ale směřujeme ke studiu dalšího vzrušujícího způsobu řešení jednoduché lineární regresní rovnice. Setkat!
Stochastický gradient klesání
Abychom rychle pochopili princip fungování stochastického gradientního sestupu, je lepší určit jeho rozdíly od běžného gradientního sestupu. My, v případě gradientu klesání, v rovnicích derivací z и použil součty hodnot všech funkcí a pravdivých odpovědí dostupných ve vzorku (tj. součty všech и ). Při sestupu stochastického gradientu nepoužijeme všechny hodnoty přítomné ve vzorku, ale místo toho pseudonáhodně vybereme takzvaný index vzorku a použijeme jeho hodnoty.
Pokud je například index určen jako číslo 3 (tři), pak vezmeme hodnoty и , pak dosadíme hodnoty do derivačních rovnic a určíme nové souřadnice. Poté, co určíme souřadnice, opět pseudonáhodně určíme index vzorku, dosadíme hodnoty odpovídající indexu do parciálních diferenciálních rovnic a určíme souřadnice novým způsobem и atd. dokud konvergence nezezelená. Na první pohled se nemusí zdát, že by to vůbec mohlo fungovat, ale je to tak. Je pravda, že stojí za zmínku, že chyba neklesá s každým krokem, ale tendence tu určitě je.
Jaké jsou výhody stochastického gradientového sestupu oproti konvenčnímu? Pokud je velikost našeho vzorku velmi velká a měří se v desítkách tisíc hodnot, pak je mnohem jednodušší zpracovat jich, řekněme, náhodných tisíc, než celý vzorek. Zde vstupuje do hry stochastický gradientní sestup. V našem případě si samozřejmě nevšimneme velkého rozdílu.
Podívejme se na kód.
Kód pro stochastický gradient sestup
# определим функцию стох.град.шага
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])
Pozorně se podíváme na koeficienty a přistihneme se, jak si klademe otázku „Jak je to možné? Získali jsme jiné hodnoty koeficientů и . Možná stochastický gradientní sestup našel pro rovnici optimálnější parametry? Bohužel ne. Stačí se podívat na součet čtverců odchylek a vidět, že s novými hodnotami koeficientů je chyba větší. Nespěcháme k zoufalství. Vytvořme graf změny chyby.
Kód pro vynesení součtu kvadrátů odchylek v sestupu 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 „Součet druhých mocnin při sestupu stochastickým gradientem“
Při pohledu na harmonogram vše zapadá a nyní vše napravíme.
Tak, co se stalo? Stalo se následující. Když náhodně vybereme měsíc, pak se náš algoritmus snaží snížit chyby ve výpočtu příjmů pro vybraný měsíc. Poté vybereme další měsíc a výpočet zopakujeme, ale snížíme chybu pro druhý vybraný měsíc. Nyní si pamatujte, že první dva měsíce se výrazně odchylují od linie jednoduché lineární regresní rovnice. To znamená, že když je vybrán kterýkoli z těchto dvou měsíců, snížením chyby každého z nich náš algoritmus vážně zvýší chybu pro celý vzorek. Tak co dělat? Odpověď je jednoduchá: musíte snížit krok klesání. Snížením kroku sestupu totiž chyba také přestane „skákat“ nahoru a dolů. Nebo spíše chyba „skákání“ se nezastaví, ale neudělá to tak rychle :) Pojďme zkontrolovat.
Kód pro spuštění SGD s menšími přírůstky
# запустим функцию, уменьшив шаг в 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 „Součet druhých mocnin při sestupu stochastickým gradientem (80 tisíc kroků)“
Koeficienty se zlepšily, ale stále nejsou ideální. Hypoteticky to lze takto napravit. Vybíráme například v posledních 1000 iteracích hodnoty koeficientů, se kterými byla provedena minimální chyba. Je pravda, že k tomu budeme muset zapsat i hodnoty samotných koeficientů. Nebudeme to dělat, ale raději dbáme na harmonogram. Vypadá to hladce a zdá se, že chyba rovnoměrně klesá. Ve skutečnosti to není pravda. Podívejme se na prvních 1000 iterací a porovnejme je s posledními.
Kód pro graf SGD (prvních 1000 kroků)
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 „Součet druhých mocnin SGD (prvních 1000 kroků)“
Graf č. 8 „Součet druhých mocnin SGD (posledních 1000 kroků)“
Na samém začátku klesání pozorujeme celkem rovnoměrný a strmý pokles chyby. V posledních iteracích vidíme, že chyba se pohybuje kolem a kolem hodnoty 1,475 a v některých okamžicích se této optimální hodnotě dokonce rovná, ale pak stále stoupá... Opakuji, můžete si zapsat hodnoty koeficienty и a poté vyberte ty, u kterých je chyba minimální. Měli jsme však vážnější problém: museli jsme udělat 80 tisíc kroků (viz kód), abychom se dostali k hodnotám blízkým optimální. A to je již v rozporu s myšlenkou úspory výpočetního času se stochastickým sestupem gradientu vzhledem k sestupu gradientu. Co lze opravit a zlepšit? Není těžké si všimnout, že v prvních iteracích s jistotou klesáme, a proto bychom měli v prvních iteracích nechat velký krok a snižovat krok, jak postupujeme vpřed. V tomto článku to neuděláme – už je příliš dlouhý. Kdo chce, může si sám myslet, jak na to, není to těžké :)
Nyní provedeme stochastický gradientní sestup pomocí knihovny nemotorný (a nechtějme klopýtat o kameny, které jsme identifikovali dříve)
Kód pro sestup 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 se ukázaly být téměř stejné jako při sestupu bez použití nemotorný. To je však logické.
Pojďme zjistit, jak dlouho nám trvaly stochastické gradientní sestupy.
Kód pro určení doby výpočtu SGD (80 tisíc kroků)
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 dále do lesa, tím jsou mraky tmavší: nejlepší výsledek opět ukazuje „vlastně napsaný“ vzorec. To vše naznačuje, že musí existovat ještě rafinovanější způsoby využití knihovny nemotorný, které skutečně urychlují výpočetní operace. V tomto článku se o nich nedozvíme. Ve volném čase bude o čem přemýšlet :)
Shrňte
Než to shrnu, rád bych odpověděl na otázku, která s největší pravděpodobností vyvstala od našeho milého čtenáře. Proč vlastně takové „mučení“ sestupy, proč potřebujeme chodit nahoru a dolů z hory (většinou dolů), abychom našli vzácnou nížinu, máme-li v rukou tak výkonný a jednoduchý přístroj, v formou analytického řešení, které nás okamžitě teleportuje na správné místo?
Odpověď na tuto otázku leží na povrchu. Nyní jsme se podívali na velmi jednoduchý příklad, ve kterém je pravdivá odpověď záleží na jednom znamení . To se v životě často nevidí, takže si představme, že máme 2, 30, 50 nebo více znamení. Přidejme k tomu tisíce, nebo dokonce desetitisíce hodnot pro každý atribut. V tomto případě nemusí analytické řešení test obstát a selhat. Gradient sestup a jeho variace nás zase pomalu, ale jistě přiblíží k cíli – minimu funkce. A nebojte se rychlosti – pravděpodobně se podíváme na způsoby, které nám umožní nastavit a regulovat délku kroku (tedy rychlost).
A nyní již samotné krátké shrnutí.
Zaprvé doufám, že materiál uvedený v článku pomůže začínajícím „datovým vědcům“ pochopit, jak řešit jednoduché (nejen) lineární regresní rovnice.
Za druhé jsme se podívali na několik způsobů, jak rovnici vyřešit. Nyní, v závislosti na situaci, můžeme vybrat ten, který se nejlépe hodí k vyřešení problému.
Zatřetí jsme viděli sílu dalšího nastavení, konkrétně délku kroku sestupu. Tento parametr nelze zanedbat. Jak je uvedeno výše, aby se snížily náklady na výpočty, měla by se během sestupu měnit délka kroku.
Za čtvrté, v našem případě „doma psané“ funkce vykazovaly nejlepší časové výsledky pro výpočty. Je to pravděpodobně způsobeno ne nejprofesionálnějším využitím možností knihovny nemotorný. Ale budiž, následující závěr se nabízí. Jednak se někdy vyplatí zpochybňovat ustálené názory a jednak se ne vždy vyplatí vše komplikovat – naopak někdy je jednodušší způsob řešení problému efektivnější. A protože naším cílem bylo analyzovat tři přístupy k řešení jednoduché lineární regresní rovnice, stačilo nám použití „vlastně psaných“ funkcí.
Literatura (nebo něco takového)
1. Lineární regrese
2. Metoda nejmenších čtverců
3. Derivát
4. Přechod
5. Gradientní klesání
6. Knihovna NumPy
Zdroj: www.habr.com