🥇De vergelijking van eenvoudige lineaire regressie oplossen

Het artikel bespreekt verschillende manieren om de wiskundige vergelijking van een eenvoudige (gepaarde) regressielijn te bepalen.

Alle methoden voor het oplossen van de hier besproken vergelijking zijn gebaseerd op de kleinste kwadratenmethode. Laten we de methoden als volgt aanduiden:

Analytische oplossing
Gradiënt afdaling
Stochastische gradiëntafdaling

Voor elk van de methoden voor het oplossen van de vergelijking van een rechte lijn biedt het artikel verschillende functies, die voornamelijk zijn onderverdeeld in functies die zijn geschreven zonder de bibliotheek te gebruiken NumPy en degenen die gebruiken voor berekeningen NumPy. Er wordt aangenomen dat vakkundig gebruik NumPy zal de computerkosten verlagen.

Alle code in het artikel is in de taal geschreven python 2.7 met Jupyter Notebook. De broncode en het bestand met voorbeeldgegevens worden op geplaatst Github

Het artikel is meer gericht op zowel beginners als op degenen die al geleidelijk begonnen zijn met het beheersen van de studie van een zeer breed gedeelte van de kunstmatige intelligentie: machinaal leren.

Om de stof te illustreren gebruiken we een heel eenvoudig voorbeeld.

Voorbeeld voorwaarden

We hebben vijf waarden die afhankelijkheid kenmerken Y van X (Tabel nr. 1):

Tabel nr. 1 “Voorbeeldvoorwaarden”

We gaan ervan uit dat de waarden is de maand van het jaar, en — omzet deze maand. Met andere woorden, de omzet is afhankelijk van de maand van het jaar, en - het enige teken waarvan de inkomsten afhankelijk zijn.

Het voorbeeld is zo, zowel vanuit het oogpunt van de voorwaardelijke afhankelijkheid van de omzet van de maand van het jaar, als vanuit het oogpunt van het aantal waarden - er zijn er maar heel weinig. Een dergelijke vereenvoudiging zal het echter mogelijk maken, zoals ze zeggen, om, niet altijd met gemak, het materiaal uit te leggen dat beginners assimileren. En ook de eenvoud van de cijfers zal degenen die het voorbeeld op papier willen oplossen, mogelijk maken zonder noemenswaardige arbeidskosten.

Laten we aannemen dat de in het voorbeeld gegeven afhankelijkheid vrij goed kan worden benaderd door de wiskundige vergelijking van een eenvoudige (gepaarde) regressielijn van de vorm:

waar is de maand waarin de inkomsten zijn ontvangen, — inkomsten die overeenkomen met de maand, и zijn de regressiecoëfficiënten van de geschatte lijn.

Merk op dat de coëfficiënt vaak de helling of gradiënt van de geschatte lijn genoemd; vertegenwoordigt het bedrag waarmee de wanneer het verandert .

Uiteraard is het onze taak in het voorbeeld om dergelijke coëfficiënten in de vergelijking te selecteren и , waarbij de afwijkingen van onze berekende omzetwaarden per maand van de echte antwoorden, d.w.z. waarden gepresenteerd in het voorbeeld zullen minimaal zijn.

Kleinste kwadratenmethode

Volgens de kleinste kwadratenmethode moet de afwijking worden berekend door deze te kwadrateren. Met deze techniek kunt u voorkomen dat afwijkingen wederzijds worden geannuleerd als ze tegengestelde tekenen hebben. Als er in één geval bijvoorbeeld een afwijking is +5 (plus vijf), en in de andere -5 (min vijf), dan heffen de som van de afwijkingen elkaar op en bedragen 0 (nul). Het is mogelijk om de afwijking niet te kwadrateren, maar om de eigenschap van de modulus te gebruiken, en dan zullen alle afwijkingen positief zijn en zich ophopen. We zullen niet in detail op dit punt ingaan, maar eenvoudigweg aangeven dat het voor het gemak van berekeningen gebruikelijk is om de afwijking te kwadrateren.

Zo ziet de formule eruit waarmee we de kleinste som van kwadratische afwijkingen (fouten) gaan bepalen:

waar is een functie van de benadering van ware antwoorden (dat wil zeggen, de inkomsten die we hebben berekend),

zijn de echte antwoorden (inkomsten vermeld in de steekproef),

is de steekproefindex (nummer van de maand waarin de afwijking wordt vastgesteld)

