🥇Die Gleichung der einfachen linearen Regression lösen

Der Artikel diskutiert verschiedene Möglichkeiten, die mathematische Gleichung einer einfachen (gepaarten) Regressionsgeraden zu bestimmen.

Alle hier besprochenen 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 stellt der Artikel verschiedene Funktionen bereit, die hauptsächlich in solche unterteilt werden, die ohne Verwendung der Bibliothek geschrieben werden NumPy und diejenigen, die für Berechnungen verwendet werden NumPy. Es wird angenommen, dass der geschickte Einsatz NumPy wird die Rechenkosten senken.

Der gesamte im Artikel angegebene Code ist in der Sprache geschrieben Python 2.7 mit Jupyter Notizbuch. Der Quellcode und die Datei mit Beispieldaten werden auf veröffentlicht Github

Der Artikel richtet sich eher an Anfänger als auch an diejenigen, die bereits allmählich begonnen haben, das Studium eines sehr breiten Abschnitts der künstlichen Intelligenz – maschinelles Lernen – zu beherrschen.

Zur Veranschaulichung des Materials verwenden wir ein sehr einfaches Beispiel.

Beispielbedingungen

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

Tabelle Nr. 1 „Beispielbedingungen“

Wir gehen davon aus, dass die Werte ist der Monat des Jahres und – Umsatz in diesem Monat. Mit anderen Worten, der Umsatz hängt vom Monat des Jahres ab und - das einzige Zeichen, von dem der Umsatz abhängt.

Das Beispiel ist mittelmäßig, sowohl im Hinblick auf die bedingte Abhängigkeit des Umsatzes vom Monat des Jahres als auch im Hinblick auf die Anzahl der Werte – es gibt nur sehr wenige davon. Eine solche Vereinfachung wird es jedoch, wie man so sagt, ermöglichen, den Stoff, den sich Anfänger aneignen, nicht immer leicht zu erklären. Und auch die Einfachheit der Zahlen ermöglicht es denjenigen, die das Beispiel auf Papier ohne nennenswerten Arbeitsaufwand lösen möchten.

Nehmen wir an, dass die im Beispiel angegebene Abhängigkeit recht gut durch die mathematische Gleichung einer einfachen (paarigen) Regressionsgeraden der Form angenähert werden kann:

wo ist der Monat, in dem die Einnahmen eingegangen sind, — Einnahmen entsprechend dem Monat, и sind die Regressionskoeffizienten der geschätzten Linie.

Beachten Sie, dass der Koeffizient oft als Steigung oder Gefälle der geschätzten Linie bezeichnet; stellt den Betrag dar, um den die wenn es sich ändert .

Offensichtlich besteht unsere Aufgabe im Beispiel darin, solche Koeffizienten in der Gleichung auszuwählen и , bei dem die Abweichungen unserer berechneten Umsatzwerte pro Monat von den wahren Antworten, d.h. Die in der Probe dargestellten Werte sind minimal.

Methode der kleinsten Quadrate

Nach der Methode der kleinsten Quadrate sollte die Abweichung durch Quadrieren berechnet werden. Mit dieser Technik können Sie die gegenseitige Aufhebung von Abweichungen vermeiden, wenn diese entgegengesetzte Vorzeichen haben. Wenn beispielsweise in einem Fall die Abweichung beträgt +5 (plus fünf) und im anderen -5 (minus fünf), dann heben sich die Abweichungen in der Summe auf und betragen 0 (Null). Es ist möglich, die Abweichung nicht zu quadrieren, sondern die Eigenschaft des Moduls zu nutzen, und dann sind alle Abweichungen positiv und summieren sich. Wir werden auf diesen Punkt nicht näher eingehen, sondern lediglich darauf hinweisen, dass es zur Vereinfachung der Berechnungen üblich ist, die Abweichung zu quadrieren.

So sieht die Formel aus, mit der wir die kleinste Summe der quadratischen Abweichungen (Fehler) ermitteln:

wo ist eine Funktion der Annäherung an wahre Antworten (d. h. der von uns berechnete Umsatz),

sind die wahren Antworten (in der Stichprobe angegebener Umsatz),

ist der Stichprobenindex (Nummer des Monats, in dem die Abweichung ermittelt wird)

Lassen Sie uns die Funktion differenzieren, die partiellen Differentialgleichungen definieren und bereit sein, mit der analytischen Lösung fortzufahren. Aber lassen Sie uns zunächst einen kurzen Exkurs darüber machen, was Differentiation ist, und uns an die geometrische Bedeutung der Ableitung erinnern.

