Paglutas ng simpleng linear regression equation

Tinatalakay ng artikulo ang ilang paraan upang matukoy ang mathematical equation ng isang simpleng (pinares) na linya ng regression.

Ang lahat ng paraan ng paglutas ng equation na tinalakay dito ay batay sa least squares method. Tukuyin natin ang mga pamamaraan tulad ng sumusunod:

• Analitikal na solusyon
• Gradient Descent
• Stochastic gradient descent

Para sa bawat paraan ng paglutas ng equation ng isang tuwid na linya, ang artikulo ay nagbibigay ng iba't ibang mga function, na higit sa lahat ay nahahati sa mga nakasulat nang hindi gumagamit ng library. numpy at ang mga ginagamit para sa mga kalkulasyon numpy. Ito ay pinaniniwalaan na mahusay na paggamit numpy bawasan ang mga gastos sa pag-compute.

Lahat ng code na ibinigay sa artikulo ay nakasulat sa wika sawa 2.7 gamit Jupiter Notebook. Naka-post ang source code at file na may sample na data Github

Ang artikulo ay mas naglalayong kapwa sa mga nagsisimula at sa mga unti-unting nagsimulang makabisado ang pag-aaral ng isang napakalawak na seksyon sa artificial intelligence - machine learning.

Upang ilarawan ang materyal, gumagamit kami ng isang napakasimpleng halimbawa.

Halimbawang kondisyon

Mayroon kaming limang mga halaga na nagpapakilala sa pagtitiwala Y mula sa X (Talahanayan Blg. 1):

Talahanayan Blg. 1 "Mga halimbawang kondisyon"

Ipagpalagay namin na ang mga halaga ay ang buwan ng taon, at — kita ngayong buwan. Sa madaling salita, ang kita ay nakasalalay sa buwan ng taon, at - ang tanging palatandaan kung saan nakasalalay ang kita.

Ang halimbawa ay kaya-kaya, kapwa mula sa punto ng view ng kondisyon na pag-asa ng kita sa buwan ng taon, at mula sa punto ng view ng bilang ng mga halaga - napakakaunti sa kanila. Gayunpaman, ang gayong pagpapagaan ay gagawing posible, gaya ng sinasabi nila, na ipaliwanag, hindi palaging madali, ang materyal na pinagsasama-sama ng mga nagsisimula. At din ang pagiging simple ng mga numero ay magpapahintulot sa mga nais na malutas ang halimbawa sa papel nang walang makabuluhang gastos sa paggawa.

Ipagpalagay natin na ang pag-asa na ibinigay sa halimbawa ay maaaring matantya nang maayos sa pamamagitan ng matematikal na equation ng isang simpleng (pinares) na linya ng regression ng form:

saan ay ang buwan kung saan natanggap ang kita, — kita na naaayon sa buwan, и ay ang mga regression coefficient ng tinantyang linya.

Tandaan na ang coefficient madalas na tinatawag na slope o gradient ng tinantyang linya; kumakatawan sa halaga kung saan ang kapag nagbago ito .

Malinaw, ang aming gawain sa halimbawa ay ang pumili ng mga naturang coefficient sa equation и , kung saan ang mga paglihis ng aming kinakalkula na mga halaga ng kita ayon sa buwan mula sa mga totoong sagot, i.e. ang mga halaga na ipinakita sa sample ay magiging minimal.

Pinakamababang parisukat na pamamaraan

Ayon sa paraan ng least squares, dapat kalkulahin ang deviation sa pamamagitan ng pag-squaring nito. Ang pamamaraan na ito ay nagpapahintulot sa iyo na maiwasan ang magkaparehong pagkansela ng mga paglihis kung mayroon silang magkasalungat na mga palatandaan. Halimbawa, kung sa isang kaso, ang paglihis ay +5 (plus lima), at sa isa pa -5 (minus lima), pagkatapos ay ang kabuuan ng mga paglihis ay magkakansela sa isa't isa at magiging 0 (zero). Posibleng hindi i-square ang deviation, ngunit gamitin ang property ng modulus at pagkatapos ang lahat ng deviations ay magiging positibo at maiipon. Hindi namin tatalakayin nang detalyado ang puntong ito, ngunit ipahiwatig lamang na para sa kaginhawaan ng mga kalkulasyon, kaugalian na i-square ang paglihis.