Laten we de functie differentiëren, de partiële differentiaalvergelijkingen definiëren en klaar zijn om verder te gaan met de analytische oplossing. Maar laten we eerst een korte excursie maken over wat differentiatie is en de geometrische betekenis van de afgeleide onthouden.

Differentiatie

Differentiatie is de bewerking waarbij de afgeleide van een functie wordt gevonden.

Waar wordt het derivaat voor gebruikt? De afgeleide van een functie karakteriseert de veranderingssnelheid van de functie en vertelt ons de richting ervan. Als de afgeleide op een bepaald punt positief is, neemt de functie toe; anders neemt de functie af. En hoe groter de waarde van de absolute afgeleide, hoe hoger de veranderingssnelheid van de functiewaarden, en hoe steiler de helling van de functiegrafiek.

Onder de omstandigheden van een Cartesiaans coördinatensysteem is de waarde van de afgeleide op het punt M(0,0) bijvoorbeeld gelijk aan + 25 betekent dat op een gegeven moment de waarde wordt verschoven naar rechts door een conventionele eenheid, waarde stijgt met 25 conventionele eenheden. Op de grafiek lijkt het een vrij steile stijging van de waarden vanaf een bepaald punt.

Een ander voorbeeld. De afgeleide waarde is gelijk -0,1 betekent dat wanneer ontheemd per conventionele eenheid, waarde neemt af met slechts 0,1 conventionele eenheid. Tegelijkertijd kunnen we op de grafiek van de functie een nauwelijks merkbare neerwaartse helling waarnemen. Als we een analogie maken met een berg, is het alsof we heel langzaam een zachte helling afdalen vanaf een berg, in tegenstelling tot het vorige voorbeeld, waar we hele steile toppen moesten beklimmen :)

Dus na het differentiëren van de functie door overmacht и definiëren we partiële differentiaalvergelijkingen van de eerste orde. Na het bepalen van de vergelijkingen, zullen we een systeem van twee vergelijkingen ontvangen, door deze op te lossen zullen we dergelijke waarden van de coëfficiënten kunnen selecteren и , waarbij de waarden van de overeenkomstige afgeleiden op bepaalde punten met een zeer, zeer kleine hoeveelheid veranderen, en in het geval van een analytische oplossing helemaal niet. Met andere woorden, de foutfunctie op de gevonden coëfficiënten zal een minimum bereiken, omdat de waarden van de partiële afgeleiden op deze punten gelijk zullen zijn aan nul.

Dus, volgens de differentiatieregels, de partiële afgeleide vergelijking van de eerste orde met betrekking tot de coëfficiënt zal de vorm aannemen:

Partiële afgeleide vergelijking van de eerste orde met betrekking tot zal de vorm aannemen:

Als gevolg hiervan ontvingen we een systeem van vergelijkingen met een vrij eenvoudige analytische oplossing:

begin{vergelijking*}
begin{gevallen}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0

somlimieten_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — somlimieten_{i=1}^ny_i) = 0
eind{gevallen}
einde{vergelijking*}

Voordat we de vergelijking oplossen, gaan we eerst laden, controleren of het laden correct is en de gegevens formatteren.

Gegevens laden en formatteren

Opgemerkt moet worden dat vanwege het feit dat we voor de analytische oplossing, en vervolgens voor gradiënt en stochastische gradiëntafdaling, de code in twee varianten zullen gebruiken: met behulp van de bibliotheek NumPy en zonder het te gebruiken, hebben we de juiste gegevensopmaak nodig (zie code).

Code voor het laden en verwerken van gegevens

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

Visualisatie

Nu we eerst de gegevens hebben geladen, ten tweede de juistheid van het laden hebben gecontroleerd en uiteindelijk de gegevens hebben geformatteerd, zullen we de eerste visualisatie uitvoeren. De methode die hiervoor vaak wordt gebruikt is paarplot Bibliotheek zeegeborene. In ons voorbeeld heeft het vanwege de beperkte aantallen geen zin om de bibliotheek te gebruiken zeegeborene. Wij maken gebruik van de reguliere bibliotheek matplotlib en kijk maar eens naar het spreidingsdiagram.

Scatterplot-code

print 'График №1 "Зависимость выручки от месяца года"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()

Grafiek nr. 1 “Afhankelijkheid van de omzet van de maand van het jaar”

Analytische oplossing

