🥇Vienkāršas lineārās regresijas vienādojuma atrisināšana

Rakstā aplūkoti vairāki veidi, kā noteikt vienkāršas (pāra) regresijas taisnes matemātisko vienādojumu.

Visas šeit aplūkotās vienādojuma risināšanas metodes ir balstītas uz mazāko kvadrātu metodi. Apzīmēsim metodes šādi:

Analītisks risinājums
Gradienta nolaišanās
Stohastiskā gradienta nolaišanās

Katrai taisnes vienādojuma risināšanas metodei rakstā ir sniegtas dažādas funkcijas, kuras galvenokārt tiek sadalītas tajās, kuras tiek rakstītas, neizmantojot bibliotēku nejutīgs un tie, kurus izmanto aprēķiniem nejutīgs. Tiek uzskatīts, ka prasmīga izmantošana nejutīgs samazinās skaitļošanas izmaksas.

Viss rakstā norādītais kods ir rakstīts valodā python 2.7 izmantojot Jupyter piezīmju grāmatiņa. Avota kods un fails ar datu paraugiem ir ievietoti vietnē Github

Raksts vairāk paredzēts gan iesācējiem, gan tiem, kuri jau pamazām sākuši apgūt ļoti plašu mākslīgā intelekta sadaļu – mašīnmācīšanos.

Lai ilustrētu materiālu, mēs izmantojam ļoti vienkāršu piemēru.

Piemēri nosacījumi

Mums ir piecas vērtības, kas raksturo atkarību Y no X (Tabula Nr. 1):

Tabula Nr.1 “Piemēri nosacījumi”

Mēs pieņemsim, ka vērtības ir gada mēnesis, un — ieņēmumi šajā mēnesī. Citiem vārdiem sakot, ieņēmumi ir atkarīgi no gada mēneša un - vienīgā zīme, no kuras atkarīgi ieņēmumi.

Piemērs ir tik un tā, gan no ieņēmumu nosacītās atkarības viedokļa no gada mēneša, gan no vērtību skaita viedokļa - to ir ļoti maz. Tomēr šāds vienkāršojums ļaus, kā saka, ne vienmēr viegli izskaidrot materiālu, ko iesācēji asimilē. Un arī skaitļu vienkāršība ļaus tiem, kas vēlas, atrisināt piemēru uz papīra bez ievērojamām darbaspēka izmaksām.

Pieņemsim, ka piemērā doto atkarību var diezgan labi tuvināt ar formas vienkāršas (pāra) regresijas taisnes matemātisko vienādojumu:

kur ir mēnesis, kurā tika saņemti ieņēmumi, — ieņēmumi, kas atbilst mēnesim, и ir aprēķinātās līnijas regresijas koeficienti.

Ņemiet vērā, ka koeficients bieži sauc par aplēstās līnijas slīpumu vai gradientu; apzīmē summu, par kādu kad tas mainās .

Acīmredzot mūsu uzdevums piemērā ir izvēlēties šādus koeficientus vienādojumā и , kurā mūsu aprēķināto ieņēmumu vērtību novirzes pa mēnešiem no patiesajām atbildēm, t.i. paraugā uzrādītās vērtības būs minimālas.

Mazākā kvadrāta metode

Saskaņā ar mazāko kvadrātu metodi novirze jāaprēķina, sadalot to kvadrātā. Šis paņēmiens ļauj izvairīties no savstarpējas noviržu atcelšanas, ja tām ir pretējas pazīmes. Piemēram, ja vienā gadījumā novirze ir +5 (plus pieci), un otrā -5 (mīnus pieci), tad noviržu summa viena otru izslēgs un sasniegs 0 (nulle). Novirzi var nelikt kvadrātā, bet izmantot moduļa īpašību un tad visas novirzes būs pozitīvas un uzkrāsies. Mēs nekavēsimies pie šī punkta sīkāk, bet vienkārši norādīsim, ka aprēķinu ērtībai ir ierasts novirzes kvadrātā.

Šādi izskatās formula, ar kuru mēs noteiksim mazāko noviržu (kļūdu) kvadrātu summu:

kur ir patieso atbilžu (tas ir, mūsu aprēķināto ieņēmumu) tuvinājuma funkcija,

vai atbildes ir patiesas (ieņēmumi norādīti izlasē),

ir izlases indekss (tā mēneša numurs, kurā tiek noteikta novirze)