Ito ang hitsura ng formula kung saan tutukuyin natin ang pinakamaliit na kabuuan ng mga squared deviations (mga error):

saan ay isang function ng approximation ng mga totoong sagot (iyon ay, ang kita na aming kinakalkula),

ay ang mga totoong sagot (kitang ibinigay sa sample),

ay ang sample index (bilang ng buwan kung saan tinutukoy ang paglihis)

Ibahin natin ang function, tukuyin ang mga partial differential equation, at maging handa na magpatuloy sa analytical solution. Ngunit una, magsagawa tayo ng maikling iskursiyon tungkol sa kung ano ang pagkakaiba at tandaan ang geometric na kahulugan ng hinalaw.

Differentiation

Ang differentiation ay ang operasyon ng paghahanap ng derivative ng isang function.

Ano ang gamit ng derivative? Ang derivative ng isang function ay nagpapakilala sa rate ng pagbabago ng function at nagsasabi sa amin ng direksyon nito. Kung ang derivative sa isang partikular na punto ay positibo, kung gayon ang function ay tataas; kung hindi, ang function ay bumababa. At kung mas malaki ang halaga ng absolute derivative, mas mataas ang rate ng pagbabago ng mga value ng function, gayundin ang mas matarik na slope ng function graph.

Halimbawa, sa ilalim ng mga kondisyon ng isang Cartesian coordinate system, ang halaga ng derivative sa puntong M(0,0) ay katumbas ng + 25 nangangahulugan na sa isang naibigay na punto, kapag ang halaga ay inilipat sa kanan sa pamamagitan ng isang karaniwang yunit, halaga tumaas ng 25 kumbensyonal na yunit. Sa graph, mukhang isang medyo matarik na pagtaas sa mga halaga mula sa isang naibigay na punto.

Isa pang halimbawa. Ang derivative value ay pantay -0,1 nangangahulugan na kapag inilipat bawat isang karaniwang yunit, halaga bumababa lamang ng 0,1 conventional unit. Kasabay nito, sa graph ng function, maaari nating obserbahan ang isang bahagya na kapansin-pansing pababang slope. Ang pagguhit ng isang pagkakatulad sa isang bundok, para kaming napakabagal na bumababa sa isang banayad na dalisdis mula sa isang bundok, hindi tulad ng nakaraang halimbawa, kung saan kailangan naming umakyat sa mga napakatarik na taluktok :)

Kaya, pagkatapos ng pagkakaiba-iba ng function sa pamamagitan ng logro и , tinutukoy namin ang 1st order na partial differential equation. Matapos matukoy ang mga equation, makakatanggap kami ng isang sistema ng dalawang equation, sa pamamagitan ng paglutas kung saan magagawa naming pumili ng mga naturang halaga ng mga coefficient и , kung saan ang mga halaga ng kaukulang mga derivative sa mga ibinigay na punto ay nagbabago ng napakaliit na halaga, at sa kaso ng isang analytical na solusyon ay hindi nagbabago. Sa madaling salita, ang pag-andar ng error sa mga nahanap na coefficient ay aabot sa isang minimum, dahil ang mga halaga ng mga partial derivatives sa mga puntong ito ay magiging katumbas ng zero.

Kaya, ayon sa mga patakaran ng pagkita ng kaibhan, ang bahagyang derivative equation ng 1st order na may paggalang sa coefficient kukuha ng form:

1st order partial derivative equation na may kinalaman sa kukuha ng form:

Bilang resulta, nakatanggap kami ng isang sistema ng mga equation na may medyo simpleng analytical na solusyon:

simulan{equation*}
magsimula{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
wakas{cases}
wakas{equation*}

Bago lutasin ang equation, mag-preload tayo, suriin kung tama ang paglo-load, at i-format ang data.

Naglo-load at nag-format ng data

Dapat tandaan na dahil sa katotohanan na para sa analytical na solusyon, at pagkatapos ay para sa gradient at stochastic gradient descent, gagamitin namin ang code sa dalawang variation: gamit ang library numpy at nang hindi ginagamit ito, kakailanganin namin ang naaangkop na pag-format ng data (tingnan ang code).

Paglo-load ng data at pagpoproseso ng code

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

Pagpapakita

Ngayon, pagkatapos namin, una, na-load ang data, pangalawa, nasuri ang tama ng pag-load at sa wakas ay na-format ang data, isasagawa namin ang unang visualization. Ang pamamaraang kadalasang ginagamit para dito ay pairplot mga aklatan taga-dagat. Sa aming halimbawa, dahil sa limitadong bilang, walang saysay ang paggamit ng aklatan taga-dagat. Gagamitin namin ang regular na aklatan matplotlib at tingnan mo na lang ang scatterplot.

Scatterplot code

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

Tsart Blg. 1 "Pag-asa ng kita sa buwan ng taon"

Analitikal na solusyon

Gamitin natin ang pinakakaraniwang tool sa python at lutasin ang sistema ng mga equation:

simulan{equation*}
magsimula{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
wakas{cases}
wakas{equation*}

Ayon sa panuntunan ni Cramer makikita natin ang pangkalahatang determinant, gayundin ang mga determinant sa pamamagitan ng at , pagkatapos nito, hinahati ang determinant sa pamamagitan ng sa pangkalahatang determinant - hanapin ang koepisyent , katulad na nakikita natin ang koepisyent .

Analytical solution code

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

Narito ang nakuha namin:

Kaya, ang mga halaga ng mga coefficient ay natagpuan, ang kabuuan ng mga squared deviations ay naitatag. Gumuhit tayo ng isang tuwid na linya sa scattering histogram alinsunod sa mga nakitang coefficient.

Code ng linya ng regression

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

Tsart Blg. 2 "Tama at kalkuladong mga sagot"

Maaari mong tingnan ang deviation graph para sa bawat buwan. Sa aming kaso, hindi kami makakakuha ng anumang makabuluhang praktikal na halaga mula dito, ngunit matutugunan namin ang aming pag-usisa tungkol sa kung gaano kahusay na nailalarawan ng simpleng linear regression equation ang pagtitiwala ng kita sa buwan ng taon.

Code ng paglihis ng tsart

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

Tsart Blg. 3 "Mga Paglihis, %"

Hindi perpekto, ngunit natapos namin ang aming gawain.

Sumulat tayo ng isang function na, upang matukoy ang mga coefficient и gumagamit ng aklatan numpy, mas tiyak, magsusulat kami ng dalawang function: ang isa ay gumagamit ng isang pseudoinverse matrix (hindi inirerekomenda sa pagsasanay, dahil ang proseso ay computationally kumplikado at hindi matatag), ang isa ay gumagamit ng isang matrix equation.

Analytical Solution Code (NumPy)

# для начала добавим столбец с не изменяющимся значением в 1. 
# Данный столбец нужен для того, чтобы не обрабатывать отдельно коэффицент a
vector_1 = np.ones((x_np.shape[0],1))
x_np = table_zero[['x']].values # на всякий случай приведем в первичный формат вектор x_np
x_np = np.hstack((vector_1,x_np))

# проверим то, что все сделали правильно
print vector_1[0:3]
print x_np[0:3]
print '***************************************'
print

# напишем функцию, которая определяет значения коэффициентов a и b с использованием псевдообратной матрицы
def pseudoinverse_matrix(X, y):
    # задаем явный формат матрицы признаков
    X = np.matrix(X)
    # определяем транспонированную матрицу
    XT = X.T
    # определяем квадратную матрицу
    XTX = XT*X
    # определяем псевдообратную матрицу
    inv = np.linalg.pinv(XTX)
    # задаем явный формат матрицы ответов
    y = np.matrix(y)
    # находим вектор весов
    return (inv*XT)*y

# запустим функцию
ab_np = pseudoinverse_matrix(x_np, y_np)
print ab_np
print '***************************************'
print

# напишем функцию, которая использует для решения матричное уравнение
def matrix_equation(X,y):
    a = np.dot(X.T, X)
    b = np.dot(X.T, y)
    return np.linalg.solve(a, b)

# запустим функцию
ab_np = matrix_equation(x_np,y_np)
print ab_np

Ihambing natin ang oras na ginugol sa pagtukoy ng mga koepisyent и , alinsunod sa 3 ipinakitang pamamaraan.

Code para sa pagkalkula ng oras ng pagkalkula

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)

Sa kaunting data, lalabas ang isang function na "self-written", na hahanapin ang mga coefficient gamit ang paraan ng Cramer.

Ngayon ay maaari kang lumipat sa iba pang mga paraan upang makahanap ng mga coefficient и .

Gradient Descent

Una, tukuyin natin kung ano ang gradient. Sa madaling salita, ang gradient ay isang segment na nagpapahiwatig ng direksyon ng maximum na paglaki ng isang function. Sa pamamagitan ng pagkakatulad sa pag-akyat sa isang bundok, kung saan nakaharap ang gradient ay kung saan ang pinakamatarik na pag-akyat sa tuktok ng bundok. Ang pagbuo ng halimbawa sa bundok, naaalala natin na sa katunayan kailangan natin ang pinakamatarik na pagbaba upang maabot ang mababang lupain sa lalong madaling panahon, iyon ay, ang pinakamaliit - ang lugar kung saan ang paggana ay hindi tumataas o bumababa. Sa puntong ito ang derivative ay magiging katumbas ng zero. Samakatuwid, hindi namin kailangan ng isang gradient, ngunit isang antigradient. Upang mahanap ang antigradient kailangan mo lamang i-multiply ang gradient sa pamamagitan ng -1 (minus one).

Bigyang-pansin natin ang katotohanan na ang isang function ay maaaring magkaroon ng ilang minima, at kapag bumaba sa isa sa mga ito gamit ang algorithm na iminungkahi sa ibaba, hindi na tayo makakahanap ng isa pang minimum, na maaaring mas mababa kaysa sa nakita. Magpahinga tayo, hindi ito banta sa atin! Sa aming kaso kami ay nakikitungo sa isang solong minimum, dahil ang aming function sa graph ay isang regular na parabola. At gaya ng dapat alam nating lahat mula sa ating kursong matematika sa paaralan, ang isang parabola ay may isang minimum lamang.

Matapos naming malaman kung bakit kailangan namin ng isang gradient, at gayundin na ang gradient ay isang segment, iyon ay, isang vector na may ibinigay na mga coordinate, na tiyak na parehong mga coefficient. и maaari nating ipatupad ang gradient descent.

Bago magsimula, iminumungkahi kong basahin lamang ang ilang pangungusap tungkol sa descent algorithm:

• Tinutukoy namin sa isang pseudo-random na paraan ang mga coefficient ng mga coefficient и . Sa aming halimbawa, tutukuyin namin ang mga coefficient malapit sa zero. Ito ay isang karaniwang kasanayan, ngunit ang bawat kaso ay maaaring may sariling kasanayan.
• Mula sa coordinate ibawas ang halaga ng 1st order partial derivative sa punto . Kaya, kung ang derivative ay positibo, ang function ay tumataas. Samakatuwid, sa pamamagitan ng pagbabawas ng halaga ng derivative, lilipat tayo sa kabaligtaran na direksyon ng paglago, iyon ay, sa direksyon ng pagbaba. Kung ang derivative ay negatibo, ang function sa puntong ito ay bumababa at sa pamamagitan ng pagbabawas ng halaga ng derivative ay lumilipat tayo sa direksyon ng pagbaba.
• Nagsasagawa kami ng katulad na operasyon kasama ang coordinate : ibawas ang halaga ng partial derivative sa punto .
• Upang hindi tumalon sa pinakamababa at lumipad sa malalim na espasyo, kinakailangan upang itakda ang laki ng hakbang sa direksyon ng pagbaba. Sa pangkalahatan, maaari kang magsulat ng isang buong artikulo tungkol sa kung paano itakda nang tama ang hakbang at kung paano ito baguhin sa panahon ng proseso ng pagbaba upang mabawasan ang mga gastos sa pagkalkula. Ngunit ngayon mayroon kaming isang bahagyang naiibang gawain sa unahan namin, at itatatag namin ang laki ng hakbang gamit ang siyentipikong pamamaraan ng "poke" o, tulad ng sinasabi nila sa karaniwang parlance, empirically.
• Sa sandaling kami ay mula sa ibinigay na mga coordinate и ibawas ang mga halaga ng mga derivatives, nakakakuha kami ng mga bagong coordinate и . Ginagawa namin ang susunod na hakbang (pagbabawas), na mula sa kinakalkula na mga coordinate. At kaya ang pag-ikot ay nagsisimula nang paulit-ulit, hanggang sa makamit ang kinakailangang tagpo.

Lahat! Ngayon ay handa na kaming maghanap ng pinakamalalim na bangin ng Mariana Trench. Magsimula na tayo.

Code para sa gradient descent

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

Sumisid kami sa pinakailalim ng Mariana Trench at doon nakita namin ang lahat ng parehong halaga ng coefficient и , na kung ano mismo ang inaasahan.

Sumisid ulit tayo, sa pagkakataong ito, mapupuno na ng iba pang teknolohiya ang ating deep-sea vehicle, katulad ng library. numpy.

Code para sa gradient descent (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

Mga halaga ng koepisyent и hindi nababago.

Tingnan natin kung paano nagbago ang error sa panahon ng gradient descent, ibig sabihin, kung paano nagbago ang kabuuan ng mga squared deviations sa bawat hakbang.

Code para sa pag-plot ng mga kabuuan ng mga squared deviations

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

Graph No. 4 "Kabuuan ng mga squared deviation sa panahon ng gradient descent"

Sa graph nakita namin na sa bawat hakbang ay bumababa ang error, at pagkatapos ng isang tiyak na bilang ng mga pag-ulit ay naobserbahan namin ang halos pahalang na linya.

Panghuli, tantyahin natin ang pagkakaiba sa oras ng pagpapatupad ng code:

Code upang matukoy ang oras ng pagkalkula ng gradient descent

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)

Marahil ay may ginagawa tayong mali, ngunit muli ito ay isang simpleng "home-written" na function na hindi gumagamit ng library numpy higit ang pagganap sa oras ng pagkalkula ng isang function gamit ang library numpy.

Ngunit hindi kami nakatayo, ngunit lumilipat patungo sa pag-aaral ng isa pang kapana-panabik na paraan upang malutas ang simpleng linear regression equation. Kilalanin kami!

Stochastic gradient descent

Upang mabilis na maunawaan ang prinsipyo ng pagpapatakbo ng stochastic gradient descent, mas mahusay na matukoy ang mga pagkakaiba nito mula sa ordinaryong gradient descent. Kami, sa kaso ng gradient descent, sa mga equation ng derivatives ng и ginamit ang mga kabuuan ng mga halaga ng lahat ng mga tampok at mga totoong sagot na magagamit sa sample (iyon ay, ang mga kabuuan ng lahat и ). Sa stochastic gradient descent, hindi namin gagamitin ang lahat ng value na naroroon sa sample, ngunit sa halip, pseudo-random na piliin ang tinatawag na sample index at gamitin ang mga value nito.

Halimbawa, kung ang index ay tinutukoy na numero 3 (tatlo), pagkatapos ay kukunin namin ang mga halaga и , pagkatapos ay pinapalitan namin ang mga halaga sa mga derivative equation at tinutukoy ang mga bagong coordinate. Pagkatapos, nang matukoy ang mga coordinate, muli naming pseudo-random na tinutukoy ang sample index, pinapalitan ang mga halaga na tumutugma sa index sa mga partial differential equation, at tinutukoy ang mga coordinate sa isang bagong paraan и atbp. hanggang sa maging berde ang convergence. Sa unang sulyap, maaaring hindi ito maaaring gumana, ngunit ito ay gumagana. Ito ay totoo na ito ay nagkakahalaga ng noting na ang error ay hindi bumababa sa bawat hakbang, ngunit may tiyak na isang ugali.

Ano ang mga pakinabang ng stochastic gradient descent kumpara sa conventional? Kung ang aming sample size ay napakalaki at sinusukat sa sampu-sampung libong mga halaga, kung gayon ito ay mas madaling iproseso, halimbawa, isang random na libo sa mga ito, kaysa sa buong sample. Dito pumapasok ang stochastic gradient descent. Sa aming kaso, siyempre, hindi namin mapapansin ang malaking pagkakaiba.

Tingnan natin ang code.

Code para sa stochastic gradient descent

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

Tinitingnan naming mabuti ang mga coefficient at nahuhuli namin ang aming sarili na nagtatanong ng tanong na "Paano ito mangyayari?" Nakakuha kami ng iba pang mga halaga ng coefficient и . Siguro ang stochastic gradient descent ay nakahanap ng mas pinakamainam na mga parameter para sa equation? Sa kasamaang palad hindi. Ito ay sapat na upang tingnan ang kabuuan ng mga squared deviations at makita na sa mga bagong halaga ng mga coefficient, ang error ay mas malaki. Hindi kami nagmamadaling mawalan ng pag-asa. Bumuo tayo ng graph ng pagbabago ng error.

Code para sa pag-plot ng kabuuan ng mga squared deviation sa stochastic gradient descent

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

Graph No. 5 "Kabuuan ng mga squared deviation sa panahon ng stochastic gradient descent"

Sa pagtingin sa iskedyul, ang lahat ay nahuhulog sa lugar at ngayon ay aayusin namin ang lahat.

So anong nangyari? Ang mga sumusunod ay nangyari. Kapag random kaming pumili ng isang buwan, para sa napiling buwan na hinahangad ng aming algorithm na bawasan ang error sa pagkalkula ng kita. Pagkatapos ay pumili kami ng isa pang buwan at ulitin ang pagkalkula, ngunit binabawasan namin ang error para sa ikalawang napiling buwan. Ngayon tandaan na ang unang dalawang buwan ay lumihis nang malaki mula sa linya ng simpleng linear regression equation. Nangangahulugan ito na kapag napili ang alinman sa dalawang buwang ito, sa pamamagitan ng pagbabawas ng error ng bawat isa sa kanila, seryosong pinapataas ng aming algorithm ang error para sa buong sample. Ano ang gagawin? Ang sagot ay simple: kailangan mong bawasan ang hakbang ng pagbaba. Pagkatapos ng lahat, sa pamamagitan ng pagbabawas ng hakbang ng pagbaba, ang error ay titigil din sa "paglukso" pataas at pababa. O sa halip, ang error sa "paglukso" ay hindi titigil, ngunit hindi ito magagawa nang napakabilis :) Suriin natin.

