Lösen der Gleichung der einfachen linearen Regression

Der Artikel betrachtet mehrere Methoden zur Bestimmung der mathematischen Gleichung einer einfachen (paarweisen) Regression.

Alle hier betrachteten Methoden zur Lösung der Gleichung basieren auf der Methode der kleinsten Quadrate. Bezeichnen wir die Methoden wie folgt:

• Analytische Lösung
• Gradientenabstieg
• Stochastischer Gradientenabstieg

Für jede Methode zur Lösung der Geradengleichung werden im Artikel verschiedene Funktionen vorgestellt, die hauptsächlich in zwei Kategorien unterteilt werden: diejenigen, die ohne die Verwendung einer Bibliothek geschrieben sind, NumPy und diejenigen, die für Berechnungen auf NumPyzurückgreifen. Es wird angenommen, dass eine geschickte Nutzung NumPy die Berechnungskosten senken kann.

Der gesamte im Artikel angegebene Code ist in der Sprache python 2.7 unter Verwendung von Jupyter Notebookgeschrieben. Der Quellcode und die Datei mit den Beispieldaten sind auf Github

verfügbar. Der Artikel richtet sich sowohl an Anfänger als auch an diejenigen, die bereits begonnen haben, das umfangreiche Thema im Bereich Künstliche Intelligenz – maschinelles Lernen – zu erkunden.

Zur Veranschaulichung des Materials verwenden wir ein sehr einfaches Beispiel.

Die Bedingungen des Beispiels

Wir haben fünf Werte, die die Abhängigkeit charakterisieren Y ab X (Tabelle Nr. 1):

Tabelle Nr. 1 „Bedingungen des Beispiels“

Wir nehmen an, dass die Werte – das Monat des Jahres sind und – die Einnahmen in diesem Monat. Mit anderen Worten, die Einnahmen sind abhängig vom Monat des Jahres, und – ist das einzige Merkmal, von dem die Einnahmen abhängen.

Das Beispiel ist nicht besonders gut, sowohl hinsichtlich der angenommenen Abhängigkeit der Einnahmen vom Monat des Jahres als auch hinsichtlich der Anzahl der Werte – sie sind sehr gering. Allerdings erlaubt diese Vereinfachung, das Material, das von Anfängern nicht immer leicht aufgenommen wird, anschaulich zu erklären. Die Einfachheit der Zahlen ermöglicht es auch, dass diejenigen, die möchten, das Beispiel ohne großen Aufwand „auf Papier“ lösen können.

Angenommen, die im Beispiel angegebene Abhängigkeit kann durch eine mathematische Gleichung der einfachen (paarweisen) Regression vom Typ gut approximiert werden:

wo – das ist der Monat, in dem die Einnahmen erzielt wurden, – die Einnahmen, die dem Monat entsprechen, und – die Koeffizienten der geschätzten Linie.

Es sei darauf hingewiesen, dass der Koeffizient häufig als Anstieg oder Gradienten der geschätzten Linie bezeichnet wird; er stellt den Wert dar, um den sich bei einer Änderung verändert. .

Offensichtlich besteht unsere Aufgabe im Beispiel darin, in der Gleichung solche Koeffizienten zu finden und , bei denen die Abweichungen unserer berechneten Umsatzwerte pro Monat von den tatsächlichen Antworten, d. h. den Werten, die in der Stichprobe präsentiert werden, minimal sind.

Methode der kleinsten Quadrate

Gemäß der Methode der kleinsten Quadrate muss die Abweichung berechnet werden, indem sie quadriert wird. Diese Vorgehensweise vermeidet die wechselseitige Tilgung von Abweichungen, falls sie entgegengesetzte Vorzeichen haben. Zum Beispiel, wenn in einem Fall die Abweichung +5 (plus fünf) beträgt, und im anderen -5 (minus fünf), dann wird die Summe der Abweichungen sich gegenseitig aufheben und 0 (null) betragen. Man kann die Abweichung auch nicht quadrieren, sondern die Eigenschaft des Betrags nutzen, und dann werden alle Abweichungen positiv sein und sich anhäufen. Wir werden nicht näher auf diesen Punkt eingehen, sondern lediglich anmerken, dass es zur Vereinfachung der Berechnungen üblich ist, die Abweichung zu quadrieren.

So sieht die Formel aus, mit der wir die kleinste Summe der Quadrate der Abweichungen (Fehler) bestimmen werden:

wo — das ist die Funktion der Approximation der tatsächlichen Antworten (d. h. unser berechneter Umsatz),

— das sind die tatsächlichen Antworten (der in der Stichprobe bereitgestellte Umsatz),

— das ist der Stichprobenindex (die Monat-Nummer, in der die Abweichung bestimmt wird)

Wir werden die Funktion differenzieren, die Gleichungen der partiellen Ableitungen bestimmen und sind bereit, zur analytischen Lösung überzugehen. Aber zunächst wollen wir einen kleinen Exkurs darüber machen, was Differenzierung ist und die geometrische Bedeutung der Ableitung in Erinnerung rufen.

Differenzierung

Differenzierung ist die Operation zur Bestimmung der Ableitung einer Funktion.

Wozu dient die Ableitung? Die Ableitung einer Funktion charakterisiert die Änderungsrate der Funktion und zeigt uns ihre Richtung an. Wenn die Ableitung an einem bestimmten Punkt positiv ist, nimmt die Funktion zu, andernfalls nimmt die Funktion ab. Je größer der Betrag der Ableitung ist, desto höher ist die Änderungsrate der Funktionswerte und desto steiler ist der Winkel der Funktion.

Zum Beispiel bedeutet in einem kartesischen Koordinatensystem der Wert der Ableitung an dem Punkt M(0,0) +25 , dass an dem gegebenen Punkt bei einer Verschiebung des Wertes nach rechts um eine Einheit der Wert steigt um 25 Einheiten. Auf dem Diagramm sieht das aus wie ein recht steiler Anstieg der Werte von einem bestimmten Punkt.

Ein weiteres Beispiel. Der Wert der Ableitung ist gleich -0,1 bedeutet, dass sich bei einer Verschiebung um eine Einheit der Wert nur um 0,1 Einheit verringert. Dabei können wir im Graphen der Funktion eine kaum erkennbare Neigung nach unten beobachten. Wenn wir eine Analogie zu einem Berg ziehen, steigen wir sehr langsam einen sanften Hang des Berges hinab, im Gegensatz zum vorherigen Beispiel, wo wir sehr steile Gipfel erklimmen mussten :)

Somit, indem wir die Funktion differenzieren nach den Koeffizienten und , definieren wir die Gleichungen der partiellen Ableitungen 1. Ordnung. Nach der Bestimmung der Gleichungen erhalten wir ein System aus zwei Gleichungen, dessen Lösung es uns ermöglicht, solche Koeffizientenwerte zu finden, und , bei denen die Werte der entsprechenden Ableitungen an den gegebenen Punkten sich nur um eine sehr kleine Menge ändern, und im Falle einer analytischen Lösung überhaupt nicht ändern. Mit anderen Worten, die Fehlerfunktion erreicht für die gefundenen Koeffizienten ihr Minimum, da die Werte der partiellen Ableitungen an diesen Punkten gleich null sind.

Also, gemäß den Differenzierungsregeln ist die Gleichung der partiellen Ableitung 1. Ordnung nach dem Koeffizienten nimmt die Form an:

die Gleichung der partiellen Ableitung 1. Ordnung nach nimmt die Form an:

Insgesamt haben wir ein Gleichungssystem erhalten, das eine ziemlich einfache analytische Lösung hat:

begin{equation*}
begin{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
end{cases}
end{equation*}

Bevor wir die Gleichung lösen, laden wir zuerst die Daten hoch, überprüfen die Richtigkeit des Uploads und formatieren die Daten.

Daten hochladen und formatieren

Es ist wichtig zu beachten, dass wir für die analytische Lösung und später für den Gradienten- und stochastischen Gradientenabstieg den Code in zwei Varianten verwenden werden: mit und ohne die Bibliothek, weshalb wir eine entsprechende Datenformatierung benötigen (siehe Code). NumPy Code zum Hochladen und Verarbeiten von Daten

Jetzt, nachdem wir erstens die Daten hochgeladen, zweitens die Richtigkeit des Uploads überprüft und schließlich die Daten formatiert haben, werden wir die erste Visualisierung durchführen. Oft wird dafür die Methode

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

Visualisierung

pairplot Seaborn Bibliothek Seaborn. In unserem Beispiel macht es aufgrund der begrenzten Ziffern keinen Sinn, eine Bibliothek zu verwenden. Seaborn. Wir verwenden eine ganz normale Bibliothek. Matplotlib und schauen uns nur das Streudiagramm an.

Der Code für das Streudiagramm.

print 'Grafik Nr. 1 "Abhängigkeit des Umsatzes vom Monat des Jahres"'

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

Grafik Nr. 1 «Abhängigkeit des Umsatzes vom Monat des Jahres»

Analytische Lösung

Wir verwenden die ganz gewöhnlichen Werkzeuge in python und lösen das Gleichungssystem:

begin{equation*}
begin{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
end{cases}
end{equation*}

Nach Cramers Regel bestimmen wir die allgemeine Determinante sowie die Determinanten nach und nach , danach teilen wir die Determinante nach durch die allgemeine Determinante — und finden den Koeffizienten , analog finden wir den Koeffizienten .

Der Code für die analytische Lösung.

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

Das ist, was wir erhalten haben:

Somit sind die Werte der Koeffizienten gefunden, die Summe der Abweichungsquadrate festgestellt. Lassen Sie uns in das Streudiagramm eine gerade Linie gemäß den gefundenen Koeffizienten zeichnen.

Der Code für die Regressionslinie.

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

Grafik Nr. 2 «Richtige und berechnete Antworten»

Man kann sich die Grafik der Abweichungen für jeden Monat ansehen. In unserem Fall werden wir daraus keinen bedeutenden praktischen Wert ziehen, aber wir befriedigen die Neugier darüber, wie gut die Gleichung der einfachen linearen Regression die Abhängigkeit des Umsatzes vom Monat des Jahres beschreibt.

Der Code für die Grafik der Abweichungen.

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

Grafik Nr. 3 «Abweichungen, %»

Nicht perfekt, aber unsere Aufgabe haben wir erfüllt.

Lassen Sie uns eine Funktion schreiben, die zur Bestimmung der Koeffizienten und die Bibliothek verwendet. NumPy, genauer gesagt — wir schreiben zwei Funktionen: eine mit der Verwendung der Pseudoinversen Matrix (nicht empfohlen in der Praxis, da der Prozess rechenintensiv und instabil ist), die andere mit der Verwendung der Matrizen-Gleichung.

Der Code für die analytische Lösung (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

Wir vergleichen die Zeit, die für die Bestimmung der Koeffizienten verwendet wurde und , gemäß den 3 präsentierten Methoden.

Der Code zur Berechnung der Berechnungszeiten.

print ' 33[1m' + ' 33[4m' + "Die Laufzeit zur Berechnung der Koeffizienten ohne Verwendung der NumPy-Bibliothek:" + ' 33[0m'
% timeit ab_us = Kramer_methode(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Die Laufzeit zur Berechnung der Koeffizienten unter Verwendung der pseudo-inversen Matrix:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Die Laufzeit zur Berechnung der Koeffizienten mit Hilfe der Matrixgleichung:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Bei einer geringen Datenmenge ist die "selbstgeschriebene" Funktion, die Koeffizienten nach der Kramerschen Methode findet, überlegen.

Jetzt können wir zu anderen Methoden der Koeffizientenfindung übergehen. und .

Gradientenabstieg

Zunächst bestimmen wir, was ein Gradient ist. Einfach gesagt, ist der Gradient ein Vektor, der die Richtung des maximalen Anstiegs einer Funktion anzeigt. Vergleichbar mit dem Aufstieg auf einen Berg: Dorthin, wo der Gradient zeigt, gibt es den steilsten Anstieg zur Spitze des Berges. Wenn wir das Beispiel mit dem Berg fortsetzen, müssen wir uns daran erinnern, dass wir eigentlich den steilsten Abstieg benötigen, um schnellstmöglich das Tal zu erreichen, das heißt das Minimum — der Punkt, an dem die Funktion weder steigt noch fällt. An diesem Punkt ist die ableitende Funktion gleich null. Folglich benötigen wir nicht den Gradient, sondern den Antigradienten. Um den Antigradienten zu finden, müssen wir lediglich den Gradient mit -1 (minus eins) multiplizieren.

Beachten wir, dass die Funktion mehrere Minima haben kann, und wenn wir nach dem folgenden Algorithmus in eines von ihnen gefallen sind, können wir möglicherweise ein anderes Minimum nicht finden, das tiefer liegt als das gefundene. Entspannt euch, das passiert uns nicht! In unserem Fall haben wir es mit einem einzigen Minimum zu tun, da unsere Funktion im Diagramm eine gewöhnliche Parabel darstellt. Und wie wir alle aus dem Mathematikunterricht wissen sollten — eine Parabel hat nur ein Minimum.

Nachdem wir geklärt haben, warum wir den Gradient benötigen und dass der Gradient ein Vektor mit bestimmten Koordinaten ist, der genau die gesuchten Koeffizienten darstellt, und können wir den Gradientenabstieg implementieren.

Vor dem Start schlage ich vor, buchstäblich ein paar Sätze über den Abwärtsalgorithmus zu lesen:

• Wir bestimmen die Koeffizientenkoordinaten pseudorandom. und . In unserem Beispiel werden wir die Koeffizienten in der Nähe von Null bestimmen. Dies ist eine gängige Praxis, jedoch kann es für jeden Fall eine eigene Vorgehensweise geben.
• Von der Koordinate ziehen wir den Wert der partiellen Ableitung erster Ordnung an dem Punkt ab . So, wenn die Ableitung positiv ist, steigt die Funktion. Folglich bewegen wir uns, indem wir den Wert der Ableitung abziehen, in die entgegengesetzte Richtung des Anstiegs, also in Richtung des Abfalls. Ist die Ableitung negativ, dann sinkt die Funktion an diesem Punkt, und indem wir den Wert der Ableitung abziehen, bewegen wir uns in Richtung des Abfalls.
• Wir führen eine ähnliche Operation mit der Koordinate durch: wir ziehen den Wert der partiellen Ableitung an dem Punkt ab .
• Um nicht das Minimum zu überspringen und in den Tiefen des Weltraums zu verschwinden, ist es notwendig, die Schrittgröße in Richtung des Abfalls festzulegen. Insgesamt könnte man einen ganzen Artikel darüber schreiben, wie man die Schrittgröße korrekt festlegt und wie man sie während des Abstiegs ändert, um die Rechenkosten zu senken. Doch jetzt stehen wir vor einer etwas anderen Aufgabe und wir werden durch die wissenschaftliche Methode des „Tastens“ oder wie man umgangssprachlich sagt, auf empirische Weise, die Schrittgröße festlegen.
• Nachdem wir von den vorgegebenen Koordinaten und die Werte der Ableitungen abgezogen haben, erhalten wir neue Koordinaten und . Wir machen den nächsten Schritt (Abzug) bereits von den berechneten Koordinaten. Und so wird der Zyklus immer wieder neu gestartet, bis die erforderliche Konvergenz erreicht ist.

Alles! Jetzt sind wir bereit, auf die Suche nach dem tiefsten Punkt des Marianengrabens zu gehen. Fangen wir an.

Code für den Gradientenabstieg

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

Wir sind auf dem Grund des Marianengrabens angekommen und haben dort all die gleichen Koeffizienten gefunden und , was eigentlich zu erwarten war.

Wir werden einen weiteren Tauchgang machen, aber dieses Mal wird die Ausrüstung unseres Tiefseeapparates andere Technologien umfassen, nämlich die Bibliothek NumPy.

Code für den Gradientenabstieg (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

Die Werte der Koeffizienten und bleiben unverändert.

Lassen Sie uns ansehen, wie sich der Fehler während des Gradientenabstiegs änderte, d.h. wie sich die Summe der quadrierten Abweichungen mit jedem Schritt veränderte.

Code für das Diagramm der Summen der quadrierten Abweichungen

print 'Diagramm №4 "Summe der quadrierten Abweichungen schrittweise"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Schritte (Iteration)', size=16)
plt.ylabel('Summe der quadrierten Abweichungen', size=16)
plt.show()

Grafik Nr. 4 „Summe der Quadrate der Abweichungen beim Gradientenabstieg“

Im Diagramm sehen wir, dass der Fehler mit jedem Schritt abnimmt und nach einer bestimmten Anzahl von Iterationen eine fast horizontale Linie beobachtet werden kann.

Zuletzt bewerten wir den Unterschied in der Ausführungszeit des Codes:

Code zur Bestimmung der Berechnungszeit des Gradientenabstiegs

print ' 33[1m' + ' 33[4m' + "Die Ausführungszeit des Gradientenabstiegs ohne die Verwendung der NumPy-Bibliothek:" + ' 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' + "Die Ausführungszeit des Gradientenabstiegs mit der Verwendung der NumPy-Bibliothek:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Vielleicht machen wir etwas falsch, aber wieder einmal übertrifft eine einfache „selbstgeschriebene“ Funktion, die die Bibliothek nicht verwendet, im Hinblick auf die Ausführungszeit die Funktion, die die Bibliothek nutzt. NumPy Wir stehen jedoch nicht still, sondern bewegen uns in Richtung des Studiums einer weiteren faszinierenden Methode zur Lösung der Gleichung der einfachen linearen Regression. Herzlich willkommen! NumPy.

Um das Prinzip des stochastischen Gradientenabstiegs schneller zu verstehen, ist es besser, seine Unterschiede zum gewöhnlichen Gradientenabstieg zu klären. Wir verwendeten im Fall des Gradientenabstiegs in den Ableitungsformeln

Stochastischer Gradientenabstieg

die Summen der Werte aller Merkmale und der wahren Antworten, die in der Stichprobe vorhanden sind (also die Summen aller und ). Beim stochastischen Gradientenabstieg verwenden wir nicht alle Werte in der Stichprobe, sondern wählen stattdessen zufällig einen sogenannten Stichprobenindex aus und verwenden dessen Werte. und Wenn beispielsweise der Index die Nummer 3 (drei) ergibt, nehmen wir die Werte

, setzen die Werte in die Ableitungsformeln ein und bestimmen neue Koordinaten. Dann, nachdem wir die Koordinaten bestimmt haben, wählen wir erneut zufällig einen Stichprobenindex aus, setzen die entsprechenden Werte in die Gleichungen der partiellen Ableitungen ein und bestimmen erneut die Koordinaten und usw. bis zur Konvergenz. Auf den ersten Blick könnte es scheinen, dass dies überhaupt funktionieren kann; allerdings funktioniert es. Es ist jedoch anzumerken, dass der Fehler nicht mit jedem Schritt abnimmt, aber der Trend ist eindeutig vorhanden. und und so weiter bis zur Konvergenz. Auf den ersten Blick mag es scheinen, als könnte das überhaupt funktionieren, aber es funktioniert. Es ist zwar zu beachten, dass sich der Fehler nicht mit jedem Schritt verringert, aber die Tendenz ist ohne Zweifel vorhanden.

Was sind die Vorteile des stochastischen Gradientenabstiegs im Vergleich zum üblichen Verfahren? Wenn wir eine sehr große Stichprobe haben, die aus Zehntausenden von Werten besteht, ist es viel einfacher, beispielsweise tausend zufällige Werte daraus zu verarbeiten, anstatt die gesamte Stichprobe zu verwenden. In diesem Fall wird der stochastische Gradientabstieg aktiviert. In unserem Fall werden wir natürlich keinen großen Unterschied feststellen.

Schauen wir uns den Code an.

Code für den stochastischen Gradientenabstieg

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

Betrachten wir die Koeffizienten genau und fragen uns: „Wie kann das sein?“. Wir haben andere Koeffizientenwerte erhalten. und Vielleicht hat der stochastische Gradientabstieg optimalere Parameter für die Gleichung gefunden? Leider nein. Es genügt, die Summe der quadratischen Abweichungen zu betrachten und zu sehen, dass die Fehler bei den neuen Koeffizienten höher geworden sind. Lassen Sie uns nicht verzweifeln. Wir werden ein Diagramm zeichnen, um die Fehlerentwicklung darzustellen.

Code für das Diagramm der Summe der quadratischen Abweichungen beim stochastischen Gradientenabstieg

print 'Diagramm Nr. 5 "Summe der quadratischen Abweichungen schrittweise"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Schritte (Iteration)', size=16)
plt.ylabel('Summe der quadratischen Abweichungen', size=16)
plt.show()

Diagramm Nr. 5 „Summe der quadratischen Abweichungen beim stochastischen Gradientenabstieg“

Nachdem wir das Diagramm betrachtet haben, wird alles klar und wir werden jetzt alles korrigieren.

Was ist also passiert? Folgendes ist passiert. Wenn wir zufällig einen Monat auswählen, dann versucht unser Algorithmus, den Fehler bei der Berechnung der Einnahmen für genau diesen Monat zu minimieren. Dann wählen wir einen anderen Monat und wiederholen die Berechnung, legen jedoch den Fokus auf die Minimierung des Fehlers für den zweiten gewählten Monat. Und jetzt erinnern wir uns, dass die ersten beiden Monate erheblich von der Geradengleichung der einfachen linearen Regression abweichen. Das bedeutet, dass, wenn einer dieser beiden Monate gewählt wird, das Minimieren des Fehlers für jeden von ihnen den Fehler für die gesamte Stichprobe erheblich erhöht. Was sollten wir also tun? Die Antwort ist einfach: Wir müssen die Lernrate verringern. Denn wenn wir die Schritte verringern, wird der Fehler auch nicht mehr so stark „springen“, also nicht mehr so heftig nach oben und unten schwanken. Genauer gesagt, der Fehler wird weiterhin springen, aber nicht mehr so schnell :) Lassen Sie uns das überprüfen.

Code für den SGD-Start mit kleinerem Schritt

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

Diagramm Nr. 6 „Summe der quadratischen Abweichungen beim stochastischen Gradientenabstieg (80.000 Schritte)“

Die Werte der Koeffizienten haben sich verbessert, sind aber immer noch nicht ideal. Hypothetisch könnten wir das folgendermaßen beheben. Wir wählen zum Beispiel die Werte der Koeffizienten aus den letzten 1000 Iterationen aus, bei denen der minimale Fehler aufgetreten ist. Allerdings müssen wir dazu auch die tatsächlichen Werte der Koeffizienten aufzeichnen. Wir werden das nicht tun, sondern uns lieber auf das Diagramm konzentrieren. Es sieht glatt aus, und der Fehler scheint gleichmäßig zu sinken. In Wirklichkeit ist das jedoch nicht der Fall. Schauen wir uns die ersten 1000 Iterationen an und vergleichen wir sie mit den letzten.

Code für das Diagramm SGD (erste 1000 Schritte)

print 'Diagramm Nr. 7 "Summe der quadrierten Abweichungen Schritt für Schritt. Erste 1000 Iterationen"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Schritte (Iteration)', size=16)
plt.ylabel('Summe der quadrierten Abweichungen', size=16)
plt.show()

print 'Diagramm Nr. 7 "Summe der quadrierten Abweichungen Schritt für Schritt. Letzte 1000 Iterationen"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Schritte (Iteration)', size=16)
plt.ylabel('Summe der quadrierten Abweichungen', size=16)
plt.show()