Atšķirsim funkciju, definēsim daļējos diferenciālvienādojumus un būsim gatavi pāriet uz analītisko risinājumu. Bet vispirms veiksim īsu ekskursiju par to, kas ir diferencēšana, un atcerēsimies atvasinājuma ģeometrisko nozīmi.

Diferencēšana

Diferenciācija ir funkcija, lai atrastu funkcijas atvasinājumu.

Kam tiek izmantots atvasinājums? Funkcijas atvasinājums raksturo funkcijas izmaiņu ātrumu un norāda tās virzienu. Ja atvasinājums noteiktā punktā ir pozitīvs, funkcija palielinās, pretējā gadījumā funkcija samazinās. Un jo lielāka ir absolūtā atvasinājuma vērtība, jo lielāks ir funkciju vērtību izmaiņu ātrums, kā arī stāvāks ir funkcijas grafika slīpums.

Piemēram, Dekarta koordinātu sistēmas apstākļos atvasinājuma vērtība punktā M(0,0) ir vienāda ar +25 nozīmē, ka noteiktā punktā, kad vērtība tiek nobīdīta pa labi ar nosacītu vienību, vērtību palielinās par 25 parastajām vienībām. Diagrammā tas izskatās kā diezgan straujš vērtību pieaugums no dotā punkta.

Vēl viens piemērs. Atvasinātā vērtība ir vienāda -0,1 nozīmē, ka tad, kad tas ir pārvietots uz vienu nosacīto vienību, vērtību samazinās tikai par 0,1 parasto vienību. Tajā pašā laikā funkcijas grafikā mēs varam novērot tikko pamanāmu lejupvērstu slīpumu. Velkot analoģiju ar kalnu, ir tā, it kā mēs ļoti lēni lejam no kalna pa lēzenu nogāzi, atšķirībā no iepriekšējā piemēra, kur mums bija jākāpj ļoti stāvās virsotnēs :)

Tādējādi pēc funkcijas diferencēšanas pēc izredzēm и , mēs definējam 1. kārtas daļējos diferenciālvienādojumus. Pēc vienādojumu noteikšanas iegūsim divu vienādojumu sistēmu, kuru risinot varēsim izvēlēties šādas koeficientu vērtības и , kurā attiecīgo atvasinājumu vērtības noteiktos punktos mainās ļoti, ļoti maz, un analītiskā risinājuma gadījumā nemainās vispār. Citiem vārdiem sakot, kļūdas funkcija pie atrastajiem koeficientiem sasniegs minimumu, jo daļējo atvasinājumu vērtības šajos punktos būs vienādas ar nulli.

Tātad saskaņā ar diferenciācijas noteikumiem pirmās kārtas daļējais atvasinātais vienādojums attiecībā pret koeficientu būs šādā formā:

1. kārtas daļēja atvasinājuma vienādojums attiecībā uz būs šādā formā:

Rezultātā mēs saņēmām vienādojumu sistēmu, kurai ir diezgan vienkāršs analītisks risinājums:

sākt{vienādojums*}
sākt{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
beigas{gadījumi}
beigu{vienādojums*}

Pirms vienādojuma atrisināšanas veiksim priekšielādēšanu, pārbaudīsim, vai ielāde ir pareiza, un formatēsim datus.

Datu ielāde un formatēšana

Jāatzīmē, ka, ņemot vērā to, ka analītiskajam risinājumam un pēc tam gradienta un stohastiskā gradienta nolaišanās kodu izmantosim divās variācijās: izmantojot bibliotēku. nejutīgs un neizmantojot to, tad mums būs nepieciešams atbilstošs datu formatējums (skat. kodu).

Datu ielādes un apstrādes kods

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

Tagad, kad, pirmkārt, būsim ielādējuši datus, otrkārt, pārbaudījuši ielādes pareizību un beidzot formatējuši datus, mēs veiksim pirmo vizualizāciju. Šim nolūkam bieži izmantotā metode ir pāru gabals bibliotēkas Jūras rags. Mūsu piemērā ierobežotā skaita dēļ bibliotēku izmantot nav jēgas Jūras rags. Mēs izmantosim parasto bibliotēku matplotlib un vienkārši paskaties uz izkliedes diagrammu.

Izkliedes kods

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

Diagramma Nr.1 “Ieņēmumu atkarība no gada mēneša”

Analītisks risinājums

Izmantosim visizplatītākos rīkus pitons un atrisiniet vienādojumu sistēmu:

sākt{vienādojums*}
sākt{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
beigas{gadījumi}
beigu{vienādojums*}

Saskaņā ar Krāmera likumu mēs atradīsim vispārīgo noteicēju, kā arī determinantus pēc un pēc , pēc kura determinantu dalot ar uz vispārīgo determinantu - atrodiet koeficientu , līdzīgi mēs atrodam koeficientu .

Analītiskā risinājuma kods

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

Lūk, ko mēs saņēmām:

Tātad ir atrastas koeficientu vērtības, ir noteikta noviržu kvadrātā summa. Izkliedes histogrammā zīmēsim taisnu līniju atbilstoši atrastajiem koeficientiem.

Regresijas līnijas kods

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

Diagramma Nr.2 “Pareizas un aprēķinātas atbildes”

Varat apskatīt katra mēneša noviržu grafiku. Mūsu gadījumā mēs no tā neiegūsim būtisku praktisku vērtību, bet apmierināsim savu ziņkāri par to, cik labi vienkāršais lineārās regresijas vienādojums raksturo ieņēmumu atkarību no gada mēneša.

Noviržu diagrammas kods

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

Diagramma Nr. 3 “Novirzes, %”

Nav ideāls, bet mēs paveicām savu uzdevumu.

Uzrakstīsim funkciju, kas, lai noteiktu koeficientus и izmanto bibliotēku nejutīgs, precīzāk, rakstīsim divas funkcijas: vienu izmantojot pseidoinverso matricu (praksē nav ieteicams, jo process ir skaitļošanas ziņā sarežģīts un nestabils), otru izmantojot matricas vienādojumu.

Analītiskā risinājuma kods (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

Salīdzināsim koeficientu noteikšanai patērēto laiku и , saskaņā ar 3 uzrādītajām metodēm.

Kods aprēķina laika aprēķināšanai

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)

Ar nelielu datu apjomu priekšā iznāk “pašrakstīta” funkcija, kas atrod koeficientus, izmantojot Krāmera metodi.

Tagad varat pāriet uz citiem veidiem, kā atrast koeficientus и .

Gradienta nolaišanās

Vispirms definēsim, kas ir gradients. Vienkārši sakot, gradients ir segments, kas norāda funkcijas maksimālās pieauguma virzienu. Pēc analoģijas ar kāpšanu kalnā, kur slīpuma sejas ir tur, kur ir stāvākais kāpums kalna virsotnē. Izstrādājot piemēru ar kalnu, mēs atceramies, ka patiesībā mums ir nepieciešams stāvākais nobrauciens, lai pēc iespējas ātrāk sasniegtu zemieni, tas ir, minimums - vieta, kur funkcija nepalielinās vai nesamazinās. Šajā brīdī atvasinājums būs vienāds ar nulli. Tāpēc mums nav vajadzīgs gradients, bet gan antigradients. Lai atrastu antigradientu, jums vienkārši jāreizina gradients ar -1 (mīnus viens).

Pievērsīsim uzmanību tam, ka funkcijai var būt vairāki minimumi, un, nolaižoties vienā no tiem, izmantojot tālāk piedāvāto algoritmu, mēs nevarēsim atrast citu minimumu, kas var būt zemāks par atrasto. Atpūšamies, tas mums nedraud! Mūsu gadījumā mums ir darīšana ar vienu minimumu, jo mūsu funkcija grafikā ir regulāra parabola. Un kā mums visiem ļoti labi vajadzētu zināt no mūsu skolas matemātikas kursa, parabolai ir tikai viens minimums.

Pēc tam, kad mēs uzzinājām, kāpēc mums vajadzīgs gradients, kā arī to, ka gradients ir segments, tas ir, vektors ar noteiktām koordinātām, kas ir tieši tādi paši koeficienti и mēs varam īstenot gradienta nolaišanos.

Pirms darba sākšanas iesaku izlasīt tikai dažus teikumus par nolaišanās algoritmu:

Mēs pseidogadījuma veidā nosakām koeficientu koordinātas и . Mūsu piemērā mēs noteiksim koeficientus, kas ir tuvu nullei. Tā ir izplatīta prakse, taču katram gadījumam var būt sava prakse.
No koordinātas punktā atņem 1. kārtas daļējā atvasinājuma vērtību . Tātad, ja atvasinājums ir pozitīvs, funkcija palielinās. Tāpēc, atņemot atvasinājuma vērtību, mēs virzīsimies pretējā pieauguma virzienā, tas ir, lejupslīdes virzienā. Ja atvasinājums ir negatīvs, tad funkcija šajā punktā samazinās un, atņemot atvasinājuma vērtību, mēs virzāmies lejupešanas virzienā.
Mēs veicam līdzīgu darbību ar koordinātu : atņem daļējā atvasinājuma vērtību punktā .
Lai nepārlēktu pāri minimumam un neaizlidotu dziļā kosmosā, ir jāiestata soļa lielums nolaišanās virzienā. Kopumā jūs varētu uzrakstīt veselu rakstu par to, kā pareizi iestatīt soli un kā to mainīt nolaišanās procesa laikā, lai samazinātu skaitļošanas izmaksas. Bet tagad mūs gaida nedaudz cits uzdevums, un soļa lielumu noteiksim, izmantojot zinātnisko metodi “poke” jeb, kā mēdz teikt, empīriski.
Kad esam no dotajām koordinātēm и atņemot atvasinājumu vērtības, mēs iegūstam jaunas koordinātas и . Mēs veicam nākamo soli (atņemšanu), jau no aprēķinātajām koordinātām. Un tā cikls sākas atkal un atkal, līdz tiek sasniegta vajadzīgā konverģence.

Visi! Tagad esam gatavi doties meklēt Marianas tranšejas dziļāko aizu. Sāksim.

Gradienta nolaišanās kods

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

Mēs nirām līdz pašai Marianas tranšejas apakšai, un tur atradām visas tās pašas koeficientu vērtības и , kas bija tieši tas, kas bija gaidāms.

Nirsim vēlreiz, tikai šoreiz mūsu dziļjūras transportlīdzeklis būs piepildīts ar citām tehnoloģijām, proti, bibliotēku nejutīgs.

Gradienta nolaišanās kods (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

Koeficientu vērtības и nemaināms.

Apskatīsim, kā mainījās kļūda gradienta nolaišanās laikā, tas ir, kā mainījās noviržu summa kvadrātā ar katru soli.

Kods kvadrātu noviržu summu uzzīmēšanai

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

Diagramma Nr. 4 “Noviržu summa kvadrātā gradienta nolaišanās laikā”

Grafikā redzam, ka ar katru soli kļūda samazinās, un pēc noteikta iterāciju skaita novērojam gandrīz horizontālu līniju.

Visbeidzot, novērtēsim koda izpildes laika atšķirību:

Kods, lai noteiktu gradienta nolaišanās aprēķina laiku

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)

Varbūt mēs darām kaut ko nepareizi, bet atkal tā ir vienkārša “mājās rakstīta” funkcija, kas neizmanto bibliotēku nejutīgs pārspēj funkcijas aprēķina laiku, izmantojot bibliotēku nejutīgs.

Bet mēs nestāvam uz vietas, bet virzāmies uz vēl vienu aizraujošu veidu, kā atrisināt vienkāršu lineārās regresijas vienādojumu. Iepazīstieties!

Stohastiskā gradienta nolaišanās

Lai ātri saprastu stohastiskā gradienta nolaišanās darbības principu, labāk ir noteikt tā atšķirības no parastās gradienta nolaišanās. Mēs, gradienta nolaišanās gadījumā, atvasinājumu vienādojumos и izmantoja visu paraugā pieejamo pazīmju un patieso atbilžu vērtību summas (tas ir, visu и ). Stohastiskā gradienta nolaišanās gadījumā mēs neizmantosim visas izlasē esošās vērtības, bet tā vietā pseido nejauši atlasīsim tā saukto izlases indeksu un izmantosim tā vērtības.

Piemēram, ja indekss ir noteikts kā numurs 3 (trīs), tad mēs ņemam vērtības и , tad vērtības aizstājam atvasinātajos vienādojumos un nosakām jaunas koordinātas. Pēc tam, nosakot koordinātas, mēs atkal pseido nejauši nosakām izlases indeksu, aizstājam indeksam atbilstošās vērtības daļējos diferenciālvienādojumos un nosakām koordinātas jaunā veidā. и utt. līdz konverģence kļūst zaļa. No pirmā acu uzmetiena var šķist, ka tas nemaz nedarbosies, bet tas tā ir. Tiesa, ir vērts atzīmēt, ka kļūda nesamazinās ar katru soli, taču tendence noteikti ir.

Kādas ir stohastiskā gradienta nolaišanās priekšrocības salīdzinājumā ar parasto? Ja mūsu izlases lielums ir ļoti liels un mērāms desmitos tūkstošu vērtību, tad ir daudz vieglāk apstrādāt, teiksim, nejaušu tūkstoti no tiem, nevis visu izlasi. Šeit tiek izmantota stohastiskā gradienta nolaišanās. Mūsu gadījumā, protams, mēs nepamanīsim lielu atšķirību.

Apskatīsim kodu.

Stohastiskā gradienta nolaišanās kods

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

Mēs uzmanīgi skatāmies uz koeficientiem un pieķeram sevi uzdodam jautājumu "Kā tas var būt?" Mēs saņēmām citas koeficientu vērtības и . Varbūt stohastiskā gradienta nolaišanās ir atradusi optimālākus vienādojuma parametrus? Diemžēl nē. Pietiek aplūkot noviržu kvadrātu summu un redzēt, ka ar jaunām koeficientu vērtībām kļūda ir lielāka. Mēs nesteidzamies izmisumā. Izveidosim kļūdu izmaiņu grafiku.

Kods kvadrātveida noviržu summas attēlošanai stohastiskā gradienta nolaišanās gadījumā

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

Grafiks Nr. 5 “Noviržu summa kvadrātā stohastiskā gradienta nolaišanās laikā”

Skatoties grafiku, viss nostājas savās vietās un tagad visu salabosim.

Kas tad notika? Notika sekojošais. Ja mēs nejauši izvēlamies mēnesi, mūsu algoritms cenšas samazināt ieņēmumu aprēķināšanas kļūdu attiecībā uz atlasīto mēnesi. Pēc tam izvēlamies citu mēnesi un atkārtojam aprēķinu, bet otrajam izvēlētajam mēnesim kļūdu samazinām. Tagad atcerieties, ka pirmie divi mēneši ievērojami atšķiras no vienkāršā lineārās regresijas vienādojuma līnijas. Tas nozīmē, ka, izvēloties kādu no šiem diviem mēnešiem, samazinot katra no tiem kļūdu, mūsu algoritms nopietni palielina kļūdu visai izlasei. Tātad, ko darīt? Atbilde ir vienkārša: jums jāsamazina nolaišanās solis. Galu galā, samazinot nolaišanās pakāpi, kļūda arī pārtrauks “lēkt” uz augšu un uz leju. Pareizāk sakot, “lēkšanas” kļūda neapstāsies, bet ne tik ātri :) Pārbaudīsim.

Kods, lai palaistu SGD ar mazāku soli

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

Grafiks Nr. 6 “Noviržu kvadrātā summa stohastiskā gradienta nolaišanās laikā (80 tūkstoši soļu)”

Koeficienti ir uzlabojušies, taču joprojām nav ideāli. Hipotētiski to var labot šādā veidā. Mēs izvēlamies, piemēram, pēdējās 1000 iterācijās to koeficientu vērtības, ar kuriem tika pieļauta minimālā kļūda. Tiesa, šim nolūkam mums būs arī jāpieraksta pašu koeficientu vērtības. Mēs to nedarīsim, bet drīzāk pievērsīsim uzmanību grafikam. Tas izskatās gluds, un šķiet, ka kļūda vienmērīgi samazinās. Patiesībā tā nav taisnība. Apskatīsim pirmās 1000 iterācijas un salīdzināsim tās ar pēdējām.

SGD diagrammas kods (pirmie 1000 soļi)

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

Grafiks Nr. 7 “Noviržu kvadrātā summa SGD (pirmie 1000 soļi)”

8. grafiks “Noviržu kvadrātā summa SGD (pēdējie 1000 soļi)”

Pašā nolaišanās sākumā mēs novērojam diezgan vienmērīgu un strauju kļūdu samazināšanos. Pēdējās iterācijās mēs redzam, ka kļūda iet ap un ap vērtību 1,475 un dažos brīžos pat ir vienāda ar šo optimālo vērtību, bet tad tā joprojām palielinās... Es atkārtoju, jūs varat pierakstīt vērtības koeficienti и un pēc tam atlasiet tos, kuriem kļūda ir minimāla. Tomēr mums bija nopietnāka problēma: mums bija jāveic 80 tūkstoši soļu (skatiet kodu), lai iegūtu vērtības, kas tuvas optimālajām. Un tas jau ir pretrunā ar ideju ietaupīt aprēķina laiku, izmantojot stohastisku gradienta nolaišanos attiecībā pret gradienta nolaišanos. Ko var labot un uzlabot? Nav grūti pamanīt, ka pirmajās iterācijās mēs pārliecinoši ejam uz leju, tāpēc pirmajās iterācijās jāatstāj liels solis un, virzoties uz priekšu, solis jāsamazina. Šajā rakstā mēs to nedarīsim - tas jau ir pārāk garš. Tie, kas vēlas, var paši padomāt, kā to izdarīt, tas nav grūti :)