Code para patakbuhin ang SGD na may mas maliliit na increment

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

Graph No. 6 "Kabuuan ng mga squared deviations sa panahon ng stochastic gradient descent (80 thousand steps)"

Ang mga coefficient ay bumuti, ngunit hindi pa rin perpekto. Hypothetically, ito ay maaaring itama sa ganitong paraan. Pinipili namin, halimbawa, sa huling 1000 na pag-ulit ang mga halaga ng mga coefficient kung saan nagawa ang pinakamababang error. Totoo, para dito kailangan din nating isulat ang mga halaga ng mga coefficient mismo. Hindi namin ito gagawin, bagkus ay bigyang pansin ang iskedyul. Mukhang makinis at ang error ay tila bumaba nang pantay. Sa totoo lang hindi ito totoo. Tingnan natin ang unang 1000 pag-ulit at ihambing ang mga ito sa huli.

Code para sa SGD chart (unang 1000 hakbang)

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

Graph No. 7 "Kabuuan ng mga squared deviations SGD (unang 1000 hakbang)"

Graph No. 8 "Kabuuan ng mga squared deviations SGD (huling 1000 hakbang)"

Sa pinakadulo simula ng pagbaba, napansin namin ang isang medyo pare-pareho at matarik na pagbaba ng error. Sa mga huling pag-ulit, nakikita natin na ang error ay umiikot sa halagang 1,475 at sa ilang sandali ay katumbas ng pinakamainam na halaga na ito, ngunit pagkatapos ay tumataas pa rin ito... Uulitin ko, maaari mong isulat ang mga halaga ng coefficients и , at pagkatapos ay piliin ang mga kung saan ang error ay minimal. Gayunpaman, nagkaroon kami ng mas malubhang problema: kinailangan naming gumawa ng 80 libong hakbang (tingnan ang code) upang makuha ang mga halaga na malapit sa pinakamainam. At ito ay sumasalungat sa ideya ng pag-save ng oras ng pagkalkula gamit ang stochastic gradient descent na may kaugnayan sa gradient descent. Ano ang maaaring itama at pagbutihin? Hindi mahirap mapansin na sa mga unang pag-ulit ay kumpiyansa tayong bababa at, samakatuwid, dapat tayong mag-iwan ng malaking hakbang sa mga unang pag-ulit at bawasan ang hakbang habang sumusulong tayo. Hindi namin ito gagawin sa artikulong ito - ito ay masyadong mahaba. Ang mga nagnanais ay maaaring mag-isip para sa kanilang sarili kung paano gawin ito, hindi ito mahirap :)

