🥇Az egyszerű lineáris regresszió egyenletének megoldása

A cikk számos módszert tárgyal egy egyszerű (páros) regressziós egyenes matematikai egyenletének meghatározására.

Az itt tárgyalt egyenlet megoldásának valamennyi módszere a legkisebb négyzetek módszerén alapul. Jelöljük a módszereket a következőképpen:

Analitikai megoldás
Gradiens Descent
Sztochasztikus gradiens süllyedés

Az egyenes egyenlet megoldásának mindegyik módszeréhez a cikk különféle függvényeket kínál, amelyek főként a könyvtár használata nélkül íródnak. numpy és azokat, amelyeket számításokhoz használnak numpy. Úgy tartják, hogy ügyes használat numpy csökkenti a számítási költségeket.

A cikkben szereplő összes kód a nyelven van írva Python 2.7 felhasználásával Jupyter Jegyzetfüzet. A forráskódot és a mintaadatokat tartalmazó fájlt a következő helyen teszik közzé Github

A cikk inkább a kezdőknek és azoknak szól, akik már fokozatosan elkezdték elsajátítani a mesterséges intelligencia egy nagyon széles szakaszát - a gépi tanulást.

Az anyag illusztrálására egy nagyon egyszerű példát használunk.

Példa feltételek

Öt értékünk van, amelyek a függőséget jellemzik Y -tól X (1. számú táblázat):

1. számú táblázat „Példafeltételek”

Feltételezzük, hogy az értékek az év hónapja, és — bevétel ebben a hónapban. Más szóval, a bevétel az év hónapjától függ, és - az egyetlen jel, amelytől a bevétel függ.

A példa ilyen-olyan, mind a bevétel feltételes függése az év hónapjától, mind az értékek száma szempontjából - nagyon kevés van belőlük. Egy ilyen egyszerűsítés azonban lehetővé teszi, hogy – amint mondják – a kezdők által megszerzett anyagot nem mindig könnyen elmagyarázzuk. És a számok egyszerűsége lehetővé teszi azoknak, akik jelentős munkaerőköltségek nélkül akarják megoldani a példát papíron.

Tételezzük fel, hogy a példában megadott függés jól közelíthető egy egyszerű (páros) regressziós egyenes matematikai egyenletével:

ahol az a hónap, amelyben a bevétel megérkezett, — a hónapnak megfelelő bevétel, и a becsült egyenes regressziós együtthatói.

Vegye figyelembe, hogy az együttható gyakran a becsült vonal meredekségének vagy gradiensének nevezik; azt az összeget jelenti, amellyel a amikor megváltozik .

Nyilvánvalóan az a feladatunk a példában, hogy ilyen együtthatókat válasszunk az egyenletben и , amelynél a számított bevételi értékeink havi eltérései az igaz válaszoktól, pl. a mintában szereplő értékek minimálisak lesznek.

Legkisebb négyzet alakú módszer

A legkisebb négyzetek módszere szerint az eltérést négyzetre emelve kell kiszámítani. Ez a technika lehetővé teszi az eltérések kölcsönös törlésének elkerülését, ha ellentétes előjelűek. Például, ha egy esetben az eltérés az +5 (plusz öt), és a másikban -5 (mínusz öt), akkor az eltérések összege kioltja egymást, és 0 (nulla) lesz. Nem kell négyzetre emelni az eltérést, hanem használjuk a modulus tulajdonságot, és akkor az összes eltérés pozitív lesz és felhalmozódik. Ezen a ponton nem térünk ki részletesen, hanem egyszerűen jelezzük, hogy a számítások megkönnyítése érdekében az eltérést négyzetre szokás emelni.

Így néz ki a képlet, amellyel meghatározzuk az eltérések (hibák) legkisebb négyzetes összegét:

ahol az igaz válaszok (vagyis az általunk kiszámított bevétel) közelítésének függvénye,

igazak-e a válaszok (a mintában szereplő bevétel),

a minta indexe (annak a hónapnak a száma, amelyben az eltérést meghatározták)