Laten we de meest voorkomende hulpmiddelen gebruiken python en los het stelsel vergelijkingen op:

begin{vergelijking*}
begin{gevallen}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0

somlimieten_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — somlimieten_{i=1}^ny_i) = 0
eind{gevallen}
einde{vergelijking*}

Volgens de regel van Cramer we zullen de algemene determinant vinden, evenals de determinanten door en , waarna de determinant wordt gedeeld door naar de algemene determinant - vind de coëfficiënt , op dezelfde manier vinden we de coëfficiënt .

Analytische oplossingscode

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

Dit is wat we hebben:

Dus de waarden van de coëfficiënten zijn gevonden, de som van de kwadratische afwijkingen is vastgesteld. Laten we een rechte lijn tekenen op het verstrooiingshistogram in overeenstemming met de gevonden coëfficiënten.

Regressielijncode

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

Grafiek nr. 2 ‘Juiste en berekende antwoorden’

U kunt voor elke maand de afwijkingsgrafiek bekijken. In ons geval zullen we er geen significante praktische waarde aan ontlenen, maar we zullen onze nieuwsgierigheid bevredigen over hoe goed de eenvoudige lineaire regressievergelijking de afhankelijkheid van de omzet van de maand van het jaar karakteriseert.

Afwijkingskaartcode

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

Grafiek nr. 3 “Afwijkingen, %”

Niet perfect, maar we hebben onze taak volbracht.

Laten we een functie schrijven die de coëfficiënten bepaalt и maakt gebruik van de bibliotheek NumPyOm preciezer te zijn, we zullen twee functies schrijven: één met behulp van een pseudo-inverse matrix (in de praktijk niet aanbevolen, omdat het proces computationeel complex en onstabiel is), de andere met behulp van een matrixvergelijking.

Analytische oplossingscode (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

Laten we de tijd vergelijken die is besteed aan het bepalen van de coëfficiënten и , in overeenstemming met de 3 gepresenteerde methoden.

Code voor het berekenen van de rekentijd

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)

Met een kleine hoeveelheid gegevens komt er een ‘zelfgeschreven’ functie naar voren, die de coëfficiënten vindt met behulp van de methode van Cramer.

Nu kunt u verder gaan met andere manieren om coëfficiënten te vinden и .

Gradiënt afdaling

Laten we eerst definiëren wat een verloop is. Simpel gezegd is de gradiënt een segment dat de richting van maximale groei van een functie aangeeft. Naar analogie met het beklimmen van een berg: waar de hellingshoek zich bevindt, bevindt zich de steilste klim naar de top van de berg. Als we het voorbeeld met de berg ontwikkelen, herinneren we ons dat we in feite de steilste afdaling nodig hebben om het laagland zo snel mogelijk te bereiken, dat wil zeggen het minimum - de plaats waar de functie niet toeneemt of afneemt. Op dit punt is de afgeleide gelijk aan nul. Daarom hebben we geen gradiënt nodig, maar een antigradiënt. Om de antigradiënt te vinden hoef je alleen maar de gradiënt te vermenigvuldigen met -1 (min een).

Laten we aandacht besteden aan het feit dat een functie meerdere minima kan hebben, en als we naar een van deze minima zijn afgedaald met behulp van het hieronder voorgestelde algoritme, zullen we geen ander minimum kunnen vinden, dat mogelijk lager is dan het gevonden minimum. Laten we ontspannen, dit is geen bedreiging voor ons! In ons geval hebben we te maken met één enkel minimum, aangezien onze functie in de grafiek staat een regelmatige parabool. En zoals we allemaal heel goed zouden moeten weten uit onze wiskundecursus op school, heeft een parabool slechts één minimum.

Nadat we ontdekten waarom we een gradiënt nodig hadden, en ook dat de gradiënt een segment is, dat wil zeggen een vector met gegeven coördinaten, die precies dezelfde coëfficiënten zijn и we kunnen gradiëntafdaling implementeren.

Voordat ik begin, stel ik voor om slechts een paar zinnen over het afdalingsalgoritme te lezen:

We bepalen op een pseudo-willekeurige manier de coördinaten van de coëfficiënten и . In ons voorbeeld zullen we coëfficiënten definiëren die dichtbij nul liggen. Dit is een gebruikelijke praktijk, maar elk geval kan zijn eigen praktijk hebben.
Van coördinaat trek de waarde van de partiële afgeleide van de eerste orde op het punt af . Dus als de afgeleide positief is, neemt de functie toe. Daarom zullen we, door de waarde van de afgeleide af te trekken, in de tegenovergestelde richting van de groei bewegen, dat wil zeggen in de richting van afdaling. Als de afgeleide negatief is, neemt de functie op dit punt af en door de waarde van de afgeleide af te trekken gaan we in de richting van de afdaling.
We voeren een soortgelijke operatie uit met de coördinaat : trek de waarde van de partiële afgeleide op het punt af .
Om niet over het minimum te springen en de diepe ruimte in te vliegen, is het noodzakelijk om de stapgrootte in de richting van de afdaling in te stellen. Over het algemeen zou u een heel artikel kunnen schrijven over hoe u de stap correct instelt en hoe u deze tijdens het afdalingsproces kunt wijzigen om de rekenkosten te verlagen. Maar nu hebben we een iets andere taak voor de boeg, en we zullen de stapgrootte vaststellen met behulp van de wetenschappelijke methode van "porren" of, zoals ze in het gewone taalgebruik zeggen, empirisch.
Zodra we van de opgegeven coördinaten zijn и trek de waarden van de afgeleiden af, we krijgen nieuwe coördinaten и . We nemen de volgende stap (aftrekken), al vanaf de berekende coördinaten. En zo begint de cyclus steeds opnieuw, totdat de vereiste convergentie is bereikt.

Alle! Nu zijn we klaar om op zoek te gaan naar de diepste kloof van de Marianentrog. Laten we beginnen.

Codeer voor gradiëntafdaling

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

We doken naar de bodem van de Marianentrog en daar vonden we allemaal dezelfde coëfficiëntwaarden и , en dat was precies wat te verwachten was.

Laten we nog een duik nemen, maar deze keer zal ons diepzeevoertuig gevuld zijn met andere technologieën, namelijk een bibliotheek NumPy.

Code voor gradiëntafdaling (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

Coëfficiënt waarden и onveranderlijk.

Laten we eens kijken hoe de fout veranderde tijdens de gradiëntdaling, dat wil zeggen hoe de som van de kwadratische afwijkingen bij elke stap veranderde.

Code voor het plotten van de som van gekwadrateerde afwijkingen

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

Grafiek nr. 4 “Som van kwadratische afwijkingen tijdens gradiëntdaling”

In de grafiek zien we dat bij elke stap de fout afneemt, en na een bepaald aantal iteraties zien we een bijna horizontale lijn.

Laten we tot slot het verschil in code-uitvoeringstijd schatten:

Code om de berekeningstijd van de gradiëntdaling te bepalen

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)

Misschien doen we iets verkeerd, maar nogmaals, het is een eenvoudige “zelfgeschreven” functie die geen gebruik maakt van de bibliotheek NumPy presteert beter dan de berekeningstijd van een functie die de bibliotheek gebruikt NumPy.

Maar we staan niet stil, maar zijn op weg naar het bestuderen van een andere opwindende manier om de eenvoudige lineaire regressievergelijking op te lossen. Ontmoeten!

Stochastische gradiëntafdaling

Om het werkingsprincipe van stochastische gradiëntdaling snel te begrijpen, is het beter om de verschillen met gewone gradiëntdaling te bepalen. Wij, in het geval van gradiëntafdaling, in de vergelijkingen van afgeleiden van и gebruikte de som van de waarden van alle kenmerken en echte antwoorden die beschikbaar zijn in de steekproef (dat wil zeggen de som van alles и ). Bij stochastische gradiëntafdaling zullen we niet alle waarden gebruiken die in de steekproef aanwezig zijn, maar in plaats daarvan pseudo-willekeurig de zogenaamde voorbeeldindex selecteren en de waarden ervan gebruiken.

Als bijvoorbeeld wordt vastgesteld dat de index nummer 3 (drie) is, nemen we de waarden и , dan vervangen we de waarden in de afgeleide vergelijkingen en bepalen we nieuwe coördinaten. Vervolgens bepalen we, nadat we de coördinaten hebben bepaald, opnieuw pseudo-willekeurig de voorbeeldindex, vervangen we de waarden die overeenkomen met de index in de partiële differentiaalvergelijkingen en bepalen we de coördinaten op een nieuwe manier и enz. totdat de convergentie groen wordt. Op het eerste gezicht lijkt het misschien niet zo dat dit helemaal zou kunnen werken, maar dat is wel zo. Het is waar dat het de moeite waard is om op te merken dat de fout niet bij elke stap afneemt, maar er is zeker een tendens.