Ngayon, gawin natin ang stochastic gradient descent gamit ang library numpy (at huwag tayong matisod sa mga batong natukoy natin kanina)

Code para sa Stochastic Gradient Descent (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

Ang mga halaga ay naging halos kapareho ng kapag bumaba nang hindi gumagamit numpy. Gayunpaman, ito ay lohikal.

Alamin natin kung gaano katagal tayo inabot ng stochastic gradient descents.

Code para sa pagtukoy ng oras ng pagkalkula ng SGD (80 libong hakbang)

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)

Ang mas malayo sa kagubatan, mas madilim ang mga ulap: muli, ang "self-written" na formula ay nagpapakita ng pinakamahusay na resulta. Ang lahat ng ito ay nagmumungkahi na dapat mayroong mas banayad na mga paraan upang magamit ang aklatan numpy, na talagang nagpapabilis sa mga pagpapatakbo ng pagkalkula. Sa artikulong ito hindi natin malalaman ang tungkol sa kanila. Mayroong isang bagay na dapat isipin sa iyong libreng oras :)

Ibigay namin ang buod

Bago buod, nais kong sagutin ang isang katanungan na malamang na lumitaw mula sa aming mahal na mambabasa. Bakit, sa katunayan, ang gayong "pahirap" na may pagbaba, bakit kailangan nating maglakad pataas at pababa ng bundok (karamihan pababa) upang mahanap ang treasured lowland, kung nasa ating mga kamay ang napakalakas at simpleng kagamitan, sa anyo ng isang analytical na solusyon, na agad na nag-teleport sa atin sa Tamang lugar?