Diagramm Nr. 7 „Summe der quadrierten Abweichungen SGD (erste 1000 Schritte)“

Diagramm Nr. 8 „Summe der quadrierten Abweichungen SGD (letzte 1000 Schritte)“

Zu Beginn des Abstiegs beobachten wir eine relativ gleichmäßige und steile Abnahme des Fehlers. In den letzten Iterationen sehen wir, dass der Fehler um den Wert von 1,475 schwankt und in einigen Momenten sogar diesem optimalen Wert entspricht, aber dann dennoch wieder ansteigt… Ich wiederhole, man könnte die Werte der Koeffizienten aufzeichnen und , und dann die auswählen, bei denen der Fehler minimal ist. Allerdings haben wir ein ernsthafteres Problem: Wir mussten 80.000 Schritte (siehe Code) machen, um Werte zu erhalten, die nahe an den optimalen liegen. Das widerspricht bereits der Idee, die Rechenzeit beim stochastischen Gradientenabstieg im Vergleich zum Gradientenabstieg zu sparen. Was könnte verbessert und behoben werden? Es ist nicht schwer zu erkennen, dass wir in den ersten Iterationen stetig abwärts gehen und daher sollten wir einen großen Schritt in den ersten Iterationen beibehalten und mit fortschreitendem Fortschritt die Schrittweite verringern. Das werden wir in diesem Artikel jedoch nicht machen — er hat sich bereits lange genug hingezogen. Interessierte können selbst darüber nachdenken, wie man das umsetzen kann; es ist nicht kompliziert 🙂