Wat zijn de voordelen van stochastische gradiëntdaling ten opzichte van conventionele afdalingen? Als onze steekproefomvang erg groot is en wordt gemeten in tienduizenden waarden, dan is het veel gemakkelijker om bijvoorbeeld een willekeurige duizend van deze waarden te verwerken, in plaats van de hele steekproef. Dit is waar stochastische gradiëntafdaling een rol speelt. In ons geval zullen we uiteraard niet veel verschil merken.

Laten we naar de code kijken.

Codeer voor stochastische gradiëntafdaling

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

We kijken zorgvuldig naar de coëfficiënten en betrappen onszelf erop dat we de vraag stellen: “Hoe kan dit?” We hebben andere coëfficiëntwaarden и . Misschien heeft de stochastische gradiëntafdaling meer optimale parameters voor de vergelijking gevonden? Helaas niet. Het volstaat om naar de som van de kwadratische afwijkingen te kijken en te zien dat bij nieuwe waarden van de coëfficiënten de fout groter is. Wij hebben geen haast om te wanhopen. Laten we een grafiek maken van de foutverandering.

Code voor het uitzetten van de som van kwadratische afwijkingen in stochastische gradiëntafdaling

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

Grafiek nr. 5 “Som van kwadratische afwijkingen tijdens stochastische gradiëntdaling”

Kijkend naar het schema valt alles op zijn plaats en nu gaan we alles in orde maken.

Dus wat gebeurde er? Het volgende gebeurde. Wanneer we willekeurig een maand selecteren, probeert ons algoritme voor de geselecteerde maand de fout bij het berekenen van de omzet te verminderen. Vervolgens selecteren we een andere maand en herhalen we de berekening, maar we verminderen de fout voor de tweede geselecteerde maand. Bedenk nu dat de eerste twee maanden aanzienlijk afwijken van de lijn van de eenvoudige lineaire regressievergelijking. Dit betekent dat wanneer een van deze twee maanden wordt geselecteerd, ons algoritme de fout voor de hele steekproef aanzienlijk vergroot, door de fout van elk van deze maanden te verkleinen. Dus wat te doen? Het antwoord is simpel: je moet de afdalingsstap verkleinen. Door de afdalingsstap te verkleinen, stopt de fout immers ook met het op en neer ‘springen’. Of beter gezegd, de “springende” fout stopt niet, maar het zal niet zo snel gebeuren :) Laten we eens kijken.

Code om SGD met kleinere stappen uit te voeren

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

Grafiek nr. 6 “Som van kwadratische afwijkingen tijdens stochastische gradiëntdaling (80 stappen)”

De coëfficiënten zijn verbeterd, maar zijn nog steeds niet ideaal. Hypothetisch gezien kan dit op deze manier worden gecorrigeerd. We selecteren bijvoorbeeld in de laatste 1000 iteraties de waarden van de coëfficiënten waarmee de minimale fout is gemaakt. Toegegeven, hiervoor zullen we ook de waarden van de coëfficiënten zelf moeten opschrijven. Wij zullen dit niet doen, maar wel op het schema letten. Het ziet er soepel uit en de fout lijkt gelijkmatig af te nemen. Eigenlijk is dit niet waar. Laten we naar de eerste 1000 iteraties kijken en deze vergelijken met de laatste.

Code voor SGD-grafiek (eerste 1000 stappen)

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

Grafiek nr. 7 “Som van kwadratische afwijkingen SGD (eerste 1000 stappen)”

Grafiek nr. 8 “Som van kwadratische afwijkingen SGD (laatste 1000 stappen)”

Helemaal aan het begin van de afdaling zien we een tamelijk uniforme en steile afname van de fouten. In de laatste iteraties zien we dat de fout rond en rond de waarde van 1,475 gaat en op sommige momenten zelfs gelijk is aan deze optimale waarde, maar dan gaat hij nog steeds omhoog... Ik herhaal, je kunt de waarden van de coëfficiënten и en selecteer vervolgens de fouten waarbij de fout minimaal is. We hadden echter een ernstiger probleem: we moesten 80 stappen zetten (zie code) om waarden bijna optimaal te krijgen. En dit is al in tegenspraak met het idee om rekentijd te besparen met stochastische gradiëntdaling ten opzichte van gradiëntdaling. Wat kan worden gecorrigeerd en verbeterd? Het is niet moeilijk om op te merken dat we in de eerste iteraties vol vertrouwen naar beneden gaan en dat we daarom in de eerste iteraties een grote stap moeten laten en de stap moeten verkleinen naarmate we verder komen. We zullen dit in dit artikel niet doen - het is al te lang. Degenen die dat willen, kunnen zelf bedenken hoe ze dit moeten doen, het is niet moeilijk :)