Ang sagot sa tanong na ito ay nasa ibabaw. Ngayon ay tumingin kami sa isang napakasimpleng halimbawa, kung saan ang tunay na sagot ay depende sa isang sign . Hindi mo ito madalas makita sa buhay, kaya isipin natin na mayroon tayong 2, 30, 50 o higit pang mga palatandaan. Dagdagan natin ito ng libu-libo, o kahit sampu-sampung libong halaga para sa bawat katangian. Sa kasong ito, ang analytical na solusyon ay maaaring hindi makatiis sa pagsubok at mabigo. Sa turn, ang gradient descent at ang mga variation nito ay dahan-dahan ngunit tiyak na maglalapit sa atin sa layunin - ang minimum ng function. At huwag mag-alala tungkol sa bilis - malamang na titingnan namin ang mga paraan na magbibigay-daan sa amin na itakda at ayusin ang haba ng hakbang (iyon ay, bilis).

At ngayon ang aktwal na maikling buod.

Una, umaasa ako na ang materyal na ipinakita sa artikulo ay makakatulong sa pagsisimula ng "mga data scientist" sa pag-unawa kung paano lutasin ang mga simple (at hindi lamang) mga linear regression equation.

Pangalawa, tumingin kami sa ilang mga paraan upang malutas ang equation. Ngayon, depende sa sitwasyon, maaari nating piliin ang isa na pinakaangkop upang malutas ang problema.

