Risolviamo l'equazione della regressione lineare semplice

L'articolo esplora diversi metodi per determinare l'equazione matematica della retta di regressione semplice (bivaria).

Tutti i metodi discussi qui per risolvere l'equazione si basano sul metodo dei minimi quadrati. Denominiamo i metodi come segue:

• Soluzione analitica
• Discesa del gradiente
• Discesa del gradiente stocastico

Per ciascuno dei metodi di risoluzione dell'equazione della retta, l'articolo fornisce diverse funzioni, che si dividono principalmente in quelle scritte senza utilizzare la libreria NumPy e quelle che utilizzano per i calcoli NumPy. Si considera che un uso abile NumPy consentirà di ridurre i costi di calcolo.

Tutto il codice fornito nell'articolo è scritto in lingua python 2.7 utilizzando Jupyter Notebook. Il codice sorgente e il file con i dati del campione sono disponibili su GitHub

L'articolo è principalmente orientato sia ai principianti sia a coloro che hanno già cominciato a esplorare un'ampia area dell'intelligenza artificiale: il machine learning.

Per illustrare il materiale utilizziamo un esempio molto semplice.

Condizioni dell'esempio

Abbiamo cinque valori che caratterizzano la dipendenza Y di X (Tabella n. 1):

Tabella n. 1 «Condizioni dell'esempio»

Consideriamo che i valori rappresentano i mesi dell'anno, e che è il fatturato di quel mese. In altre parole, il fatturato dipende dal mese dell'anno, e è l'unico indicatore dal quale dipende il fatturato.

L'esempio è un po' scarso, sia dal punto di vista della dipendenza condizionale del fatturato dal mese dell'anno, sia per quanto riguarda il numero di valori: ce ne sono davvero pochi. Tuttavia, questa semplificazione permetterà, per così dire, di spiegare, non sempre facilmente, il materiale che i principianti devono assimilare. Inoltre, la semplicità dei numeri consentirà a chi lo desidera di risolvere l'esempio su 'carta'.

Supponiamo che la dipendenza presentata nell'esempio possa essere abbastanza bene approssimata da un'equazione matematica della retta di regressione semplice (bivaria) di tipo:

dove – è il mese in cui è stato ottenuto il fatturato, – è il fatturato corrispondente al mese, e – sono i coefficienti della retta di regressione stimata.

Si osservi che il coefficiente è spesso chiamato coefficiente angolare o gradiente della retta stimata; rappresenta l'entità con cui cambierà in caso di modifica .

È evidente che il nostro compito nell'esempio è scegliere i coefficienti nell'equazione e , tali che le deviazioni dei nostri valori calcolati di fatturato mensile dalle risposte reali, cioè i valori forniti nel campione, siano minime.

Il metodo dei minimi quadrati

Secondo il metodo dei minimi quadrati, la deviazione deve essere calcolata elevandola al quadrato. Questo approccio consente di evitare la compensazione reciproca delle deviazioni, nel caso in cui esse abbiano segni opposti. Ad esempio, se in un caso la deviazione è +5 (più cinque), mentre nell'altro -5 (meno cinque), la somma delle deviazioni si annullerà e risulterà 0 (zero). Non è necessario elevare la deviazione al quadrato, ma si può utilizzare la proprietà del valore assoluto e così tutte le deviazioni saranno positive e si accumuleranno. Non ci soffermeremo su questo punto in dettaglio, ma semplicemente segnaleremo che, per comodità nei calcoli, si tende a elevare la deviazione al quadrato.

Ecco come appare la formula, con cui definiremo la somma minima dei quadrati delle deviazioni (errori):

dove — è la funzione di approssimazione delle risposte reali (cioè il fatturato calcolato da noi),

— sono le risposte reali (il fatturato fornito nel campione),

— è l'indice del campione (il numero del mese in cui viene determinata la deviazione)

Deriviamo la funzione, definiamo le equazioni delle derivate parziali e saremo pronti per passare alla soluzione analitica. Ma prima faremo una breve excursus su cosa sia la differenziazione e ricorderemo il significato geometrico della derivata.

Differenziazione

La differenziazione è l'operazione di determinare la derivata di una funzione.