Differenciáljuk a függvényt, definiáljuk a parciális differenciálegyenleteket, és készen állunk az analitikus megoldásra. Először azonban vessünk egy rövid kirándulást arról, hogy mi a differenciálás, és emlékezzünk a származék geometriai jelentésére.

Különbségtétel

A differenciálás egy függvény deriváltjának megtalálásának művelete.

Mire használják a származékot? Egy függvény deriváltja jellemzi a függvény változási sebességét és megmondja annak irányát. Ha a derivált egy adott pontban pozitív, akkor a függvény növekszik, ellenkező esetben a függvény csökken. És minél nagyobb az abszolút derivált értéke, annál nagyobb a függvényértékek változási sebessége, valamint annál meredekebb a függvénygráf meredeksége.

Például egy derékszögű koordinátarendszer feltételei között a derivált értéke az M(0,0) pontban egyenlő +25 azt jelenti, hogy egy adott ponton, amikor az érték eltolódik jobbra egyezményes mértékegységgel, értékkel 25 hagyományos egységgel nő. A grafikonon az értékek meglehetősen meredek emelkedése látszik adott pontból.

Egy másik példa. A derivált értéke egyenlő -0,1 azt jelenti, hogy amikor eltoljuk egy egyezményes egységre, értékre csak 0,1 hagyományos egységgel csökken. Ugyanakkor a függvény grafikonján alig észrevehető lejtőt figyelhetünk meg. A hegyekkel való analógiát levonva olyan, mintha nagyon lassan egy enyhe lejtőn ereszkednénk le a hegyről, ellentétben az előző példával, ahol nagyon meredek csúcsokat kellett megmásznunk :)

Így a függvény megkülönböztetése után esélyek szerint и , 1. rendű parciális differenciálegyenleteket definiálunk. Az egyenletek meghatározása után egy két egyenletrendszert kapunk, melynek megoldásával kiválaszthatjuk az együtthatók ilyen értékeit и , amelynél a megfelelő deriváltak értékei adott pontokban nagyon-nagyon csekély mértékben változnak, analitikus megoldás esetén pedig egyáltalán nem. Más szavakkal, a talált együtthatók hibafüggvénye eléri a minimumot, mivel a parciális deriváltak értéke ezeken a pontokon nulla lesz.

Tehát a differenciálás szabályai szerint az együtthatóhoz viszonyított 1. rendű parciális derivált egyenlet a következő formában lesz:

1. rendű parciális derivált egyenlet tekintetében a következő formában lesz:

Ennek eredményeként egy egyenletrendszert kaptunk, amelynek meglehetősen egyszerű analitikai megoldása van:

kezdődik{egyenlet*}
kezd{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
vége{esetek}
vége{egyenlet*}

Az egyenlet megoldása előtt töltsünk elő, ellenőrizzük a betöltés helyességét, és formázzuk meg az adatokat.

Adatok betöltése és formázása

Megjegyzendő, hogy az analitikai megoldáshoz, majd a gradiens és a sztochasztikus gradiens süllyedéshez a kódot két változatban fogjuk használni: a könyvtár használatával. numpy és használata nélkül, akkor megfelelő adatformázásra lesz szükségünk (lásd a kódot).

Adatbetöltési és feldolgozási kód

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

Megjelenítés

Most, miután először feltöltöttük az adatokat, másodsorban ellenőriztük a betöltés helyességét és végül formáztuk az adatokat, elvégezzük az első vizualizációt. Az erre gyakran használt módszer az párparcella könyvtár Tengeren született. Példánkban a korlátozott számok miatt nincs értelme a könyvtár használatának Tengeren született. A szokásos könyvtárat fogjuk használni matplotlib és nézd csak a szóródást.

Szórványkód

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

1. számú ábra „A bevétel függése az év hónapjától”

Analitikai megoldás

Használjuk a leggyakoribb eszközöket piton és oldja meg az egyenletrendszert:

kezdődik{egyenlet*}
kezd{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
vége{esetek}
vége{egyenlet*}

Cramer szabálya szerint megtaláljuk az általános determinánst, valamint a determinánsokat az által és , ami után a determinánst elosztva ezzel az általános determinánshoz - keresse meg az együtthatót , hasonlóképpen megtaláljuk az együtthatót .

Analitikai megoldás kódja

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

Íme, amit kaptunk:

Tehát az együtthatók értékeit megtaláltuk, az eltérések négyzetes összegét megállapítottuk. A szórási hisztogramon húzzunk egy egyenest a talált együtthatóknak megfelelően.

Regressziós vonal kódja

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

2. számú táblázat „Helyes és kiszámított válaszok”

Megtekintheti az egyes hónapok eltérési grafikonját. Esetünkben jelentős gyakorlati értéket nem merítünk belőle, de kielégítjük a kíváncsiságunkat, hogy az egyszerű lineáris regressziós egyenlet mennyire jellemzi a bevétel év hónaptól való függését.

Eltérési diagram kódja

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

3. számú diagram „Eltérések, %”

Nem tökéletes, de teljesítettük a feladatunkat.

Írjunk egy függvényt, amely az együtthatók meghatározásához и használja a könyvtárat numpy, pontosabban két függvényt fogunk felírni: az egyiket pszeudoinverz mátrix segítségével (a gyakorlatban nem ajánlott, mivel a folyamat számításilag bonyolult és instabil), a másikat mátrixegyenlet segítségével.

Analitikai megoldás kódja (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

Hasonlítsuk össze az együtthatók meghatározására fordított időt и , a bemutatott 3 módszer szerint.

Kód a számítási idő kiszámításához

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)

Kis mennyiségű adattal előjön egy „önírt” függvény, amely Cramer módszerével keresi meg az együtthatókat.

Most áttérhet más módszerekre az együtthatók megtalálására и .

Gradiens Descent

Először is határozzuk meg, mi az a gradiens. Egyszerűen fogalmazva, a gradiens egy szegmens, amely egy függvény maximális növekedési irányát jelzi. A hegymászás analógiájára, ahol a lejtős oldalak a legmeredekebb mászás a hegy tetejére. A hegyi példát továbbfejlesztve emlékezünk arra, hogy valójában a legmeredekebb ereszkedésre van szükség ahhoz, hogy a lehető leggyorsabban elérjük az alföldet, vagyis a minimumot - azt a helyet, ahol a funkció nem nő vagy csökken. Ezen a ponton a derivált egyenlő lesz nullával. Ezért nem gradiensre van szükségünk, hanem antigradiensre. Az antigradiens megtalálásához csak meg kell szoroznia a gradienst ezzel -1 (mínusz egy).

Figyeljünk arra, hogy egy függvénynek több minimuma is lehet, és az alábbiakban javasolt algoritmus segítségével ezek közül az egyikbe leereszkedve nem fogunk találni másik minimumot, amely alacsonyabb lehet a találtnál. Lazítsunk, ez nem fenyeget minket! Esetünkben egyetlen minimummal van dolgunk, hiszen a funkciónk a grafikonon egy szabályos parabola. És amint azt mindannyian nagyon jól tudnunk kell az iskolai matematika kurzusunkból, a parabolának csak egy minimuma van.

Miután megtudtuk, miért van szükségünk gradiensre, és azt is, hogy a gradiens egy szakasz, vagyis egy vektor adott koordinátákkal, amelyek pontosan ugyanazok az együtthatók и gradiens süllyedést tudunk megvalósítani.

Mielőtt elkezdené, javaslom, hogy olvasson el néhány mondatot az ereszkedési algoritmusról:

Álvéletlen módon határozzuk meg az együtthatók koordinátáit и . Példánkban nullához közeli együtthatókat fogunk meghatározni. Ez egy általános gyakorlat, de minden esetnek megvan a saját gyakorlata.
Koordináta alapján vonjuk le az I. rendű parciális derivált értékét a pontban . Tehát, ha a derivált pozitív, akkor a függvény növekszik. Ezért a derivált értékének levonásával a növekedéssel ellentétes, azaz a süllyedés irányába fogunk elmozdulni. Ha a derivált negatív, akkor a függvény ezen a ponton csökken, és a derivált értékét kivonva a süllyedés irányába haladunk.
Hasonló műveletet végzünk a koordinátával : vonja ki a részleges derivált értékét a pontban .
Annak érdekében, hogy ne ugorja át a minimumot és ne repüljön a mélyűrbe, be kell állítani a lépésméretet a süllyedés irányába. Általánosságban elmondható, hogy egy egész cikket írhat arról, hogyan állítsa be helyesen a lépést, és hogyan változtassa meg azt az ereszkedési folyamat során a számítási költségek csökkentése érdekében. Most azonban egy kicsit más feladat vár ránk, és a lépésméretet a „piszkálás” tudományos módszerével, vagy ahogy a köznyelvben mondják, empirikusan határozzuk meg.
Ha már a megadott koordinátáktól vagyunk и kivonjuk a deriváltak értékeit, új koordinátákat kapunk и . Megtesszük a következő lépést (kivonás), már a számított koordinátákból. Így a ciklus újra és újra elkezdődik, amíg el nem éri a kívánt konvergenciát.

Minden! Most készen állunk a Mariana-árok legmélyebb szurdokának felkutatására. Kezdjük el.

A gradiens süllyedés kódja

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

Lemerültünk a Mariana-árok legaljára, és ott ugyanazokat az együttható értékeket találtuk и , ami pontosan várható volt.

Vegyünk még egy merülést, csak ezúttal a mélytengeri járművünk tele lesz más technológiával, nevezetesen egy könyvtárral numpy.

A gradiens süllyedés kódja (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

Együttható értékek и megváltoztathatatlan.

Nézzük meg, hogyan változott a hiba a gradiens süllyedés során, vagyis hogyan változott az eltérések négyzetes összege az egyes lépéseknél.

Kód az eltérések négyzetes összegeinek ábrázolásához

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

4. grafikon „Szórás négyzetes összege a lejtős süllyedés során”

A grafikonon azt látjuk, hogy minden lépéssel csökken a hiba, és bizonyos számú iteráció után szinte vízszintes vonalat figyelünk meg.

Végül becsüljük meg a kódvégrehajtási idő különbségét:

Kód a gradiens süllyedés számítási idejének meghatározásához

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)

Lehet, hogy valamit rosszul csinálunk, de ez megint egy egyszerű „házilag írt” függvény, amely nem használja a könyvtárat numpy felülmúlja a függvénytárat használó függvény számítási idejét numpy.

De nem állunk egy helyben, hanem egy másik izgalmas módszer tanulmányozása felé haladunk az egyszerű lineáris regressziós egyenlet megoldására. Találkozz velünk!

Sztochasztikus gradiens süllyedés

A sztochasztikus gradiens süllyedés működési elvének gyors megértése érdekében jobb meghatározni annak különbségeit a szokásos gradiens süllyedéstől. Mi a gradiens leszármazás esetén a deriváltjainak egyenleteiben и a mintában elérhető összes jellemző értékének és igaz válaszának összegét használta (azaz az összes и ). A sztochasztikus gradiens süllyedésnél nem használjuk fel a mintában lévő összes értéket, hanem pszeudo-véletlenszerűen választjuk ki az úgynevezett mintaindexet, és használjuk annak értékeit.

Például, ha az indexet 3-asnak (háromnak) határozzuk meg, akkor az értékeket vesszük и , akkor az értékeket behelyettesítjük a derivált egyenletekbe, és új koordinátákat határozunk meg. Ezután a koordináták meghatározása után ismét pszeudo-véletlenszerűen meghatározzuk a minta indexét, az indexnek megfelelő értékeket behelyettesítjük a parciális differenciálegyenletekbe, és új módon határozzuk meg a koordinátákat. и stb. amíg a konvergencia zöldre nem vált. Első pillantásra úgy tűnik, hogy ez egyáltalán nem működik, de működik. Igaz, érdemes megjegyezni, hogy a hiba nem csökken minden lépéssel, de bizonyosan van tendencia.

Milyen előnyei vannak a sztochasztikus gradiens süllyedésnek a hagyományoshoz képest? Ha a mintánk nagyon nagy, és több tízezer értékben mérjük, akkor sokkal könnyebb feldolgozni, mondjuk, egy véletlenszerű ezret, mint a teljes mintát. Itt jön képbe a sztochasztikus gradiens süllyedés. A mi esetünkben természetesen nem sok különbséget fogunk észrevenni.

Nézzük a kódot.

Sztochasztikus gradiens süllyedés kódja

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

Alaposan megnézzük az együtthatókat, és azon kapjuk magunkat, hogy feltesszük a kérdést: „Hogy lehet ez?” Más együttható értékeket kaptunk и . Lehet, hogy a sztochasztikus gradiens süllyedés optimálisabb paramétereket talált az egyenlethez? Sajnos nincs. Elég megnézni a négyzetes eltérések összegét, és látni, hogy az együtthatók új értékeivel a hiba nagyobb. Nem sietünk a kétségbeesésbe. Készítsünk grafikont a hibaváltozásról.

Kód a sztochasztikus gradiens süllyedés négyzetes eltéréseinek összegének ábrázolásához

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

5. grafikon „A sztochasztikus gradiens süllyedés közbeni eltérések négyzetes összege”

A menetrendet nézve minden a helyére kerül, és most mindent kijavítunk.

Szóval mi történt? A következő történt. Ha véletlenszerűen választunk ki egy hónapot, akkor algoritmusunk a kiválasztott hónapra igyekszik csökkenteni a bevételszámítási hibát. Ezután kiválasztunk egy másik hónapot, és megismételjük a számítást, de csökkentjük a hibát a második kiválasztott hónapra. Most ne feledje, hogy az első két hónap jelentősen eltér az egyszerű lineáris regressziós egyenlet vonalától. Ez azt jelenti, hogy e két hónap bármelyikének kiválasztásakor az algoritmusunk mindegyik hibájának csökkentésével jelentősen megnöveli a teljes minta hibáját. Szóval mit kéne tenni? A válasz egyszerű: csökkentenie kell az ereszkedési lépést. Hiszen az ereszkedési lépés csökkentésével a hiba is abbahagyja a fel-le „ugrálást”. Illetve az „ugró” hiba nem szűnik meg, de nem is teszi meg olyan gyorsan :) Nézzük.

Kód az SGD kisebb lépésekkel történő futtatásához

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

6. számú grafikon „Sztochasztikus gradiens süllyedés (80 ezer lépés) alatti eltérések négyzetes összege”

Az együtthatók javultak, de még mindig nem ideálisak. Hipotetikusan ez így korrigálható. Kiválasztjuk például az utolsó 1000 iterációban azoknak az együtthatóknak az értékeit, amelyekkel a minimális hiba történt. Igaz, ehhez maguknak az együtthatók értékét is fel kell írnunk. Ezt nem fogjuk megtenni, inkább az ütemezésre figyelünk. Simának tűnik, és a hiba egyenletesen csökken. Valójában ez nem igaz. Nézzük meg az első 1000 iterációt, és hasonlítsuk össze őket az utolsóval.

SGD diagram kódja (első 1000 lépés)

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

7. grafikon „SGD eltérések négyzetes összege (első 1000 lépés)”

8. grafikon „SGD eltérések négyzetes összege (utolsó 1000 lépés)”

Az ereszkedés legelején meglehetősen egyenletes és meredek hibacsökkenést figyelhetünk meg. Az utolsó iterációkban azt látjuk, hogy a hiba az 1,475-ös érték körül-körül mozog, és bizonyos pillanatokban még megegyezik ezzel az optimális értékkel, de akkor is felmegy... Ismétlem, felírhatja az együtthatók и , majd válassza ki azokat, amelyeknél a hiba minimális. Volt azonban egy komolyabb problémánk: 80 ezer lépést kellett megtennünk (lásd a kódot), hogy az optimálishoz közeli értékeket kapjunk. És ez már ellentmond annak az elképzelésnek, hogy a számítási időt megtakarítsuk a gradiens süllyedéshez viszonyított sztochasztikus gradiens süllyedéssel. Mit lehet javítani és javítani? Nem nehéz észrevenni, hogy az első iterációk során magabiztosan megyünk lefelé, ezért az első iterációkban nagy lépést kell hagynunk, és előrehaladva csökkenteni a lépést. Ebben a cikkben ezt nem tesszük meg – már így is túl hosszú. Aki szeretné, az elgondolkodhat magában, hogyan csinálja, nem nehéz :)