Differenzierung

Differenzierung ist die Operation, die Ableitung einer Funktion zu finden.

Wofür wird die Ableitung verwendet? Die Ableitung einer Funktion charakterisiert die Änderungsrate der Funktion und gibt uns ihre Richtung an. Wenn die Ableitung an einem bestimmten Punkt positiv ist, nimmt die Funktion zu, andernfalls nimmt sie ab. Und je größer der Wert der absoluten Ableitung ist, desto höher ist die Änderungsrate der Funktionswerte und desto steiler ist die Steigung des Funktionsgraphen.

Unter den Bedingungen eines kartesischen Koordinatensystems ist beispielsweise der Wert der Ableitung am Punkt M(0,0) gleich +25 bedeutet, dass an einem bestimmten Punkt der Wert verschoben wird nach rechts um eine konventionelle Einheit, den Wert erhöht sich um 25 konventionelle Einheiten. In der Grafik sieht es nach einem ziemlich steilen Anstieg der Werte aus von einem bestimmten Punkt aus.

Ein anderes Beispiel. Der Ableitungswert ist gleich -0,1 bedeutet das, wenn es verschoben wird pro konventioneller Einheit, Wert verringert sich lediglich um 0,1 herkömmliche Einheiten. Gleichzeitig können wir im Diagramm der Funktion eine kaum wahrnehmbare Abwärtsneigung beobachten. Um eine Analogie zu einem Berg zu ziehen: Es ist, als würden wir sehr langsam einen sanften Hang von einem Berg hinabsteigen, anders als im vorherigen Beispiel, wo wir sehr steile Gipfel erklimmen mussten :)

Also nach Differenzieren der Funktion durch Chancen и , definieren wir partielle Differentialgleichungen 1. Ordnung. Nach der Bestimmung der Gleichungen erhalten wir ein System aus zwei Gleichungen, durch deren Lösung wir solche Werte der Koeffizienten auswählen können и , bei dem sich die Werte der entsprechenden Ableitungen an bestimmten Punkten um einen sehr, sehr geringen Betrag ändern und sich im Falle einer analytischen Lösung überhaupt nicht ändern. Mit anderen Worten, die Fehlerfunktion an den gefundenen Koeffizienten erreicht ein Minimum, da die Werte der partiellen Ableitungen an diesen Punkten gleich Null sind.

Also nach den Regeln der Differentiation die partielle Ableitungsgleichung 1. Ordnung nach dem Koeffizienten wird die Form annehmen:

Partielle Ableitungsgleichung 1. Ordnung nach wird die Form annehmen:

Als Ergebnis erhielten wir ein Gleichungssystem, das eine ziemlich einfache analytische Lösung hat:

begin{Gleichung*}
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{Gleichung*}

Bevor wir die Gleichung lösen, laden wir sie vor, prüfen, ob der Ladevorgang korrekt ist, und formatieren die Daten.

Laden und Formatieren von Daten

Es ist zu beachten, dass wir für die analytische Lösung und anschließend für den Gradienten und den stochastischen Gradientenabstieg den Code in zwei Varianten verwenden werden: unter Verwendung der Bibliothek NumPy und ohne es zu verwenden, benötigen wir eine entsprechende Datenformatierung (siehe Code).

Code zum Laden und Verarbeiten von Daten

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

Nachdem wir nun erstens die Daten geladen, zweitens die Richtigkeit des Ladevorgangs überprüft und schließlich die Daten formatiert haben, führen wir die erste Visualisierung durch. Die hierfür häufig verwendete Methode ist Paarplot Bibliothek Seegeboren. In unserem Beispiel macht es aufgrund der begrenzten Anzahl keinen Sinn, die Bibliothek zu nutzen Seegeboren. Wir werden die reguläre Bibliothek nutzen Matplotlib und schauen Sie sich einfach das Streudiagramm an.

Streudiagrammcode

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

Diagramm Nr. 1 „Abhängigkeit der Einnahmen vom Monat des Jahres“

Analytische Lösung

Lassen Sie uns die gängigsten Tools verwenden python und löse das Gleichungssystem:

begin{Gleichung*}
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{Gleichung*}

Nach der Cramer-Regel Wir werden die allgemeine Determinante sowie Determinanten finden und , danach Division der Determinante durch zur allgemeinen Determinante - finden Sie den Koeffizienten , ähnlich finden wir den Koeffizienten .