A cosa serve la derivata? La derivata di una funzione caratterizza la velocità di cambiamento della funzione e indica la sua direzione. Se la derivata in un dato punto è positiva, allora la funzione cresce, nel caso opposto — la funzione decresce. E quanto maggiore è il valore assoluto della derivata, tanto più alta è la velocità di cambiamento dei valori della funzione e tanto più ripido è l'angolo di inclinazione del grafico della funzione.

Ad esempio, nelle condizioni di un sistema di coordinate cartesiano, il valore della derivata nel punto M(0,0) uguale a +25 significa che nel punto dato, spostando il valore a destra di un'unità condizionata, il valore aumenta di 25 unità arbitrarie. Nel grafico, questo appare come un angolo di salita piuttosto ripido dei valori. da un punto specificato.

Un altro esempio. Il valore della derivata pari a -0,1 significa che al variare di un'unità arbitraria, il valore diminuisce solo di 0,1 unità arbitrarie. In questo caso, nel grafico della funzione, possiamo osservare una leggera inclinazione verso il basso. Tracciando un'analogia con una montagna, sembra che stiamo scendendo molto lentamente lungo un pendio dolce, a differenza dell'esempio precedente in cui dovevamo affrontare picchi molto ripidi :)

Pertanto, eseguendo la differenziazione della funzione rispetto ai coefficienti e , definiremo le equazioni delle derivate parziali di primo ordine. Dopo aver definito le equazioni, otteniamo un sistema di due equazioni, risolvendo le quali potremo selezionare tali valori dei coefficienti e , per i quali i valori delle derivate corrispondenti nei punti specificati cambiano di una quantità molto, molto piccola, mentre nel caso di una soluzione analitica non cambiano affatto. In altre parole, la funzione di errore con i coefficienti trovati raggiungerà un minimo, poiché i valori delle derivate parziali in questi punti saranno pari a zero.

Quindi, secondo le regole di differenziazione, l'equazione della derivata parziale di primo ordine rispetto al coefficiente assumerà la seguente forma:

l'equazione della derivata parziale di primo ordine rispetto a assumerà la seguente forma:

Alla fine abbiamo ottenuto un sistema di equazioni che ha una soluzione analitica piuttosto semplice:

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*}

Prima di risolvere l'equazione, carichiamo preliminarmente, verifichiamo la correttezza del caricamento e formattiamo i dati.

Caricamento e formattazione dei dati

È importante notare che, poiché per la soluzione analitica, e successivamente per il gradiente e il gradiente stocastico, utilizzeremo il codice in due varianti: con l'uso della libreria NumPy e senza il suo utilizzo, sarà necessario il formato appropriato per i dati (vedi codice).

Codice per il caricamento e l'elaborazione dei dati

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

Visualizzazione

Ora, dopo aver, da un lato, caricato i dati, dall'altro lato, verificato la correttezza del caricamento e infine formattato i dati, procederemo alla prima visualizzazione. Spesso per questo si utilizza il metodo pairplot una libreria Seaborn. Nel nostro esempio, a causa della limitatezza dei numeri, non ha senso utilizzare una libreria Seaborn. Utilizzeremo una libreria normale Matplotlib e guarderemo solo al diagramma a dispersione.

Codice del diagramma a dispersione

print 'Grafico №1 "Dipendenza delle entrate dal mese dell'anno"'

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

Grafico №1 «Dipendenza delle entrate dal mese dell'anno»

Soluzione analitica

Utilizzeremo gli strumenti più comuni in python e risolveremo un sistema di equazioni:

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*}

Secondo la regola di Cramer troveremo il determinante complessivo e anche i determinanti per e per , dopodiché, dividendo il determinante per per il determinante complessivo, troveremo il coefficiente , analogamente troveremo il coefficiente .

Codice della soluzione analitica

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

Ecco cosa abbiamo ottenuto:

Quindi, i valori dei coefficienti sono stati trovati e la somma dei quadrati degli scostamenti è stata stabilita. Disegneremo sul diagramma a dispersione una retta in base ai coefficienti trovati.

Codice della linea di regressione

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

Grafico №2 «Risposte corrette e calcolate»