Jetzt führen wir den stochastischen Gradientenabstieg unter Verwendung der Bibliothek durch NumPy (und wir werden nicht über die Steine stolpern, die wir zuvor identifiziert haben)

Code für den stochastischen Gradientenabstieg (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

Die Werte sind fast die gleichen wie beim Abstieg ohne Verwendung NumPy. Das ist allerdings logisch.

Lass uns herausfinden, wie viel Zeit uns die stochastischen Gradientenabstiege gekostet haben.

Code zur Bestimmung der Berechnungszeit von SGD (80.000 Schritte)

print ' 33[1m' + ' 33[4m' +
"Die Ausführungszeit des stochastischen Gradientenabstiegs ohne Verwendung der NumPy-Bibliothek:"
+ ' 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' +
"Die Ausführungszeit des stochastischen Gradientenabstiegs mit der Verwendung der NumPy-Bibliothek:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

Je weiter man in den Wald geht, desto dunkler werden die Wolken: Wieder zeigt die „hausgemachte“ Formel das bessere Ergebnis. All dies lässt darauf schließen, dass es noch subtilere Methoden zur Nutzung der Bibliothek geben muss NumPy, die tatsächlich die Berechnungsoperationen beschleunigen. In diesem Artikel werden wir darüber nicht mehr erfahren. Es wird etwas sein, worüber man in der Freizeit nachdenken kann :)