Analytischer Lösungscode

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

Damit sind die Werte der Koeffizienten gefunden und die Summe der quadrierten Abweichungen ermittelt. Zeichnen wir im Streuhistogramm eine gerade Linie entsprechend den gefundenen Koeffizienten.

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

Diagramm Nr. 2 „Richtige und berechnete Antworten“

Sie können sich das Abweichungsdiagramm für jeden Monat ansehen. In unserem Fall werden wir daraus keinen nennenswerten praktischen Wert ziehen, aber wir werden unsere Neugier befriedigen, wie gut die einfache lineare Regressionsgleichung die Abhängigkeit des Umsatzes vom Monat des Jahres charakterisiert.

Abweichungsdiagrammcode

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

Diagramm Nr. 3 „Abweichungen, %“

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

Schreiben wir eine Funktion, um die Koeffizienten zu bestimmen и nutzt die Bibliothek NumPyGenauer gesagt werden wir zwei Funktionen schreiben: eine unter Verwendung einer pseudoinversen Matrix (in der Praxis nicht empfohlen, da der Prozess rechentechnisch komplex und instabil ist), die andere unter Verwendung einer Matrixgleichung.

Analytischer Lösungscode (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

Vergleichen wir die Zeit, die für die Bestimmung der Koeffizienten aufgewendet wurde и , gemäß den 3 vorgestellten Methoden.

Code zur Berechnung der Berechnungszeit

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)

Bei einer kleinen Datenmenge hat eine „selbstgeschriebene“ Funktion die Nase vorn, die die Koeffizienten mithilfe der Cramer-Methode ermittelt.

Jetzt können Sie mit anderen Methoden zum Ermitteln von Koeffizienten fortfahren и .

Gradientenabstieg

Definieren wir zunächst, was ein Farbverlauf ist. Einfach ausgedrückt ist der Gradient ein Segment, das die Richtung des maximalen Wachstums einer Funktion angibt. In Analogie zum Besteigen eines Berges ist der Aufstieg zum Gipfel dort, wo das Gefälle herrscht, am steilsten. Wenn wir das Beispiel mit dem Berg entwickeln, erinnern wir uns daran, dass wir tatsächlich den steilsten Abstieg brauchen, um so schnell wie möglich das Tiefland zu erreichen, also das Minimum – den Ort, an dem die Funktion nicht zu- oder abnimmt. Zu diesem Zeitpunkt ist die Ableitung gleich Null. Daher benötigen wir keinen Gradienten, sondern einen Antigradienten. Um den Antigradienten zu finden, müssen Sie nur den Gradienten mit multiplizieren -1 (minus eins).

Achten wir auf die Tatsache, dass eine Funktion mehrere Minima haben kann, und nachdem wir mit dem unten vorgeschlagenen Algorithmus in eines davon abgetaucht sind, können wir kein anderes Minimum finden, das möglicherweise niedriger als das gefundene ist. Entspannen wir uns, das ist keine Bedrohung für uns! In unserem Fall haben wir es mit einem einzigen Minimum zu tun, da unsere Funktion in der Grafik ist eine regelmäßige Parabel. Und wie wir alle aus unserem Schulmathematikkurs ganz genau wissen dürften, hat eine Parabel nur ein Minimum.

Nachdem wir herausgefunden hatten, warum wir einen Gradienten brauchten, und auch, dass der Gradient ein Segment ist, also ein Vektor mit gegebenen Koordinaten, die genau die gleichen Koeffizienten sind и Wir können einen Gradientenabstieg implementieren.

Bevor ich anfange, empfehle ich, nur ein paar Sätze über den Abstiegsalgorithmus zu lesen:

Wir bestimmen pseudozufällig die Koordinaten der Koeffizienten и . In unserem Beispiel definieren wir Koeffizienten nahe Null. Dies ist eine gängige Praxis, aber jeder Fall kann seine eigene Praxis haben.
Von der Koordinate Subtrahieren Sie den Wert der partiellen Ableitung 1. Ordnung an diesem Punkt . Wenn also die Ableitung positiv ist, nimmt die Funktion zu. Wenn wir also den Wert der Ableitung subtrahieren, bewegen wir uns in die entgegengesetzte Wachstumsrichtung, also in die Abstiegsrichtung. Wenn die Ableitung negativ ist, nimmt die Funktion an diesem Punkt ab und indem wir den Wert der Ableitung subtrahieren, bewegen wir uns in Abstiegsrichtung.
Eine ähnliche Operation führen wir mit der Koordinate durch : Subtrahieren Sie den Wert der partiellen Ableitung an diesem Punkt .
Um das Minimum nicht zu überschreiten und in den Weltraum zu fliegen, ist es notwendig, die Schrittweite in Abstiegsrichtung einzustellen. Im Allgemeinen könnte man einen ganzen Artikel darüber schreiben, wie man den Schritt richtig einstellt und wie man ihn während des Abstiegs ändert, um den Rechenaufwand zu reduzieren. Aber jetzt haben wir eine etwas andere Aufgabe vor uns und werden die Schrittgröße mit der wissenschaftlichen Methode „Poke“ oder, wie man im Volksmund sagt, empirisch ermitteln.
Sobald wir von den angegebenen Koordinaten entfernt sind и Wenn wir die Werte der Ableitungen subtrahieren, erhalten wir neue Koordinaten и . Wir machen den nächsten Schritt (Subtraktion), bereits von den berechneten Koordinaten. Und so beginnt der Zyklus immer wieder, bis die erforderliche Konvergenz erreicht ist.

Alle! Jetzt sind wir bereit, uns auf die Suche nach der tiefsten Schlucht des Marianengrabens zu machen. Lass uns anfangen.

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 tauchten bis zum Grund des Marianengrabens und fanden dort alle die gleichen Koeffizientenwerte и , was genau das war, was zu erwarten war.

Machen wir noch einen weiteren Tauchgang, nur wird unser Tiefseefahrzeug dieses Mal mit anderen Technologien gefüllt, nämlich einer 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

Koeffizientenwerte и unveränderlich.

Schauen wir uns an, wie sich der Fehler während des Gradientenabstiegs ändert, d. h. wie sich die Summe der quadratischen Abweichungen mit jedem Schritt ändert.

Code zum Zeichnen von Summen quadrierter Abweichungen

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

Grafik Nr. 4 „Summe der quadratischen Abweichungen beim Gefälleabstieg“

In der Grafik sehen wir, dass der Fehler mit jedem Schritt abnimmt und wir nach einer bestimmten Anzahl von Iterationen eine fast horizontale Linie beobachten.

Abschließend schätzen wir den Unterschied in der Codeausführungszeit ab:

Code zur Bestimmung der Berechnungszeit für den Gradientenabstieg

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)

Vielleicht machen wir etwas falsch, aber auch hier handelt es sich um eine einfache „selbst geschriebene“ Funktion, die die Bibliothek nicht verwendet NumPy übertrifft die Berechnungszeit einer Funktion, die die Bibliothek verwendet NumPy.

Aber wir stehen nicht still, sondern untersuchen einen weiteren spannenden Weg zur Lösung der einfachen linearen Regressionsgleichung. Treffen!

Stochastischer Gradientenabstieg

Um das Funktionsprinzip des stochastischen Gradientenabstiegs schnell zu verstehen, ist es besser, seine Unterschiede zum gewöhnlichen Gradientenabstieg zu bestimmen. Wir, im Fall des Gradientenabstiegs, in den Gleichungen der Ableitungen von и verwendete die Summen der Werte aller in der Stichprobe verfügbaren Merkmale und wahren Antworten (d. h. die Summen aller). и ). Beim stochastischen Gradientenabstieg verwenden wir nicht alle in der Stichprobe vorhandenen Werte, sondern wählen stattdessen pseudozufällig den sogenannten Stichprobenindex aus und verwenden seine Werte.

Wenn der Index beispielsweise die Zahl 3 (drei) ist, nehmen wir die Werte и , dann setzen wir die Werte in die Ableitungsgleichungen ein und bestimmen neue Koordinaten. Nachdem wir die Koordinaten bestimmt haben, bestimmen wir erneut pseudozufällig den Stichprobenindex, setzen die dem Index entsprechenden Werte in die partiellen Differentialgleichungen ein und bestimmen die Koordinaten auf neue Weise и usw. bis die Konvergenz grün wird. Auf den ersten Blick sieht es vielleicht nicht so aus, als ob das überhaupt funktionieren könnte, aber es funktioniert. Es ist zwar erwähnenswert, dass der Fehler nicht mit jedem Schritt kleiner wird, aber es gibt durchaus eine Tendenz.