Possiamo dare un'occhiata al grafico degli scostamenti per ogni mese. Nel nostro caso, non trarremo alcun valore pratico significativo, ma soddisferemo la curiosità riguardo a quanto bene l'equazione della semplice regressione lineare caratterizzi la dipendenza delle entrate dal mese dell'anno.

Codice del grafico degli scostamenti

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

Grafico №3 «Scostamenti, %»

Non è perfetto, ma abbiamo svolto il nostro compito.

Scriveremo una funzione che per determinare i coefficienti e utilizza una libreria NumPy, più precisamente, scriveremo due funzioni: una utilizzando la matrice pseudoinversa (non raccomandata in pratica, poiché il processo è computazionalmente complesso e instabile), l'altra utilizzando l'equazione matriciale.

Codice della soluzione analitica (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

Confrontiamo il tempo impiegato per determinare i coefficienti e , secondo i 3 metodi presentati.

Codice per il calcolo dei tempi di elaborazione

print ' 33[1m' + ' 33[4m' + "Tempo di calcolo dei coefficienti senza utilizzare la libreria NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Tempo di calcolo dei coefficienti utilizzando la pseudo-inversa della matrice:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Tempo di calcolo dei coefficienti utilizzando l'equazione matriciale:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Con un numero limitato di dati, la funzione «scritta a mano» risulta vincente, trovando i coefficienti mediante il metodo di Kramer.

Ora possiamo passare ad altri metodi per calcolare i coefficienti. e .

Discesa del gradiente

Per cominciare, definiamo che cos'è il gradiente. In parole semplici, il gradiente è un vettore che indica la direzione di massimo incremento di una funzione. Per analogia con una salita in montagna, dove il gradiente è rivolto nel punto di massima pendenza verso la vetta. Continuando con l'esempio della montagna, ci rendiamo conto che in realtà abbiamo bisogno della discesa più ripida per raggiungere più velocemente la valle, ovvero il minimo — il luogo in cui la funzione non cresce né decresce. In questo punto, la derivata si annulla. Di conseguenza, ci serve non il gradiente, ma l'antigradiente. Per trovare l'antigradiente è sufficiente moltiplicare il gradiente per -1 (meno uno).

Notiamo che la funzione può avere più minimi, e una volta scesi in uno di essi seguendo l'algoritmo proposto successivamente, non saremo in grado di trovare un altro minimo che potrebbe trovarsi sotto a quello già raggiunto. Rilassiamoci, non è un problema che ci riguarda! Nel nostro caso, affrontiamo un solo minimo, poiché la nostra funzione secondo il grafico è una classica parabola. E come dovremmo sapere da scuola — la parabola ha un solo minimo.

Dopo aver chiarito perché ci serve il gradiente e che il gradiente è un vettore con coordinate specifiche, che sono proprio i coefficienti e possiamo implementare la discesa del gradiente.

Prima di iniziare, consiglio di leggere brevemente alcune righe riguardo l'algoritmo di discesa:

• Definiamo casualmente le coordinate dei coefficienti e . Nel nostro esempio, definiremo i coefficienti nelle vicinanze dello zero. Questa è una pratica comune, ma per ogni caso potrebbe essere prevista una propria prassi.
• Dalla coordinata sottraiamo il valore della derivata parziale di primo ordine nel punto . Così, se la derivata è positiva, la funzione cresce. Di conseguenza, sottraendo il valore della derivata, ci muoveremo nella direzione opposta alla crescita, cioè verso il declino. Se la derivata è negativa, significa che la funzione in quel punto decresce e sottraendo il valore della derivata ci muoviamo verso il declino.
• Effettuiamo un'operazione analoga con la coordinata : sottraiamo il valore della derivata parziale nel punto .
• Per evitare di saltare il minimo e di volare verso lo spazio lontano, è necessario stabilire la dimensione del passo nella direzione del declino. In generale, si potrebbe scrivere un intero articolo su come impostare correttamente il passo e come cambiarlo durante il processo di discesa, per ridurre i costi di calcolo. Ma ora abbiamo un compito un po’ diverso e, attraverso un metodo scientifico di "tentativo" o, come si dice nel linguaggio comune, per via empirica, stabiliremo la dimensione del passo.
• Dopo aver sottratto i valori delle derivate dalle coordinate date e otteniamo nuove coordinate e . Facciamo il passo successivo (sottrazione), già dalle coordinate calcolate. E così il ciclo si avvia ancora e ancora, fino a quando non viene raggiunta la convergenza richiesta.

Tutto! Ora siamo pronti a partire alla ricerca della fossa più profonda della Fossa delle Marianne. Iniziamo.

Codice per la discesa del gradiente

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

Ci siamo immersi nel fondo della Fossa delle Marianne e lì abbiamo trovato gli stessi valori dei coefficienti e , come ci si aspettava.

Faremo un'altra immersione, solo che questa volta, l'allestimento del nostro apparecchio sottomarino utilizzerà altre tecnologie, ovvero la libreria NumPy.

Codice per la discesa del gradiente (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

I valori dei coefficienti e rimangono invariati.

Osserviamo come è cambiato l'errore durante la discesa del gradiente, cioè come è cambiata la somma dei quadrati delle deviazioni ad ogni passo.

Codice per il grafico delle somme dei quadrati delle deviazioni

print 'Grafico№4 "Somma dei quadrati delle deviazioni passo dopo passo"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Passi (Iterazione)', size=16)
plt.ylabel('Somma dei quadrati delle deviazioni', size=16)
plt.show()

Grafico n.° 4 «Somma dei quadrati delle deviazioni durante il discesa del gradiente»

Nel grafico vediamo che ad ogni passo l'errore diminuisce e, dopo un certo numero di iterazioni, osserviamo una linea praticamente orizzontale.

Per concludere, valutiamo la differenza nel tempo di esecuzione del codice:

Codice per determinare il tempo di calcolo della discesa del gradiente

print ' 33[1m' + ' 33[4m' + "Tempo di esecuzione della discesa del gradiente senza utilizzare la libreria 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' + "Tempo di esecuzione della discesa del gradiente utilizzando la libreria NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Forse stiamo facendo qualcosa di sbagliato, ma di nuovo una semplice funzione "fatta in casa" che non utilizza la libreria NumPy precede nel tempo di esecuzione la funzione che utilizza la libreria NumPy.

Ma non ci fermiamo, stiamo avanzando verso l'apprendimento di un altro affascinante modo di risolvere l'equazione della regressione lineare semplice. Presentiamo!

Discesa del gradiente stocastico

Per capire più rapidamente il principio di funzionamento della discesa del gradiente stocastico, è meglio definirne le differenze rispetto alla discesa del gradiente standard. Nel caso della discesa del gradiente, nelle equazioni delle derivate da e utilizzavamo le somme dei valori di tutte le caratteristiche e delle risposte reali presenti nel campione (cioè le somme di tutti e ). Nella discesa del gradiente stocastico non utilizzeremo tutti i valori presenti nel campione, ma invece selezioneremo casualmente quello che viene chiamato indice del campione e utilizzeremo i suoi valori.

Ad esempio, se l'indice è stato determinato come numero 3 (tre), allora prendiamo i valori e , poi sostituiamo i valori nelle equazioni delle derivate e determiniamo nuove coordinate. Successivamente, una volta definite le coordinate, determiniamo di nuovo casualmente l'indice del campione, sostituiamo i valori corrispondenti all'indice nelle equazioni delle derivate parziali e rideterminiamo le coordinate e e così via fino alla convergenza. A prima vista, potrebbe sembrare impossibile che questo funzioni, ma invece funziona. Vale la pena notare che l'errore non diminuisce ad ogni passo, ma c'è sicuramente una tendenza.

Quali sono i vantaggi del gradient descent stocastico rispetto a quello normale? Nel caso in cui abbiamo un grande campione che si misura in decine di migliaia di valori, è molto più semplice elaborare, ad esempio, una mille casuale di essi piuttosto che l'intero campione. In questo caso, inizia il gradient descent stocastico. Nel nostro caso, ovviamente, non noteremo una grande differenza.

Guardiamo il codice.

Codice per il gradient descent stocastico

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

Osserviamo attentamente i coefficienti e ci poniamo la domanda: «Ma come mai?». Abbiamo ottenuto altri valori per i coefficienti e . Forse il gradient descent stocastico ha trovato parametri più ottimali per l'equazione? Purtroppo no. Basta guardare la somma dei quadrati delle deviazioni e vedere che con i nuovi valori dei coefficienti, l'errore è maggiore. Non affrettiamoci a scoraggiarci. Costruiamo un grafico del cambiamento dell'errore.

Codice per il grafico della somma dei quadrati delle deviazioni nel gradient descent stocastico

print 'Grafico n. 5 "Somma dei quadrati delle deviazioni per step"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Passi (Iterazione)', size=16)
plt.ylabel('Somma dei quadrati delle deviazioni', size=16)
plt.show()