Most hajtsunk végre sztochasztikus gradiens süllyedést a könyvtár segítségével numpy (és ne botorkáljunk a korábban azonosított köveken)

Sztochasztikus gradiens süllyedés kódja (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

Az értékek majdnem ugyanazok, mint a használat nélküli ereszkedésnél numpy. Ez azonban logikus.

Nézzük meg, mennyi ideig tartott a sztochasztikus gradiens süllyedés.

Kód az SGD számítási idő meghatározásához (80 ezer lépés)

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)

Minél beljebb kerül az erdő, annál sötétebbek a felhők: ismét a „saját írta” képlet mutatja a legjobb eredményt. Mindez arra utal, hogy a könyvtár használatának még finomabb módjaira van szükség numpy, amelyek valóban felgyorsítják a számítási műveleteket. Ebben a cikkben nem fogunk tudni róluk. Lesz min gondolkodni a szabadidődben :)

Összefoglaljuk

Összegzés előtt egy olyan kérdésre szeretnék választ adni, amely nagy valószínűséggel kedves olvasónk részéről felmerült. Tulajdonképpen miért ilyen „kínzások” ereszkedésekkel, miért kell fel-alá járnunk a hegyen (leginkább lefelé), hogy megtaláljuk a kincses alföldet, ha egy ilyen erős és egyszerű eszköz van a kezünkben, a egy analitikus megoldás formájában, amely azonnal a megfelelő helyre teleportál minket?

