🥇Løsning af ligningen for simpel lineær regression

Artiklen diskuterer flere måder at bestemme den matematiske ligning for en simpel (parret) regressionslinje.

Alle metoder til at løse ligningen diskuteret her er baseret på mindste kvadraters metode. Lad os betegne metoderne som følger:

Analytisk løsning
Gradient nedstigning
Stokastisk gradientnedstigning

For hver metode til at løse en lige linjes ligning giver artiklen forskellige funktioner, som hovedsageligt er opdelt i dem, der er skrevet uden brug af biblioteket nusset og dem, der bruges til beregninger nusset. Считается, что умелое использование nusset vil reducere computeromkostningerne.

Весь код, приведенный в статье, написан на языке python 2.7 med Jupyter Notebook. Kildekoden og filen med eksempeldata er lagt på Github

Artiklen henvender sig mere til både begyndere og dem, der allerede så småt er begyndt at mestre studiet af et meget bredt afsnit i kunstig intelligens – maskinlæring.

For at illustrere materialet bruger vi et meget simpelt eksempel.

Eksempelforhold

Vi har fem værdier, der karakteriserer afhængighed Y fra X (Tabel nr. 1):

Tabel nr. 1 "Eksempler på forhold"

Vi vil antage, at værdierne er årets måned, og — omsætning i denne måned. Med andre ord afhænger omsætningen af årets måned, og - det eneste tegn, som indtægterne afhænger af.

Eksemplet er halvdårligt, både med hensyn til omsætningens betingede afhængighed af årets måned og ud fra antallet af værdier - der er meget få af dem. Men en sådan forenkling vil gøre det muligt, som de siger, at forklare, ikke altid med lethed, det materiale, som begyndere assimilerer. Og også tallenes enkelhed vil tillade dem, der ønsker at løse eksemplet på papir uden betydelige arbejdsomkostninger.

Предположим, что приведенная в примере зависимость, может быть достаточно хорошо аппроксимирована математическим уравнением линии простой (парной) регрессии вида:

где er den måned, hvor indtægten blev modtaget, — выручка, соответствующая месяцу, и er regressionskoefficienterne for den estimerede linje.

Отметим, что коэффициент ofte kaldet hældningen eller gradienten af den estimerede linje; repræsenterer det beløb, hvormed når det ændrer sig .

Det er klart, at vores opgave i eksemplet er at vælge sådanne koefficienter i ligningen и , hvor afvigelserne af vores beregnede omsætningsværdier pr. måned fra de sande svar, dvs. værdier præsenteret i prøven vil være minimale.

Mindste kvadratisk metode

В соответствии с методом наименьших квадратов, отклонение стоит рассчитывать, возводя его в квадрат. Подобный прием позволяет избежать взаимного погашения отклонений, в том случае, если они имеют противоположные знаки. Например, если в одном случае, отклонение составляет +5 (plus fem), og i den anden -5 (minus fem), så vil summen af afvigelserne ophæve hinanden og udgøre 0 (nul). Det er muligt ikke at kvadrere afvigelsen, men at bruge modulets egenskab og så vil alle afvigelserne være positive og akkumulere. Vi vil ikke dvæle ved dette punkt i detaljer, men blot indikere, at det af hensyn til beregningerne er sædvanligt at kvadrere afvigelsen.

Sådan ser formlen ud, hvormed vi vil bestemme den mindste sum af kvadrerede afvigelser (fejl):

где er en funktion af tilnærmelse af sande svar (det vil sige den indtjening, vi har beregnet),

er de sande svar (omsætning angivet i stikprøven),

er stikprøveindekset (nummer på den måned, hvor afvigelsen er bestemt)

Lad os differentiere funktionen, definere de partielle differentialligninger og være klar til at gå videre til den analytiske løsning. Men lad os først tage en kort udflugt om, hvad differentiering er, og huske den geometriske betydning af derivatet.

Differentiering

Differentiering er operationen med at finde den afledede af en funktion.

Hvad bruges derivatet til? Den afledede af en funktion karakteriserer funktionens ændringshastighed og fortæller os dens retning. Hvis den afledede på et givet punkt er positiv, så stiger funktionen, ellers falder funktionen. Og jo større værdien af den absolutte afledte er, jo højere er ændringshastigheden af funktionsværdierne, såvel som jo stejlere hældning af funktionsgrafen.

For eksempel, under betingelserne for et kartesisk koordinatsystem, er værdien af den afledte i punktet M(0,0) lig med + 25 betyder, at på et givet tidspunkt, når værdien forskydes til højre af en konventionel enhed, værdi stiger med 25 konventionelle enheder. På grafen ligner det en ret stejl stigning i værdierne fra et givet punkt.

