🥇 Yksinkertaisen lineaarisen regression yhtälön ratkaiseminen

Artikkelissa käsitellään useita tapoja määrittää yksinkertaisen (parillisen) regressiosuoran matemaattinen yhtälö.

Kaikki tässä käsitellyt yhtälön ratkaisumenetelmät perustuvat pienimmän neliösumman menetelmään. Merkitään menetelmät seuraavasti:

Analyyttinen ratkaisu
Gradientti laskeutuminen
Stokastinen gradienttilasku

Kullekin menetelmälle suoran yhtälön ratkaisemiseksi artikkeli tarjoaa erilaisia funktioita, jotka on jaettu pääasiassa niihin, jotka on kirjoitettu ilman kirjastoa nuhjuinen ja ne, joita käytetään laskelmissa nuhjuinen. Uskotaan, että taitava käyttö nuhjuinen alentaa laskentakustannuksia.

Kaikki artikkelissa annettu koodi on kirjoitettu kielellä python-2.7 kanssa Jupyter Notebook. Lähdekoodi ja tiedosto esimerkkitiedoilla on lähetetty Github

Artikkeli on suunnattu enemmän sekä aloittelijoille että niille, jotka ovat jo vähitellen alkaneet hallita erittäin laajaa tekoälyn osaa - koneoppimista.

Materiaalin havainnollistamiseksi käytämme hyvin yksinkertaista esimerkkiä.

Esimerkkiehdot

Meillä on viisi arvoa, jotka kuvaavat riippuvuutta Y alkaen X (taulukko nro 1):

Taulukko 1 "Esimerkkiehdot"

Oletetaan, että arvot on vuoden kuukausi ja - tulot tässä kuussa. Toisin sanoen tulot riippuvat vuoden kuukaudesta ja - ainoa merkki, josta tulot riippuvat.

Esimerkki on niin ja niin, sekä tulojen ehdollisen riippuvuuden kannalta vuoden kuukaudesta että arvojen lukumäärän kannalta - niitä on hyvin vähän. Tällainen yksinkertaistaminen mahdollistaa kuitenkin, kuten he sanovat, selittää, ei aina helposti, materiaalia, jonka aloittelijat omaksuvat. Ja myös numeroiden yksinkertaisuus antaa niille, jotka haluavat ratkaista esimerkin paperilla ilman merkittäviä työvoimakustannuksia.

Oletetaan, että esimerkissä annettu riippuvuus voidaan approksimoida melko hyvin muodon yksinkertaisen (parillisen) regressiosuoran matemaattisella yhtälöllä:

missä on kuukausi, jolloin tulot on saatu, — kuukautta vastaava tulo, и ovat estimoidun suoran regressiokertoimia.

Huomaa, että kerroin kutsutaan usein arvioidun viivan kaltevuudeksi tai kaltevuudeksi; edustaa määrää, jolla kun se muuttuu .

Ilmeisesti tehtävämme esimerkissä on valita tällaiset kertoimet yhtälöstä и , jossa laskettujen tuloarvojemme poikkeamat kuukausittain oikeista vastauksista, ts. näytteessä esitetyt arvot ovat minimaalisia.

Pienimmän neliön menetelmä

Pienimmän neliösumman menetelmän mukaan poikkeama tulee laskea neliöimällä se. Tämän tekniikan avulla voit välttää poikkeamien keskinäisen peruuntumisen, jos niillä on päinvastaisia merkkejä. Esimerkiksi jos yhdessä tapauksessa, poikkeama on +5 (plus viisi) ja toisessa -5 (miinus viisi), silloin poikkeamien summa kumoaa toisensa ja on 0 (nolla). Sinun ei tarvitse neliöida poikkeamaa, vaan käytä moduuliominaisuutta ja sitten kaikki poikkeamat ovat positiivisia ja kerääntyvät. Emme käsittele tätä kohtaa yksityiskohtaisesti, vaan osoitamme yksinkertaisesti, että laskelmien mukavuuden vuoksi on tapana neliöttää poikkeama.

Tältä näyttää kaava, jolla määritämme pienimmän neliösumman poikkeamia (virheitä):

missä on todellisten vastausten (eli laskemiemme tulojen) likiarvon funktio,

ovatko vastaukset oikeita (tulot otoksessa),

on näyteindeksi (sen kuukauden numero, jolloin poikkeama määritetään)