Tagad veiksim stohastisko gradienta nolaišanos, izmantojot bibliotēku nejutīgs (un nepaklupsim aiz akmeņiem, kurus mēs identificējām iepriekš)

Stohastiskā gradienta nolaišanās kods (NumPy)

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

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

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

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


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



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

Vērtības izrādījās gandrīz tādas pašas kā nolaižoties bez lietošanas nejutīgs. Tomēr tas ir loģiski.

Noskaidrosim, cik ilgi mums prasīja stohastiskā gradienta nolaišanās.

Kods SGD aprēķina laika noteikšanai (80 tūkstoši soļu)

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)

Jo tālāk mežā, jo tumšāki mākoņi: atkal “pašrakstītā” formula parāda labāko rezultātu. Tas viss liek domāt, ka ir jābūt vēl smalkākiem bibliotēkas izmantošanas veidiem nejutīgs, kas patiešām paātrina skaitļošanas darbības. Šajā rakstā mēs par tiem neuzzināsim. Būs par ko padomāt brīvajā laikā :)

Mēs apkopojam

Pirms apkopošanas es vēlētos atbildēt uz jautājumu, kas, visticamāk, radās no mūsu dārgā lasītāja. Kāpēc patiesībā ir tādas "spīdzināšanas" ar nolaišanos, kāpēc mums ir jāstaigā kalnā augšā un lejā (galvenokārt lejā), lai atrastu dārgo zemieni, ja mūsu rokās ir tik spēcīga un vienkārša ierīce, analītiskā risinājuma veidā, kas mūs uzreiz teleportē uz pareizo vietu?