Laten we nu een stochastische gradiëntafdaling uitvoeren met behulp van de bibliotheek NumPy (en laten we niet struikelen over de stenen die we eerder hebben geïdentificeerd)

Code voor stochastische gradiëntdaling (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

De waarden bleken vrijwel hetzelfde als bij het afdalen zonder gebruik NumPy. Dit is echter logisch.

Laten we eens kijken hoe lang stochastische gradiëntdalingen ons hebben gekost.

Code voor het bepalen van de SGD-berekeningstijd (80 stappen)

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)

Hoe verder het bos in, hoe donkerder de wolken: ook hier laat de ‘zelfgeschreven’ formule het beste resultaat zien. Dit alles suggereert dat er nog subtielere manieren moeten zijn om de bibliotheek te gebruiken NumPy, wat de rekenbewerkingen echt versnelt. In dit artikel zullen we er niets over leren. Er is vast wel iets om over na te denken in je vrije tijd :)

We vatten samen

Voordat ik samenvat, wil ik een vraag beantwoorden die hoogstwaarschijnlijk is ontstaan door onze geachte lezer. Waarom eigenlijk zo'n "marteling" met afdalingen, waarom moeten we de berg op en af lopen (meestal naar beneden) om het gekoesterde laagland te vinden, als we zo'n krachtig en eenvoudig apparaat in onze handen hebben, in de vorm van een analytische oplossing, die ons onmiddellijk naar de juiste plaats teleporteert?

Het antwoord op deze vraag ligt aan de oppervlakte. Nu hebben we naar een heel eenvoudig voorbeeld gekeken, waarin het ware antwoord luidt: hangt af van één teken . Je ziet dit niet vaak in het leven, dus laten we ons voorstellen dat we 2, 30, 50 of meer tekens hebben. Laten we hieraan duizenden of zelfs tienduizenden waarden toevoegen voor elk attribuut. In dit geval is het mogelijk dat de analytische oplossing de test niet doorstaat en faalt. Op zijn beurt zullen de gradiëntafdaling en zijn variaties ons langzaam maar zeker dichter bij het doel brengen: het minimum van de functie. En maak je geen zorgen over de snelheid: we zullen waarschijnlijk kijken naar manieren waarmee we de staplengte (dat wil zeggen snelheid) kunnen instellen en regelen.

En nu de eigenlijke korte samenvatting.

Ten eerste hoop ik dat het materiaal dat in dit artikel wordt gepresenteerd, beginnende ‘datawetenschappers’ zal helpen begrijpen hoe ze eenvoudige (en niet alleen) lineaire regressievergelijkingen kunnen oplossen.

Ten tweede hebben we verschillende manieren bekeken om de vergelijking op te lossen. Afhankelijk van de situatie kunnen we nu degene kiezen die het meest geschikt is om het probleem op te lossen.

Ten derde zagen we de kracht van aanvullende instellingen, namelijk de staplengte van de gradiëntafdaling. Deze parameter kan niet worden verwaarloosd. Zoals hierboven opgemerkt, moet, om de kosten van berekeningen te verminderen, de staplengte tijdens de afdaling worden gewijzigd.

Ten vierde lieten in ons geval ‘zelfgeschreven’ functies de beste tijdresultaten voor berekeningen zien. Dit komt waarschijnlijk doordat er niet op de meest professionele manier gebruik wordt gemaakt van de mogelijkheden van de bibliotheek NumPy. Maar hoe het ook zij, de volgende conclusie dringt zich op. Aan de ene kant is het soms de moeite waard om gevestigde meningen in twijfel te trekken, en aan de andere kant is het niet altijd de moeite waard om alles ingewikkeld te maken - integendeel, soms is een eenvoudigere manier om een probleem op te lossen effectiever. En aangezien het ons doel was om drie benaderingen te analyseren voor het oplossen van een eenvoudige lineaire regressievergelijking, was het gebruik van 'zelfgeschreven' functies voor ons voldoende.