Pangatlo, nakita namin ang kapangyarihan ng mga karagdagang setting, katulad ng gradient descent step length. Ang parameter na ito ay hindi maaaring pabayaan. Tulad ng nabanggit sa itaas, upang mabawasan ang halaga ng mga kalkulasyon, dapat baguhin ang haba ng hakbang sa panahon ng pagbaba.

Pang-apat, sa aming kaso, ang mga function na "isinulat sa bahay" ay nagpakita ng pinakamahusay na mga resulta ng oras para sa mga kalkulasyon. Ito ay marahil dahil sa hindi ang pinaka-propesyonal na paggamit ng mga kakayahan ng library numpy. Ngunit maging iyon man, ang sumusunod na konklusyon ay nagmumungkahi mismo. Sa isang banda, kung minsan ito ay nagkakahalaga ng pagtatanong sa mga itinatag na opinyon, at sa kabilang banda, hindi palaging sulit na kumplikado ang lahat - sa kabaligtaran, kung minsan ang isang mas simpleng paraan ng paglutas ng isang problema ay mas epektibo. At dahil ang aming layunin ay pag-aralan ang tatlong mga diskarte sa paglutas ng isang simpleng linear regression equation, ang paggamit ng "self-written" na mga function ay sapat na para sa amin.

Panitikan (o isang katulad nito)

1. Linear regression

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

2. Pamamaraan ng hindi bababa sa mga parisukat

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Derivative

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Gradient descent

habr.com/en/post/471458

habr.com/en/post/307312

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

6. NumPy library

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

Pinagmulan: www.habr.com