Atbilde uz šo jautājumu slēpjas virspusē. Tagad mēs esam apskatījuši ļoti vienkāršu piemēru, kurā ir patiesā atbilde atkarīgs no vienas zīmes . Jūs to dzīvē neredzat bieži, tāpēc iedomāsimies, ka mums ir 2, 30, 50 vai vairāk zīmes. Katram atribūtam pievienosim tūkstošiem vai pat desmitiem tūkstošu vērtību. Šādā gadījumā analītiskais risinājums var neizturēt pārbaudi un neizdoties. Savukārt gradienta nolaišanās un tās variācijas lēnām, bet noteikti tuvinās mūs mērķim – funkcijas minimumam. Un neuztraucieties par ātrumu — mēs, iespējams, izskatīsim veidus, kas ļaus mums iestatīt un regulēt soļa garumu (tas ir, ātrumu).

Un tagad īstais īss kopsavilkums.

Pirmkārt, es ceru, ka rakstā sniegtais materiāls palīdzēs iesācējiem “datu zinātniekiem” saprast, kā atrisināt vienkāršus (un ne tikai) lineārās regresijas vienādojumus.

Otrkārt, mēs apskatījām vairākus vienādojuma risināšanas veidus. Tagad, atkarībā no situācijas, mēs varam izvēlēties to, kas ir vispiemērotākais problēmas risināšanai.

Treškārt, mēs redzējām papildu iestatījumu spēku, proti, gradienta nolaišanās soļa garumu. Šo parametru nevar ignorēt. Kā minēts iepriekš, lai samazinātu aprēķinu izmaksas, nolaišanās laikā jāmaina pakāpiena garums.

Ceturtkārt, mūsu gadījumā “mājās rakstītās” funkcijas uzrādīja vislabākos laika rezultātus aprēķiniem. Iespējams, tas ir saistīts ar ne visprofesionālāko bibliotēkas iespēju izmantošanu nejutīgs. Bet, lai kā arī būtu, šāds secinājums liek domāt par sevi. No vienas puses, reizēm ir vērts apšaubīt iedibinātos viedokļus, no otras puses, ne vienmēr ir vērts visu sarežģīt – gluži otrādi, reizēm vienkāršāks problēmas risināšanas veids ir efektīvāks. Un tā kā mūsu mērķis bija analizēt trīs vienkāršas lineārās regresijas vienādojuma risināšanas pieejas, mums pilnīgi pietika ar “pašu rakstītu” funkciju izmantošanu.