Et andet eksempel. Den afledte værdi er lig -0,1 betyder, at når fortrængt pr. en konventionel enhed, værdi falder kun med 0,1 konventionel enhed. Samtidig kan vi på funktionens graf observere en knap mærkbar nedadgående hældning. Når man tegner en analogi med et bjerg, er det, som om vi meget langsomt går ned ad en svag skråning fra et bjerg, i modsætning til det forrige eksempel, hvor vi skulle bestige meget stejle tinder :)

Således efter differentiering af funktionen efter odds и , definerer vi 1. ordens partielle differentialligninger. Efter at have bestemt ligningerne, vil vi modtage et system af to ligninger, ved at løse hvilke vi vil være i stand til at vælge sådanne værdier af koefficienterne и , for hvilke værdierne af de tilsvarende derivater på givne punkter ændres meget, meget lille, og i tilfælde af en analytisk løsning ændres slet ikke. Med andre ord vil fejlfunktionen ved de fundne koefficienter nå et minimum, da værdierne af de partielle afledte på disse punkter vil være lig med nul.

Så ifølge reglerne for differentiering er den partielle afledte ligning af 1. orden med hensyn til koefficienten vil tage formen:

1. ordens partiel afledt ligning mhp vil tage formen:

Som et resultat modtog vi et ligningssystem, der har en ret simpel analytisk løsning:

begynde{ligning*}
begynde{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
slut{cases}
slut{ligning*}

Før vi løser ligningen, lad os forudindlæse, kontrollere, at indlæsningen er korrekt, og formatere dataene.

Indlæsning og formatering af data

Det skal bemærkes, at på grund af det faktum, at for den analytiske løsning, og efterfølgende for gradient og stokastisk gradientnedstigning, vil vi bruge koden i to variationer: ved at bruge biblioteket nusset og uden at bruge det, så har vi brug for passende dataformatering (se kode).

Dataindlæsning og behandlingskode

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

Visualisering

Nu, efter at vi for det første har indlæst dataene, for det andet kontrolleret korrektheden af indlæsningen og til sidst formateret dataene, vil vi udføre den første visualisering. Metoden der ofte bruges til dette er parplot bibliotek Søfødt. I vores eksempel er der på grund af det begrænsede antal ingen mening i at bruge biblioteket Søfødt. Мы воспользуемся обычной библиотекой Matplotlib og se bare på scatterplotten.

Scatterplot kode

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

Figur nr. 1 "Afhængighed af omsætning på årets måned"

Analytisk løsning

Lad os bruge de mest almindelige værktøjer i python og løse ligningssystemet:

begynde{ligning*}
begynde{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
slut{cases}
slut{ligning*}

Ifølge Cramers regel vi finder den generelle determinant, samt determinanter ved og , hvorefter determinanten divideres med til den generelle determinant - find koefficienten , på samme måde finder vi koefficienten .

Kode for analytisk løsning

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

Her er hvad vi fik:

Så værdierne af koefficienterne er fundet, summen af de kvadrerede afvigelser er blevet etableret. Lad os tegne en lige linje på spredningshistogrammet i overensstemmelse med de fundne koefficienter.

Regressionslinjekode

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

Diagram nr. 2 "Korrekte og beregnede svar"

Du kan se på afvigelsesgrafen for hver måned. I vores tilfælde vil vi ikke udlede nogen væsentlig praktisk værdi af det, men vi vil tilfredsstille vores nysgerrighed om, hvor godt den simple lineære regressionsligning karakteriserer omsætningens afhængighed af måneden i året.

Afvigelsesdiagramkode

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

Diagram nr. 3 "Afvigelser, %"

Ikke perfekt, men vi fuldførte vores opgave.

Lad os skrive en funktion, der bestemmer koefficienterne и bruger biblioteket nusset, mere præcist vil vi skrive to funktioner: den ene bruger en pseudoinvers matrix (anbefales ikke i praksis, da processen er beregningsmæssigt kompleks og ustabil), den anden bruger en matrixligning.

Analytisk løsningskode (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

Сравним время, которое было затрачено на определение коэффициентов и , i overensstemmelse med de 3 præsenterede metoder.

Kode til beregning af beregningstid

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)

Med en lille mængde data kommer der en "selvskrevet" funktion frem, som finder koefficienterne ved hjælp af Cramers metode.

Nu kan du gå videre til andre måder at finde koefficienter på и .

Gradient nedstigning

Lad os først definere, hvad en gradient er. Kort sagt er gradienten et segment, der angiver retningen for maksimal vækst af en funktion. I analogi med at bestige et bjerg, hvor gradienten er der, hvor den stejleste stigning til toppen af bjerget er. Ved at udvikle eksemplet med bjerget husker vi, at vi faktisk har brug for den stejleste nedstigning for at nå lavlandet så hurtigt som muligt, det vil sige minimum - det sted, hvor funktionen ikke øges eller falder. På dette tidspunkt vil den afledede være lig nul. Derfor har vi ikke brug for en gradient, men en antigradient. For at finde antigradienten skal du bare gange gradienten med -1 (minus en).

Lad os være opmærksomme på, at en funktion kan have flere minima, og efter at have gået ned i en af dem ved hjælp af algoritmen foreslået nedenfor, vil vi ikke være i stand til at finde et andet minimum, som kan være lavere end det, der blev fundet. Lad os slappe af, dette er ikke en trussel mod os! I vores tilfælde har vi at gøre med et enkelt minimum, da vores funktion på grafen er en regulær parabel. Og som vi alle burde vide udmærket fra vores skolematematikkursus, har en parabel kun et minimum.

Efter at vi fandt ud af, hvorfor vi havde brug for en gradient, og også at gradienten er et segment, det vil sige en vektor med givne koordinater, som er præcis de samme koefficienter и vi kan implementere gradientnedstigning.

Перед запуском, предлагаю прочитать буквально несколько предложений об алгоритме спуска:

Vi bestemmer på en pseudo-tilfældig måde koefficienternes koordinater и . I vores eksempel vil vi definere koefficienter nær nul. Dette er en almindelig praksis, men hver sag kan have sin egen praksis.
От координаты вычитаем значение частной производной 1-го порядка в точке . Så hvis den afledede er positiv, så øges funktionen. Derfor vil vi, ved at trække værdien af den afledte, bevæge os i den modsatte retning af vækst, det vil sige i retning af nedstigning. Hvis den afledte er negativ, så falder funktionen på dette tidspunkt, og ved at trække værdien af den afledte bevæger vi os i nedstigningsretningen.
Vi udfører en lignende operation med koordinatet : trække værdien af den partielle afledte ved punktet .
For ikke at springe over minimumet og flyve ind i det dybe rum, er det nødvendigt at indstille trinstørrelsen i retning af nedstigning. Generelt kan du skrive en hel artikel om, hvordan du indstiller trinnet korrekt, og hvordan du ændrer det under nedstigningsprocessen for at reducere beregningsomkostningerne. Men nu har vi en lidt anden opgave foran os, og vi vil etablere trinstørrelsen ved hjælp af den videnskabelige metode "poke" eller, som man siger i almindelig sprogbrug, empirisk.
Når vi er fra de givne koordinater и trække værdierne af de afledte, får vi nye koordinater и . Vi tager det næste trin (subtraktion), allerede fra de beregnede koordinater. Og så starter cyklussen igen og igen, indtil den nødvendige konvergens er opnået.

Alle! Nu er vi klar til at gå på jagt efter den dybeste kløft i Marianergraven. Lad os komme igang.

Kode for gradient nedstigning

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

Vi dykkede helt ned til bunden af Marianergraven, og der fandt vi alle de samme koefficientværdier и , hvilket er præcis hvad man kunne forvente.

Lad os tage endnu et dyk, men denne gang vil vores dybhavsfartøj være fyldt med andre teknologier, nemlig et bibliotek nusset.

Kode for gradientnedstigning (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

Koefficientværdier и uforanderlig.

Lad os se på, hvordan fejlen ændrede sig under gradientnedstigning, det vil sige hvordan summen af kvadrerede afvigelser ændrede sig med hvert trin.

Kode til at plotte summer af kvadrerede afvigelser

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

Graf nr. 4 "Summen af kvadrerede afvigelser under gradientnedstigning"

På grafen ser vi, at for hvert trin aftager fejlen, og efter et vist antal iterationer observerer vi en næsten vandret linje.

Напоследок оценим разницу во времени исполнения кода:

Kode til at bestemme tid for beregning af gradientnedstigning

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)

Måske gør vi noget forkert, men igen er det en simpel "hjemmeskrevet" funktion, der ikke bruger biblioteket nusset udkonkurrerer beregningstiden for en funktion ved hjælp af biblioteket nusset.

Men vi står ikke stille, men er på vej mod at studere endnu en spændende måde at løse den simple lineære regressionsligning på. Mød os!

Stokastisk gradientnedstigning

For hurtigt at forstå princippet om drift af stokastisk gradientnedstigning er det bedre at bestemme dets forskelle fra almindelig gradientnedstigning. Vi, i tilfælde af gradient nedstigning, i ligningerne af afledte af и brugte summen af værdierne af alle funktioner og sande svar tilgængelige i prøven (det vil sige summen af alle и ). I stokastisk gradientnedstigning vil vi ikke bruge alle de værdier, der er til stede i prøven, men i stedet vælge pseudo-tilfældigt det såkaldte prøveindeks og bruge dets værdier.

For eksempel, hvis indekset er bestemt til at være nummer 3 (tre), så tager vi værdierne и , så erstatter vi værdierne i de afledte ligninger og bestemmer nye koordinater. Derefter, efter at have bestemt koordinaterne, bestemmer vi igen pseudo-tilfældigt prøveindekset, erstatter de værdier, der svarer til indekset i de partielle differentialligninger og bestemmer koordinaterne på en ny måde и etc. indtil konvergensen bliver grøn. Ved første øjekast ser det måske ikke ud til, at dette overhovedet kunne fungere, men det gør det. Det er rigtigt, at det er værd at bemærke, at fejlen ikke falder for hvert trin, men der er bestemt en tendens.

Hvad er fordelene ved stokastisk gradientnedstigning i forhold til konventionel? Hvis vores stikprøvestørrelse er meget stor og målt i titusindvis af værdier, så er det meget nemmere at behandle f.eks. et tilfældigt tusind af dem i stedet for hele stikprøven. Det er her stokastisk gradientnedstigning kommer i spil. I vores tilfælde vil vi selvfølgelig ikke mærke den store forskel.

Lad os se på koden.

Kode for stokastisk gradientnedstigning

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

Vi ser nøje på koefficienterne og fanger os selv i at stille spørgsmålet "Hvordan kan det være?" Vi har andre koefficientværdier и . Måske har stokastisk gradientnedstigning fundet mere optimale parametre for ligningen? Desværre ikke. Det er nok at se på summen af kvadrerede afvigelser og se, at med nye værdier af koefficienterne er fejlen større. Vi har ikke travlt med at fortvivle. Lad os bygge en graf over fejlændringen.

Код для графика суммы квадратов отклонений при стохастическом градиентном спуске

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

Graf nr. 5 "Summen af kvadrerede afvigelser under stokastisk gradientnedstigning"

Ser man på tidsplanen, falder alt på plads, og nu vil vi ordne det hele.

Hvad skete der? Følgende skete. Når vi tilfældigt vælger en måned, så er det for den valgte måned, at vores algoritme søger at reducere fejlen ved beregning af omsætning. Så vælger vi endnu en måned og gentager beregningen, men vi reducerer fejlen for den anden valgte måned. Lad os nu huske, at de første to måneder afviger væsentligt fra linjen i den simple lineære regressionsligning. Det betyder, at når en af disse to måneder vælges, ved at reducere fejlen for hver af dem, øger vores algoritme alvorligt fejlen for hele stikprøven. Så hvad skal man gøre? Svaret er enkelt: du skal reducere nedstigningstrinnet. Når alt kommer til alt, vil fejlen også stoppe med at "hoppe" op og ned ved at reducere nedstigningstrinnet. Eller rettere sagt, "hopping"-fejlen stopper ikke, men den vil ikke gøre det så hurtigt :) Lad os tjekke.

Kode til at køre SGD med mindre trin

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

Graf nr. 6 "Summen af kvadrerede afvigelser under stokastisk gradientnedstigning (80 tusind trin)"

Koefficienterne er forbedret, men er stadig ikke ideelle. Hypotetisk kan dette korrigeres på denne måde. Vi vælger f.eks. i de sidste 1000 iterationer værdierne af de koefficienter, som minimumsfejlen blev lavet med. Sandt nok, for dette bliver vi også nødt til at nedskrive værdierne af selve koefficienterne. Vi vil ikke gøre dette, men være opmærksomme på tidsplanen. Det ser glat ud, og fejlen ser ud til at falde jævnt. Faktisk er dette ikke sandt. Lad os se på de første 1000 iterationer og sammenligne dem med de sidste.

Kode til SGD-diagram (første 1000 trin)

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

Graf nr. 7 "Sum af kvadrerede afvigelser SGD (første 1000 trin)"

Graf nr. 8 "Sum af kvadrerede afvigelser SGD (sidste 1000 trin)"

Allerede i begyndelsen af nedstigningen observerer vi et ret ensartet og stejlt fald i fejl. I de sidste iterationer ser vi, at fejlen går rundt og rundt om værdien af 1,475 og i nogle øjeblikke endda er lig med denne optimale værdi, men så går den stadig op... Jeg gentager, du kan nedskrive værdierne af koefficienter и , og vælg derefter dem, hvor fejlen er minimal. Vi havde dog et mere alvorligt problem: vi skulle tage 80 tusinde trin (se kode) for at få værdier tæt på det optimale. Og dette er allerede i modstrid med ideen om at spare beregningstid med stokastisk gradientnedstigning i forhold til gradientnedstigning. Hvad kan rettes og forbedres? Det er ikke svært at bemærke, at i de første iterationer går vi selvsikkert ned, og derfor bør vi efterlade et stort skridt i de første iterationer og reducere trinnet, efterhånden som vi bevæger os fremad. Det vil vi ikke gøre i denne artikel - den er allerede for lang. De, der ønsker det, kan selv tænke over, hvordan man gør dette, det er ikke svært :)