Erotetaan funktio, määritellään osittaisdifferentiaaliyhtälöt ja ollaan valmiita siirtymään analyyttiseen ratkaisuun. Mutta ensin tehdään lyhyt retki siitä, mitä differentiaatio on, ja muistetaan derivaatan geometrinen merkitys.

Erilaistuminen

Differentiointi on operaatio, jolla löydetään funktion derivaatta.

Mihin johdannaista käytetään? Funktion derivaatta luonnehtii funktion muutosnopeutta ja kertoo sen suunnan. Jos derivaatta tietyssä pisteessä on positiivinen, funktio kasvaa, muuten funktio pienenee. Ja mitä suurempi on absoluuttisen derivaatan arvo, sitä suurempi on funktion arvojen muutosnopeus ja sitä jyrkempi funktiokaavion kaltevuus.

Esimerkiksi karteesisen koordinaattijärjestelmän olosuhteissa derivaatan arvo pisteessä M(0,0) on yhtä suuri kuin + 25 tarkoittaa, että tietyssä pisteessä, kun arvoa siirretään oikealle sopimuksella, arvolla kasvaa 25 tavanomaisella yksiköllä. Kaaviossa se näyttää melko jyrkältä arvojen nousulta tietystä pisteestä.

Toinen esimerkki. Johdannainen arvo on sama -0,1 tarkoittaa sitä siirrettynä yhtä tavanomaista yksikköä kohti, arvo pienenee vain 0,1 tavanomaisella yksiköllä. Samanaikaisesti funktion kaaviossa voimme havaita tuskin havaittavan alaspäin. Vuoreen verrattaessa, tuntuu kuin laskeutuisimme hitaasti loivaa rinnettä vuorelta, toisin kuin edellisessä esimerkissä, jossa meidän piti kiivetä erittäin jyrkkiä huippuja :)

Näin ollen funktion erottamisen jälkeen kertoimilla и , määrittelemme ensimmäisen asteen osittaisdifferentiaaliyhtälöt. Yhtälöiden määrittämisen jälkeen saamme kahden yhtälön järjestelmän, jonka ratkaisemalla voimme valita tällaiset kertoimien arvot и , joiden vastaavien johdannaisten arvot muuttuvat tietyissä kohdissa hyvin, hyvin vähän, eivätkä analyyttisen ratkaisun tapauksessa muutu ollenkaan. Toisin sanoen virhefunktio löydetyissä kertoimissa saavuttaa minimin, koska osittaisten derivaattojen arvot näissä kohdissa ovat nolla.

Eli differentiaatiosääntöjen mukaan ensimmäisen asteen osittaisderivaatayhtälö suhteessa kertoimeen tulee muodossa:

1. asteen osittaisen derivaatan yhtälö suhteessa tulee muodossa:

Tuloksena saimme yhtälöjärjestelmän, jolla on melko yksinkertainen analyyttinen ratkaisu:

aloita{yhtälö*}
aloita{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
lopeta{tapaukset}
loppu{yhtälö*}

Ennen yhtälön ratkaisemista esilatataan, tarkistetaan, että lataus on oikein, ja muotoillaan tiedot.

Tietojen lataus ja muotoilu

On huomattava, että koska analyyttisen ratkaisun ja sen jälkeen gradientin ja stokastisen gradientin laskeutumisen yhteydessä käytämme koodia kahdessa muunnelmassa: kirjaston avulla. nuhjuinen ja ilman sitä, tarvitsemme asianmukaisen datan muotoilun (katso koodi).

Tietojen lataus- ja käsittelykoodi

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

Visualisointi

Nyt, kun olemme ensinnäkin ladaneet tiedot, toiseksi tarkastaneet latauksen oikeellisuuden ja lopuksi alustaneet tiedot, suoritamme ensimmäisen visualisoinnin. Usein tähän käytetty menetelmä on paripiirros kirjasto seaborn. Esimerkissämme kirjaston käyttämisessä ei ole mitään järkeä rajoitetun määrän vuoksi seaborn. Käytämme tavallista kirjastoa Matplotlib ja katso vain sirontakaaviota.

Scatterplot-koodi

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

Kaavio nro 1 "Tulojen riippuvuus vuoden kuukaudesta"

Analyyttinen ratkaisu

Käytetään yleisimpiä työkaluja pytonkäärme ja ratkaise yhtälöjärjestelmä:

aloita{yhtälö*}
aloita{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
lopeta{tapaukset}
loppu{yhtälö*}

Cramerin säännön mukaan löydämme yleisdeterminantin sekä determinantit by ja , jonka jälkeen determinantti jaetaan arvolla yleisdeterminantille - löydä kerroin , vastaavasti löydämme kertoimen .

Analyyttisen ratkaisun koodi

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

Tässä on mitä saimme:

Joten kertoimien arvot on löydetty, neliöityjen poikkeamien summa on määritetty. Piirretään sirontahistogrammiin suora viiva löydettyjen kertoimien mukaisesti.

Regressioviivakoodi

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

Kaavio nro 2 "Oikeita ja laskettuja vastauksia"

Voit tarkastella kunkin kuukauden poikkeamakaaviota. Meidän tapauksessamme emme saa siitä mitään merkittävää käytännön arvoa, mutta tyydytämme uteliaisuutemme siitä, kuinka hyvin yksinkertainen lineaarinen regressioyhtälö luonnehtii tulojen riippuvuutta vuoden kuukaudesta.

Poikkeamakaavion koodi

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

Kaavio nro 3 "Poikkeamat, %"

Ei täydellinen, mutta suoritimme tehtävämme.

Kirjoitetaan funktio, joka määrittää kertoimet и käyttää kirjastoa nuhjuinen, tarkemmin sanottuna, kirjoitamme kaksi funktiota: toinen pseudoinversomatriisilla (ei suositella käytännössä, koska prosessi on laskennallisesti monimutkainen ja epävakaa), toinen matriisiyhtälöllä.

Analyyttisen ratkaisun koodi (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

Verrataan kertoimien määrittämiseen käytettyä aikaa и 3 esitetyn menetelmän mukaisesti.

Koodi laskenta-ajan laskemiseen

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)

Pienellä datamäärällä eteen tulee "itsekirjoitettu" funktio, joka etsii kertoimet Cramerin menetelmällä.

Nyt voit siirtyä muihin tapoihin löytää kertoimet и .

Gradientti laskeutuminen

Ensin määritellään, mikä gradientti on. Yksinkertaisesti sanottuna gradientti on segmentti, joka osoittaa funktion suurimman kasvun suunnan. Vastaavasti vuorelle kiipeämisen kanssa, missä kaltevuuspinnat ovat jyrkimmällä nousulla vuoren huipulle. Kehittäen esimerkkiä vuoren kanssa muistamme, että itse asiassa tarvitsemme jyrkimmän laskun päästäksemme alangolle mahdollisimman nopeasti, eli minimiin - paikkaan, jossa toiminto ei kasva tai vähene. Tässä vaiheessa derivaatta on yhtä suuri kuin nolla. Siksi emme tarvitse gradienttia, vaan antigradienttia. Antigradientin löytämiseksi sinun tarvitsee vain kertoa gradientti -1 (miinus yksi).

Kiinnittäkäämme huomiota siihen, että funktiolla voi olla useita minimiä, ja kun olemme laskeneet yhteen niistä alla ehdotetulla algoritmilla, emme löydä toista minimiä, joka voi olla pienempi kuin löydetty. Rentoudutaan, tämä ei ole meille uhka! Meidän tapauksessamme on kyse yhdestä vähimmäismäärästä, koska funktiomme kaaviossa on säännöllinen paraabeli. Ja kuten meidän kaikkien pitäisi tietää hyvin koulumatematiikan kurssistamme, paraabelilla on vain yksi minimi.

Kun saimme selville, miksi tarvitsemme gradientin, ja myös, että gradientti on segmentti, eli vektori, jolla on tietyt koordinaatit, jotka ovat täsmälleen samat kertoimet и voimme toteuttaa kaltevuuden laskeutumisen.

Ennen kuin aloitat, suosittelen lukemaan vain muutaman lauseen laskeutumisalgoritmista:

Määritämme pseudosatunnaisella tavalla kertoimien koordinaatit и . Esimerkissämme määritetään kertoimet lähellä nollaa. Tämä on yleinen käytäntö, mutta jokaisella tapauksella voi olla oma käytäntönsä.
Koordinaatista vähennä ensimmäisen asteen osittaisen derivaatan arvo pisteessä . Joten jos derivaatta on positiivinen, funktio kasvaa. Siksi johdannaisen arvon vähentämisellä siirrymme kasvun vastakkaiseen suuntaan eli laskusuuntaan. Jos derivaatta on negatiivinen, niin funktio tässä pisteessä pienenee ja derivaatan arvoa vähentämällä siirrytään laskeutumissuuntaan.
Suoritamme samanlaisen toimenpiteen koordinaatilla : vähennä osittaisen derivaatan arvo pisteessä .
Jotta et hypätä minimin yli ja lennä syvään avaruuteen, on tarpeen asettaa askelkoko laskeutumissuuntaan. Yleisesti ottaen voit kirjoittaa kokonaisen artikkelin siitä, kuinka askel asetetaan oikein ja miten sitä muutetaan laskeutumisprosessin aikana laskennallisten kustannusten vähentämiseksi. Mutta nyt meillä on hieman erilainen tehtävä, ja määritämme askelkoon tieteellisellä "poke" -menetelmällä tai, kuten yleisessä kielessä sanotaan, empiirisesti.
Kun olemme annetuista koordinaateista и vähennä derivaattojen arvot, saamme uudet koordinaatit и . Otamme seuraavan vaiheen (vähennyslasku), jo lasketuista koordinaateista. Ja niin sykli alkaa uudestaan ja uudestaan, kunnes vaadittu konvergenssi saavutetaan.

Kaikki! Nyt olemme valmiita lähtemään etsimään Mariana-haudon syvintä rotkoa. Aloitetaan.

Koodi kaltevuuden laskulle

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

Sukelsimme Mariana-haudon pohjalle ja löysimme sieltä kaikki samat kerroinarvot и , mikä oli juuri sitä mitä oli odotettavissa.

Otetaan toinen sukellus, vain tällä kertaa syvänmeren ajoneuvomme on täynnä muita teknologioita, nimittäin kirjastoa nuhjuinen.

Kaltevuuden laskukoodi (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

Kertoimien arvot и muuttumaton.

Katsotaan kuinka virhe muuttui gradientin laskeutumisen aikana, eli kuinka neliöityjen poikkeamien summa muuttui jokaisen askeleen aikana.

Koodi neliöpoikkeamien summien piirtämiseen

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

Kaavio nro 4 "Poikkeamien neliösumma kaltevuuden laskun aikana"

Kaaviossa näemme, että jokaisella askeleella virhe pienenee ja tietyn iteraatiomäärän jälkeen havaitsemme melkein vaakasuuntaisen viivan.

Lopuksi arvioidaan ero koodin suoritusajassa:

Koodi, jolla määritetään gradientin laskeutumisaika

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)

Ehkä teemme jotain väärin, mutta taas se on yksinkertainen "kotiin kirjoitettu" toiminto, joka ei käytä kirjastoa nuhjuinen ylittää kirjastoa käyttävän funktion laskenta-ajan nuhjuinen.

Mutta emme seiso paikallaan, vaan siirrymme tutkimaan toista jännittävää tapaa ratkaista yksinkertainen lineaarinen regressioyhtälö. Tavata!

Stokastinen gradienttilasku

Stokastisen gradienttilaskeutumisen toimintaperiaatteen ymmärtämiseksi nopeasti on parempi määrittää sen erot tavallisesta gradienttilaskeutumisesta. Me, gradientin laskeutumisen tapauksessa, derivaatan yhtälöissä и käytti kaikkien otoksessa saatavilla olevien ominaisuuksien ja todellisten vastausten arvojen summia (eli kaikkien и ). Stokastisessa gradienttilaskeutumisessa emme käytä kaikkia näytteessä olevia arvoja, vaan valitsemme näennäissatunnaisesti ns. näyteindeksin ja käytämme sen arvoja.

Jos indeksiksi määritetään esimerkiksi numero 3 (kolme), otamme arvot и , sitten korvaamme arvot johdannaisyhtälöihin ja määritämme uudet koordinaatit. Sitten, kun koordinaatit on määritetty, määritämme taas näennäissatunnaisesti näyteindeksin, korvaamme indeksiä vastaavat arvot osittaisdifferentiaaliyhtälöihin ja määritämme koordinaatit uudella tavalla и jne. kunnes konvergenssi muuttuu vihreäksi. Ensi silmäyksellä saattaa tuntua, että tämä ei toimi ollenkaan, mutta toimii. On totta, että virhe ei pienene joka askeleella, mutta suuntaus on varmasti olemassa.

Mitä etuja stokastisesta gradienttilaskeutumisesta on tavanomaiseen verrattuna? Jos otoskokomme on erittäin suuri ja mitataan kymmenissä tuhansissa arvoissa, on paljon helpompi käsitellä vaikkapa satunnainen tuhat niistä kuin koko otos. Tässä tulee esiin stokastinen gradienttilasku. Meidän tapauksessamme emme tietenkään huomaa suurta eroa.

Katsotaanpa koodia.

Koodi stokastista gradienttilaskua varten

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

Tarkastelemme kertoimia huolellisesti ja huomaamme kysyvämme kysymyksen "Miten tämä voi olla?" Saimme muita kertoimia и . Ehkä stokastinen gradienttilasku on löytänyt optimaalisempia parametreja yhtälölle? Valitettavasti ei. Riittää, kun tarkastellaan neliöityjen poikkeamien summaa ja nähdään, että kertoimien uusilla arvoilla virhe on suurempi. Meillä ei ole kiirettä epätoivoon. Rakennetaan kaavio virheen muutoksesta.

Koodi stokastisen gradientin laskeutumisen neliöityjen poikkeamien summan piirtämiseen

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

Kaavio nro 5 "Stokastisen gradientin laskeutumisen neliöpoikkeamien summa"

Aikataulua katsoessa kaikki loksahtaa paikoilleen ja nyt korjaamme kaiken.

Mitä tapahtui? Tapahtui seuraavaa. Kun valitsemme satunnaisesti kuukauden, algoritmimme pyrkii vähentämään virhettä tulojen laskennassa valitun kuukauden osalta. Sitten valitsemme toisen kuukauden ja toistamme laskennan, mutta vähennämme virhettä toiselle valitulle kuukaudelle. Muista nyt, että kaksi ensimmäistä kuukautta poikkeavat merkittävästi yksinkertaisen lineaarisen regressioyhtälön viivasta. Tämä tarkoittaa, että kun mikä tahansa näistä kahdesta kuukaudesta valitaan, pienentämällä kummankin virhettä algoritmimme lisää vakavasti koko otoksen virhettä. Eli mikä neuvoksi? Vastaus on yksinkertainen: sinun on vähennettävä laskeutumisaskelta. Loppujen lopuksi, vähentämällä laskuaskelta, virhe myös lopettaa "hyppäämisen" ylös ja alas. Tai pikemminkin "hyppy" virhe ei lopu, mutta se ei tee sitä niin nopeasti :) Tarkistetaan.

Koodi ajaa SGD pienemmillä askelilla

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

Kaavio nro 6 "Stokastisen gradientin laskeutumisen aikana (80 tuhatta askelta) poikkeamien neliöityjen summa"

Kertoimet ovat parantuneet, mutta eivät vieläkään ihanteellisia. Hypoteettisesti tämä voidaan korjata tällä tavalla. Valitsemme esimerkiksi viimeisen 1000 iteraatiossa kertoimien arvot, joilla pienin virhe tehtiin. Totta, tätä varten meidän on myös kirjoitettava itse kertoimien arvot. Emme tee tätä, vaan kiinnitämme huomiota aikatauluun. Se näyttää sileältä ja virhe näyttää vähenevän tasaisesti. Itse asiassa tämä ei ole totta. Katsotaanpa ensimmäisiä 1000 iteraatiota ja verrataan niitä viimeisiin.

SGD-kaavion koodi (ensimmäiset 1000 askelta)

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

Kaavio nro 7 "Poikkeamien neliösumma SGD (ensimmäiset 1000 askelta)"

Kaavio nro 8 "Poikkeamien neliösumma SGD (viimeiset 1000 askelta)"

Aivan laskeutumisen alussa havaitsemme melko tasaisen ja jyrkän virheen pienenemisen. Viimeisissä iteraatioissa näemme, että virhe kiertää arvon 1,475 ja välillä on jopa yhtä suuri kuin tämä optimaalinen arvo, mutta sitten se silti nousee... Toistan, voit kirjoittaa muistiin kertoimet и ja valitse sitten ne, joissa virhe on minimaalinen. Meillä oli kuitenkin vakavampi ongelma: meidän piti ottaa 80 tuhatta askelta (katso koodi) saadaksemme arvot lähellä optimaalista. Ja tämä on jo ristiriidassa ajatuksen kanssa säästää laskenta-aikaa stokastisella gradienttilaskeutumisella suhteessa gradienttilaskeutumiseen. Mitä voidaan korjata ja parantaa? Ei ole vaikeaa huomata, että ensimmäisissä iteraatioissa menemme luottavaisesti alaspäin ja siksi meidän tulee jättää suuri askel ensimmäisissä iteraatioissa ja pienentää askelta eteenpäin siirtyessämme. Emme tee tätä tässä artikkelissa - se on jo liian pitkä. Halukkaat voivat itse miettiä, miten tämä tehdään, se ei ole vaikeaa :)