Grafico n. 5 «Somma dei quadrati delle deviazioni nel gradient descent stocastico»

Guardando il grafico, tutto si chiarisce e ora correggeremo tutto.

Quindi, cosa è successo? È successo quanto segue. Quando scegliamo casualmente un mese, per quel mese specifico, il nostro algoritmo cerca di ridurre l'errore nel calcolo del fatturato. Poi scegliamo un altro mese e ripetiamo il calcolo, ma ora riduciamo l'errore per il secondo mese scelto. E ora ricordiamo che i primi due mesi si discostano notevolmente dalla linea dell'equazione della regressione lineare semplice. Ciò significa che quando si sceglie uno di questi due mesi, riducendo l'errore di ciascuno di essi, il nostro algoritmo aumenta notevolmente l'errore sull'intero campione. Cosa fare? La risposta è semplice: dobbiamo ridurre il passo di discesa. Infatti, diminuendo il passo di discesa, l'errore smetterà di

Codice per avviare SGD con un passo minore

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

Grafico n. 6 «Somma dei quadrati delle deviazioni nel gradient descent stocastico (80 mila passi)»

I valori dei coefficienti sono migliorati, ma non sono comunque ideali. Ipotesi, si potrebbe sistemare in questo modo. Scegliamo, ad esempio, nei recenti 1000 iterazioni, i valori dei coefficienti con cui è stato commesso il minimo errore. Tuttavia, per fare ciò, dovremo registrare anche i valori stessi dei coefficienti. Non lo faremo, ma rivolgiamo la nostra attenzione al grafico. Appare liscio, e l'errore sembra diminuire uniformemente. In realtà non è così. Esaminiamo le prime 1000 iterazioni e confrontiamole con le ultime.