Lad os nu udføre stokastisk gradientnedstigning ved hjælp af biblioteket nusset (og lad os ikke snuble over de sten, vi identificerede tidligere)

Kode for Stokastisk Gradient Descent (NumPy)

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

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

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

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


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



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

Værdierne viste sig at være næsten de samme, som når de falder uden brug nusset. Dette er dog logisk.

Lad os finde ud af, hvor lang tid det tog os med nedstigninger med stokastiske gradienter.

Kode til bestemmelse af SGD-beregningstid (80 tusind trin)

print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска без использования библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print '***************************************'
print

print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска с использованием библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

Jo længere ind i skoven, jo mørkere er skyerne: igen viser den "selvskrevne" formel det bedste resultat. Alt dette tyder på, at der skal være endnu mere subtile måder at bruge biblioteket på nusset, som virkelig fremskynder beregningsoperationer. I denne artikel lærer vi ikke om dem. Der vil være noget at tænke over i din fritid :)

Vi opsummerer

Før jeg opsummerer, vil jeg gerne svare på et spørgsmål, som højst sandsynligt er opstået fra vores kære læser. Hvorfor, i virkeligheden, sådan en "tortur" med nedstigninger, hvorfor skal vi gå op og ned ad bjerget (for det meste ned) for at finde det skattede lavland, hvis vi har en sådan kraftfuld og enkel enhed i hænderne form for en analytisk løsning, som øjeblikkeligt teleporterer os til det rigtige sted?