Suoritetaan nyt stokastinen gradienttilasku kirjaston avulla nuhjuinen (äläkä kompastele aiemmin tunnistamiimme kiviin)

Stokastisen gradientin laskeutumisen koodi (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

Arvot osoittautuivat lähes samoiksi kuin ilman käyttöä laskeutuessa nuhjuinen. Tämä on kuitenkin loogista.

Selvitetään kuinka kauan stokastiset gradienttilaskeutumiset veivät.

Koodi SGD-laskentaajan määrittämiseksi (80 tuhatta askelta)

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)

Mitä pidemmälle metsään, sitä tummemmat pilvet ovat: jälleen ”itsekirjoitettu” kaava näyttää parhaan tuloksen. Kaikki tämä viittaa siihen, että kirjaston käyttöön täytyy olla vielä hienovaraisempia tapoja nuhjuinen, jotka todella nopeuttavat laskentatoimia. Tässä artikkelissa emme opi niistä. Vapaa-ajalla on ajateltavaa :)

Me tiivistämme

Ennen kuin teen yhteenvedon, haluaisin vastata kysymykseen, joka todennäköisimmin nousi rakkaalle lukijallemme. Miksi itse asiassa tällainen "kidutus" laskeutumisilla, miksi meidän täytyy kävellä ylös ja alas vuorelta (enimmäkseen alas) löytääksemme arvokkaan alangon, jos meillä on käsissämme niin tehokas ja yksinkertainen laite, analyyttisen ratkaisun muodossa, joka teleporttaa meidät välittömästi oikeaan paikkaan?