Lassen Sie uns zusammenfassen

Bevor wir zusammenfassen, möchten wir eine Frage beantworten, die wahrscheinlich bei unserem geschätzten Leser aufgetaucht ist. Warum also solche „Qualen“ mit den Abstiegen, warum müssen wir den Berg hinauf und hinunter (meist hinunter) gehen, um die begehrte Senke zu finden, wenn wir ein so mächtiges und einfaches Gerät in Form einer analytischen Lösung in den Händen halten, das uns sofort an den gewünschten Ort teleportiert?

Die Antwort auf diese Frage liegt auf der Hand. Wir haben jetzt ein sehr einfaches Beispiel behandelt, in dem die wahre Antwort von einem Merkmal abhängt . Solch eine Situation begegnet man im Leben nicht oft, also stellen wir uns vor, wir haben 2, 30, 50 oder mehr Merkmale. Fügen wir dazu Tausende oder sogar Zehntausende von Werten für jedes Merkmal hinzu. In diesem Fall könnte die analytische Lösung eine Überlastung nicht standhalten und versagen. Im Gegenzug wird der Gradientenabstieg und seine Variationen uns langsam, aber sicher dem Ziel – dem Minimum der Funktion – näherbringen. Und bezüglich der Geschwindigkeit brauchen Sie sich keine Sorgen zu machen – wir werden sicherlich noch Methoden durchgehen, die es uns ermöglichen, die Schrittgröße (also die Geschwindigkeit) zu bestimmen und zu regulieren.