Codice per il grafico SGD (prime 1000 fasi)

print 'Grafico n. 7 "Somma dei quadrati delle deviazioni in fasi. Prime 1000 iterazioni"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Fasi (Iterazione)', size=16)
plt.ylabel('Somma dei quadrati delle deviazioni', size=16)
plt.show()

print 'Grafico n. 7 "Somma dei quadrati delle deviazioni in fasi. Ultime 1000 iterazioni"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Fasi (Iterazione)', size=16)
plt.ylabel('Somma dei quadrati delle deviazioni', size=16)
plt.show()

Grafico n. 7 «Somma dei quadrati delle deviazioni SGD (prime 1000 fasi)»

Grafico n. 8 «Somma dei quadrati delle deviazioni SGD (ultime 1000 fasi)»

All'inizio della discesa osserviamo un notevole e uniforme calo dell'errore. Nelle ultime iterazioni vediamo che l'errore oscilla attorno a un valore di 1,475 e in alcuni momenti raggiunge anche questo valore ottimale, ma poi comunque tende a risalire… Ripeto, si potrebbero registrare i valori dei coefficienti e , e poi scegliere quelli per cui l'errore è minimo. Tuttavia, abbiamo avuto un problema più serio: abbiamo dovuto effettuare 80.000 passi (vedi codice) per ottenere valori vicini all'ottimale. Questo contraddice già l'idea di risparmiare tempo di calcolo con il gradiente stocastico rispetto a quello classico. Cosa si potrebbe correggere e migliorare? Non è difficile notare che nelle prime iterazioni andiamo decisamente verso il basso e, pertanto, dovremmo mantenere un grande passo all'inizio e ridurlo man mano che procediamo. Non lo faremo in questo articolo — è già abbastanza lungo. Chi è interessato può riflettere su come farlo, non è difficile 🙂