A válasz erre a kérdésre a felszínen rejlik. Most egy nagyon egyszerű példát néztünk meg, amelyben az igazi válasz az egy jeltől függ . Az életben nem gyakran látni ilyet, ezért képzeljük el, hogy 2, 30, 50 vagy több jelünk van. Adjunk ehhez hozzá több ezer vagy akár tízezer értéket minden attribútumhoz. Ebben az esetben előfordulhat, hogy az analitikai oldat nem bírja ki a tesztet, és meghiúsul. A gradiens süllyedés és annak változatai viszont lassan, de biztosan közelebb visznek a célhoz - a függvény minimumához. És ne aggódjon a sebesség miatt – valószínűleg megvizsgáljuk azokat a módokat, amelyek lehetővé teszik a lépéshossz (vagyis a sebesség) beállítását és szabályozását.

És most a tényleges rövid összefoglaló.

Először is remélem, hogy a cikkben bemutatott anyag segít a kezdő „adattudósoknak” megérteni, hogyan kell megoldani az egyszerű (és nem csak) lineáris regressziós egyenleteket.

Másodszor, az egyenlet megoldásának több módját is megvizsgáltuk. Most a helyzettől függően kiválaszthatjuk a probléma megoldására legalkalmasabbat.

Harmadszor, láttuk a további beállítások erejét, nevezetesen a gradiens süllyedés lépéshosszát. Ez a paraméter nem elhanyagolható. Amint fentebb megjegyeztük, a számítások költségeinek csökkentése érdekében a lépéshosszt módosítani kell az ereszkedés során.

Negyedszer, esetünkben az „otthon írt” függvények mutatták a legjobb időeredményeket a számításokhoz. Ennek oka valószínűleg a könyvtár lehetőségeinek nem a legprofibb kihasználása numpy. De akárhogy is legyen, a következő következtetés önmagát sugallja. Egyrészt néha érdemes megkérdőjelezni a kialakult véleményeket, másrészt nem mindig érdemes mindent bonyolítani - éppen ellenkezőleg, néha egy egyszerűbb problémamegoldási mód hatékonyabb. És mivel a célunk egy egyszerű lineáris regressziós egyenlet megoldásának három megközelítésének elemzése volt, a „saját írt” függvények használata bőven elég volt számunkra.