Und jetzt zur eigentlichen kurzen Zusammenfassung.

Zuerst hoffe ich, dass das in dem Artikel behandelte Material den angehenden „Data Scientists“ hilft, zu verstehen, wie man die Gleichungen der einfachen (und nicht nur) linearen Regression löst.

Zweitens haben wir mehrere Ansätze zur Lösung der Gleichung untersucht. Abhängig von der Situation können wir jetzt diejenige auswählen, die am besten geeignet ist, um das gegebene Problem zu lösen.

Drittens haben wir die Kraft zusätzlicher Anpassungen, insbesondere der Schrittgröße des Gradientenabstiegs, gesehen. Diese Parameter sollte man nicht vernachlässigen. Wie oben erwähnt, sollte man zur Reduzierung der Berechnungskosten die Schrittgröße während des Abstiegs anpassen.

Zudem haben in unserem Fall die „selbstgeschriebenen“ Funktionen die besten Berechnungszeiten gezeigt. Wahrscheinlich liegt dies an einer nicht ganz professionellen Nutzung der Möglichkeiten der Bibliothek. NumPy. Doch wie dem auch sei, die folgende Schlussfolgerung drängt sich auf. Einerseits ist es manchmal sinnvoll, etablierte Meinungen in Frage zu stellen, andererseits sollte man nicht immer alles komplizierter machen – manchmal ist eine einfachere Herangehensweise effektiver. Und da unser Ziel war, drei Ansätze zur Lösung der Gleichung der einfachen linearen Regression zu analysieren, hat uns die Verwendung von „selbstgeschriebenen“ Funktionen völlig ausgereicht.

Literatur (oder etwas Ähnliches)

1. Lineare Regression

http://statistica.ru/theory/osnovy-lineynoy-regressii/

2. Methode der kleinsten Quadrate

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Ableitung

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Gradientabstieg

habr.com/ru/post/471458

habr.com/ru/post/307312

artemarakcheev.com/2017-12-31/linear_regression

6. NumPy Bibliothek

docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.linalg.solve.html

docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.pinv.html

pythonworld.ru/numpy/2.html

Quelle: habr.com