Ora eseguiamo la discesa del gradiente stocastico, utilizzando la libreria NumPy (e non ci inciampiamo sulle pietre che abbiamo identificato in precedenza)

Codice per la discesa del gradiente stocastico (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

I valori ottenuti sono stati quasi gli stessi di quelli ottenuti senza l'uso NumPy. D'altronde, è logico.

Scopriamo quanto tempo ci hanno richiesto le discese del gradiente stocastico.

Codice per determinare il tempo di calcolo del SGD (80.000 passi)

print ' 33[1m' + ' 33[4m' +
"Tempo di esecuzione della discesa del gradiente stocastico senza utilizzare la libreria 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' +
"Tempo di esecuzione della discesa del gradiente stocastico utilizzando la libreria NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

Più ci allontaniamo nel bosco, più scure sono le nuvole: di nuovo la formula 'scritta a mano' mostra il miglior risultato. Tutto ciò porta a riflettere sul fatto che devono esserci modi ancora più raffinati di utilizzare la libreria NumPy, che accelerano realmente le operazioni di calcolo. In questo articolo non avremo modo di scoprirli. Sarà un pensiero su cui riflettere nel tempo libero:)

Riassumiamo

Prima di riassumere, vorrei rispondere a una domanda che probabilmente è sorta nel nostro caro lettore. A cosa servono, in fondo, tali "sofferenze" con le discesa, perché dobbiamo salire e scendere (principalmente scendere) da una montagna per trovare la tanto desiderata valle, se nelle nostre mani abbiamo uno strumento così potente e semplice, sotto forma di soluzione analitica, che ci teletrasporta istantaneamente nel luogo desiderato?

La risposta a questa domanda è evidente. Abbiamo trattato un esempio molto semplice, in cui la risposta reale dipende da un solo attributo . Nella vita, situazioni di questo tipo non si incontrano spesso, quindi immaginiamo di avere 2, 30, 50 o più attributi. Aggiungiamo a questo migliaia, se non decine di migliaia di valori per ciascun attributo. In questo caso, la soluzione analitica potrebbe non reggere il colpo e andare in errore. D'altra parte, la discesa del gradiente e le sue variazioni ci avvicineranno lentamente ma con certezza all'obiettivo — il minimo della funzione. E per quanto riguarda la velocità, non preoccupatevi: sicuramente analizzeremo ancora una volta i modi che ci permetteranno di impostare e regolare la lunghezza del passo (cioè la velocità).

E ora, in effetti, un breve riassunto.

Innanzitutto, spero che il materiale esposto nell'articolo possa aiutare i neo "data scientist" a comprendere come risolvere le equazioni della regressione lineare semplice (e non solo).

In secondo luogo, abbiamo esaminato diversi metodi per risolvere l'equazione. Ora, a seconda della situazione, possiamo scegliere quello che si adatta meglio alla soluzione del problema proposto.

In terzo luogo, abbiamo visto la forza di ulteriori impostazioni, in particolare la lunghezza del passo del gradiente. Questo parametro non può essere trascurato. Come accennato sopra, per ridurre i costi di calcolo, è consigliabile modificare la lunghezza del passo durante la discesa.

In quarto luogo, nel nostro caso, le funzioni "personalizzate" hanno mostrato un miglior risultato temporale nei calcoli. Probabilmente, questo è legato a un uso non sempre professionale delle potenzialità della libreria. NumPy. Ma in ogni caso, la conclusione che si impone è la seguente. Da un lato, a volte è opportuno mettere in discussione le opinioni consolidate, dall'altro, non sempre bisogna complicare tutto – a volte un approccio più semplice si rivela più efficace. E dato che il nostro obiettivo era analizzare tre approcci alla risoluzione dell'equazione della regressione lineare semplice, l'uso delle funzioni "personalizzate" è stato più che sufficiente.

Letteratura (o qualcosa di simile)

1. Regressione lineare

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

2. Metodo dei minimi quadrati

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Derivata

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradiente

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Discesa del gradiente

habr.com/ru/post/471458

habr.com/ru/post/307312

artemarakcheev.com/2017-12-31/linear_regression

6. Biblioteca NumPy

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

Fonte: habr.com