Was sind die Vorteile des stochastischen Gradientenabstiegs gegenüber dem herkömmlichen? Wenn unsere Stichprobe sehr groß ist und Zehntausende von Werten umfasst, ist es viel einfacher, beispielsweise zufällig ausgewählte Tausend davon zu verarbeiten, als die gesamte Stichprobe. Hier kommt der stochastische Gradientenabstieg ins Spiel. In unserem Fall werden wir natürlich keinen großen Unterschied bemerken.

Schauen wir uns den Code an.

Code für 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])

Wir schauen uns die Koeffizienten genau an und stellen uns die Frage: „Wie kann das sein?“ Wir haben andere Koeffizientenwerte erhalten и . Vielleicht hat der stochastische Gradientenabstieg optimalere Parameter für die Gleichung gefunden? Leider nein. Es reicht aus, die Summe der quadratischen Abweichungen zu betrachten und festzustellen, dass der Fehler bei neuen Werten der Koeffizienten größer ist. Wir haben es nicht eilig, zu verzweifeln. Lassen Sie uns ein Diagramm der Fehleränderung erstellen.

Code zum Zeichnen der Summe der quadratischen Abweichungen beim stochastischen Gradientenabfall

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

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

Wenn wir uns den Zeitplan ansehen, passt alles zusammen und jetzt werden wir alles in Ordnung bringen.

Also was ist passiert? Folgendes ist passiert. Wenn wir einen Monat zufällig auswählen, versucht unser Algorithmus, den Fehler bei der Berechnung des Umsatzes für den ausgewählten Monat zu reduzieren. Dann wählen wir einen anderen Monat aus und wiederholen die Berechnung, reduzieren jedoch den Fehler für den zweiten ausgewählten Monat. Bedenken Sie nun, dass die ersten beiden Monate erheblich von der Linie der einfachen linearen Regressionsgleichung abweichen. Das bedeutet, dass unser Algorithmus bei Auswahl eines dieser beiden Monate den Fehler für die gesamte Stichprobe erheblich erhöht, indem er den Fehler jedes dieser Monate verringert. Was also tun? Die Antwort ist einfach: Sie müssen die Abstiegsstufe reduzieren. Denn durch die Reduzierung der Abstiegsstufe hört der Fehler auch auf, auf und ab zu „springen“. Oder besser gesagt, der „springende“ Fehler wird nicht aufhören, aber nicht so schnell :) Schauen wir mal nach.

Code zum Ausführen von SGD mit kleineren Schritten

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

Grafik Nr. 6 „Summe der quadratischen Abweichungen während des stochastischen Gradientenabstiegs (80 Schritte)“

Die Koeffizienten haben sich verbessert, sind aber immer noch nicht ideal. Hypothetisch kann dies auf diese Weise korrigiert werden. Wir wählen beispielsweise in den letzten 1000 Iterationen die Werte der Koeffizienten aus, mit denen der minimale Fehler gemacht wurde. Dafür müssen wir zwar auch die Werte der Koeffizienten selbst aufschreiben. Wir werden dies nicht tun, sondern auf den Zeitplan achten. Es sieht glatt aus und der Fehler scheint gleichmäßig abzunehmen. Eigentlich stimmt das nicht. Schauen wir uns die ersten 1000 Iterationen an und vergleichen sie mit den letzten.

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

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

Grafik Nr. 7 „Summe der quadratischen Abweichungen SGD (erste 1000 Schritte)“

Grafik Nr. 8 „Summe der quadratischen Abweichungen SGD (letzte 1000 Schritte)“

Gleich zu Beginn des Abstiegs beobachten wir einen ziemlich gleichmäßigen und steilen Abfall des Fehlers. In den letzten Iterationen sehen wir, dass der Fehler immer um den Wert von 1,475 herumgeht und in einigen Momenten sogar diesem optimalen Wert entspricht, aber dann steigt er immer noch ... Ich wiederhole, Sie können die Werte aufschreiben Koeffizienten и , und wählen Sie dann diejenigen aus, bei denen der Fehler minimal ist. Wir hatten jedoch ein ernsteres Problem: Wir mussten 80 Schritte unternehmen (siehe Code), um nahezu optimale Werte zu erhalten. Und dies widerspricht bereits der Idee, beim stochastischen Gradientenabstieg im Vergleich zum Gradientenabstieg Rechenzeit zu sparen. Was kann korrigiert und verbessert werden? Es ist nicht schwer zu erkennen, dass wir in den ersten Iterationen zuversichtlich nach unten gehen und deshalb sollten wir in den ersten Iterationen einen großen Schritt belassen und den Schritt im weiteren Verlauf verringern. Wir werden dies in diesem Artikel nicht tun, da er bereits zu lang ist. Wer möchte, kann sich selbst überlegen, wie das geht, es ist nicht schwer :)