Vastaus tähän kysymykseen on pinnalla. Nyt olemme tarkastelleet hyvin yksinkertaista esimerkkiä, jossa oikea vastaus on riippuu yhdestä merkistä . Et näe tätä usein elämässä, joten kuvitellaan, että meillä on 2, 30, 50 tai enemmän merkkiä. Lisätään tähän tuhansia tai jopa kymmeniä tuhansia arvoja jokaiselle attribuutille. Tässä tapauksessa analyyttinen liuos ei ehkä kestä testiä ja epäonnistuu. Gradienttilasku ja sen vaihtelut puolestaan vievät meidät hitaasti mutta varmasti lähemmäs tavoitetta - funktion minimiä. Ja älä ole huolissasi nopeudesta - tutkimme todennäköisesti tapoja, joiden avulla voimme asettaa ja säätää askelpituutta (eli nopeutta).

Ja nyt varsinainen lyhyt yhteenveto.

Ensinnäkin toivon, että artikkelissa esitetty materiaali auttaa aloittelevia "tietotieteilijöitä" ymmärtämään, kuinka ratkaista yksinkertaisia (eikä vain) lineaarisia regressioyhtälöitä.

Toiseksi tarkastelimme useita tapoja ratkaista yhtälö. Nyt voimme tilanteesta riippuen valita sen, joka sopii parhaiten ongelman ratkaisemiseen.

Kolmanneksi näimme lisäasetusten tehon, nimittäin kaltevuuden laskuaskeleen pituuden. Tätä parametria ei voi jättää huomiotta. Kuten edellä todettiin, laskelmien kustannusten vähentämiseksi askelpituutta tulisi muuttaa laskeutumisen aikana.

Neljänneksi, meidän tapauksessamme "kotikirjoitetut" funktiot osoittivat parhaat aikatulokset laskelmiin. Tämä johtuu luultavasti kirjaston ominaisuuksien ei ammattimaisimmasta käytöstä nuhjuinen. Mutta olipa tilanne miten tahansa, seuraava johtopäätös viittaa itsestään. Toisaalta joskus kannattaa kyseenalaistaa vakiintuneita mielipiteitä, toisaalta kaikkea ei aina kannata monimutkaista - päinvastoin, joskus yksinkertaisempi tapa ratkaista ongelma on tehokkaampi. Ja koska tavoitteemme oli analysoida kolmea lähestymistapaa yksinkertaisen lineaarisen regressioyhtälön ratkaisemiseen, "itsekirjoitettujen" funktioiden käyttö riitti meille.