Svaret på dette spørgsmål ligger på overfladen. Nu har vi set på et meget simpelt eksempel, hvor det sande svar er afhænger af ét tegn . Du ser det ikke ofte i livet, så lad os forestille os, at vi har 2, 30, 50 eller flere tegn. Lad os tilføje tusinder eller endda titusinder af værdier for hver egenskab. I dette tilfælde kan den analytiske opløsning muligvis ikke modstå testen og mislykkes. Til gengæld vil gradientnedstigning og dens variationer langsomt, men sikkert bringe os tættere på målet - funktionens minimum. Og bekymre dig ikke om hastighed - vi vil sandsynligvis se på måder, der vil give os mulighed for at indstille og regulere skridtlængde (det vil sige hastighed).

Og nu den egentlige korte opsummering.

For det første håber jeg, at materialet præsenteret i artiklen vil hjælpe begyndende "dataforskere" med at forstå, hvordan man løser simple (og ikke kun) lineære regressionsligninger.

For det andet så vi på flere måder at løse ligningen på. Nu kan vi, afhængigt af situationen, vælge den, der er bedst egnet til at løse problemet.

For det tredje så vi styrken ved yderligere indstillinger, nemlig gradientens nedstigningstrinlængde. Denne parameter kan ikke ignoreres. Som nævnt ovenfor, for at reducere omkostningerne ved beregninger, bør trinlængden ændres under nedstigningen.

For det fjerde viste "hjemmeskrevne" funktioner i vores tilfælde de bedste tidsresultater til beregninger. Det skyldes formentlig ikke den mest professionelle brug af bibliotekets muligheder nusset. Men uanset hvad, så tyder følgende konklusion sig selv. På den ene side er det nogle gange værd at sætte spørgsmålstegn ved etablerede meninger, og på den anden side er det ikke altid værd at komplicere alt - tværtimod er en mere enkel måde at løse et problem på mere effektiv. Og da vores mål var at analysere tre tilgange til at løse en simpel lineær regressionsligning, var brugen af "selvskrevne" funktioner ganske nok for os.