Lassen Sie uns nun mithilfe der Bibliothek einen stochastischen Gradientenabstieg durchführen NumPy (Und stolpern wir nicht über die Steine, die wir zuvor identifiziert haben)

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

Es stellte sich heraus, dass die Werte fast die gleichen waren wie beim Abstieg ohne Verwendung NumPy. Dies ist jedoch logisch.

Finden wir heraus, wie lange die stochastischen Gradientenabstiege gedauert haben.

Code zur Bestimmung der SGD-Berechnungszeit (80 Schritte)

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)

Je weiter in den Wald hinein, desto dunkler die Wolken: Auch hier zeigt die „selbst geschriebene“ Formel das beste Ergebnis. All dies deutet darauf hin, dass es noch subtilere Möglichkeiten zur Nutzung der Bibliothek geben muss NumPy, was Rechenoperationen wirklich beschleunigt. In diesem Artikel werden wir nichts über sie erfahren. In Ihrer Freizeit gibt es bestimmt etwas zum Nachdenken :)

Fassen wir zusammen

Bevor ich zusammenfasse, möchte ich eine Frage beantworten, die höchstwahrscheinlich von unserem lieben Leser gestellt wurde. Warum eigentlich so eine „Folter“ beim Abstieg, warum müssen wir den Berg hinauf und hinunter laufen (hauptsächlich hinunter), um das geschätzte Tiefland zu finden, wenn wir ein so leistungsstarkes und einfaches Gerät in der Hand haben? Form einer analytischen Lösung, die uns sofort an den richtigen Ort teleportiert?

Die Antwort auf diese Frage liegt an der Oberfläche. Jetzt haben wir uns ein sehr einfaches Beispiel angesehen, in dem die wahre Antwort lautet hängt von einem Zeichen ab . Das sieht man im Leben nicht oft, also stellen wir uns vor, wir hätten 2, 30, 50 oder mehr Zeichen. Fügen wir dazu Tausende oder sogar Zehntausende Werte für jedes Attribut hinzu. In diesem Fall kann es sein, dass die analytische Lösung dem Test nicht standhält und versagt. Der Gradientenabstieg und seine Variationen wiederum werden uns langsam aber sicher dem Ziel näher bringen – dem Minimum der Funktion. Und machen Sie sich keine Sorgen um die Geschwindigkeit – wir werden wahrscheinlich nach Möglichkeiten suchen, wie wir die Schrittlänge (also die Geschwindigkeit) einstellen und regulieren können.

Und nun die eigentliche Kurzzusammenfassung.

Erstens hoffe ich, dass das in dem Artikel vorgestellte Material Anfängern von „Datenwissenschaftlern“ dabei helfen wird, zu verstehen, wie einfache (und nicht nur) lineare Regressionsgleichungen gelöst werden können.

Zweitens haben wir uns verschiedene Möglichkeiten zur Lösung der Gleichung angesehen. Jetzt können wir je nach Situation diejenige auswählen, die am besten zur Lösung des Problems geeignet ist.

Drittens haben wir die Leistungsfähigkeit zusätzlicher Einstellungen gesehen, nämlich die Schrittlänge des Gradientenabstiegs. Dieser Parameter kann nicht vernachlässigt werden. Wie oben erwähnt, sollte die Schrittlänge während des Abstiegs geändert werden, um den Berechnungsaufwand zu reduzieren.

Viertens zeigten in unserem Fall „selbst geschriebene“ Funktionen die besten Zeitergebnisse für Berechnungen. Dies ist wahrscheinlich auf die nicht besonders professionelle Nutzung der Bibliotheksfunktionen zurückzuführen NumPy. Aber wie dem auch sei, die folgende Schlussfolgerung liegt nahe. Einerseits lohnt es sich manchmal, etablierte Meinungen zu hinterfragen, andererseits lohnt es sich nicht immer, alles zu verkomplizieren – im Gegenteil, manchmal ist eine einfachere Lösung eines Problems effektiver. Und da unser Ziel darin bestand, drei Ansätze zur Lösung einer einfachen linearen Regressionsgleichung zu analysieren, reichte uns die Verwendung „selbst geschriebener“ Funktionen völlig aus.