Risolvendo l'equazione della regressione lineare semplice

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

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

• Soluzione analitica
• Discesa del gradiente
• Discesa del gradiente stocastica

Per ciascuno dei metodi di risoluzione dell'equazione della retta, l'articolo presenta diverse funzioni, suddivise in quelle scritte senza l'uso della libreria NumPy e quelle che applicano NumPy. Si ritiene che un uso abile NumPy consentirà di ridurre i costi di calcolo.

Tutto il codice riportato nell'articolo è scritto in linguaggio python 2.7 con l'uso di Jupyter Notebook. Il codice sorgente e il file con i dati dell'insieme di dati sono disponibili su GitHub

L'articolo è principalmente rivolto sia ai principianti che a coloro che hanno già iniziato ad approcciarsi a questo vasto tema 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 da X (Tabella n. 1):

Tabella n. 1 «Condizioni dell'esempio»

Consideriamo che i valori — sono i mesi dell'anno, e — sono i ricavi di quel mese. In altre parole, i ricavi dipendono dal mese dell'anno, e — è l'unico fattore da cui dipendono i ricavi.

L'esempio è piuttosto semplice sia dal punto di vista della dipendenza condizionata dei ricavi dal mese dell'anno, sia dal punto di vista del numero di valori, che sono molto pochi. Tuttavia, questa semplificazione permetterà di spiegare, anche se non sempre con facilità, il materiale che i principianti faticano ad assimilare. Inoltre, la semplicità dei numeri consentirà a chi lo desidera di risolvere l'esempio 'su carta' senza grandi sforzi.

Supponiamo che la dipendenza presentata nell'esempio possa essere abbastanza bene approssimata dall'equazione matematica della retta di regressione semplice (lineare) nella forma:

dove — è il mese in cui sono stati ottenuti i ricavi, — è il ricavo corrispondente al mese, e — sono i coefficienti della retta di regressione stimata.

Notiamo che il coefficiente è spesso chiamato coefficiente angolare o gradiente della retta stimata; rappresenta la quantità con cui cambia nell'alterazione .

È chiaro che il nostro obiettivo nell'esempio è quello di trovare i coefficienti nell'equazione e , per i quali le deviazioni dei nostri valori calcolati delle entrate mensili dalle risposte reali, ossia i valori presentati nel campione, saranno minime.

Metodo dei minimi quadrati

Secondo il metodo dei minimi quadrati, la deviazione deve essere calcolata elevandola al quadrato. Questo approccio consente di evitare l'annullamento reciproco delle deviazioni quando hanno 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ò usare la proprietà del valore assoluto, accumulando tutte le deviazioni come positive. Non ci soffermeremo su questo aspetto in dettaglio, ma semplicemente indicheremo che, per comodità nei calcoli, è consuetudine elevare la deviazione al quadrato.

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

dove — è una funzione di approssimazione delle risposte vere (cioè il nostro fatturato calcolato),

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

— è l'indice del campione (il numero del mese in cui si verifica la determinazione della deviazione)

Deriviamo la funzione, definiamo le equazioni delle derivate parziali e prepariamoci a passare alla soluzione analitica. Ma per iniziare faremo una breve introduzione su cosa sia la derivazione e ricorderemo il significato geometrico della derivata.

Derivazione

La derivazione è l'operazione che consente di trovare la derivata di una funzione.

A cosa serve la derivata? La derivata di una funzione caratterizza la velocità di cambiamento della funzione e ci indica la sua direzione. Se la derivata in un dato punto è positiva, la funzione è in crescita, altrimenti la funzione è in diminuzione. E quanto più alto è il valore assoluto della derivata, tanto maggiore è la velocità di cambiamento dei valori della funzione, così come più ripido è l'angolo di inclinazione del grafico della funzione.

Ad esempio, nelle condizioni di un sistema di coordinate cartesiane, il valore della derivata nel punto M(0,0) pari a +25 significa che in un dato punto, spostando il valore a destra di un'unità convenzionale, valore aumenta di 25 unità convenzionali. Nel grafico questo appare come un'inclinazione piuttosto ripida dei valori a partire da un punto dato.

Un altro esempio. Il valore della derivata pari a -0,1 significa che spostandosi di un'unità convenzionale, il valore diminuisce solo di 0,1 unità convenzionali. In questo caso, nel grafico della funzione, possiamo osservare una leggera pendenza verso il basso. Se facciamo un'analogia con una montagna, è come se stessimo lentamente scendendo da un pendio leggero, a differenza del precedente esempio, dove dovevamo affrontare vette molto ripide:)

Pertanto, effettuando la derivazione della funzione rispetto ai coefficienti e , definiamo le equazioni delle derivate parziali di 1° ordine. Dopo aver definito le equazioni, otterremo un sistema di due equazioni, risolvendo il quale saremo in grado di determinare i valori dei coefficienti. e , dove i valori delle derivate parziali corrispondenti nei punti specificati cambiano di una quantità estremamente piccola, mentre nel caso della 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 forma:

l'equazione della derivata parziale di primo ordine rispetto a assumerà la 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, caricheremo, controlleremo la correttezza del caricamento e formatteremo i dati.

Caricamento e formattazione dei dati

Va notato che, dato che per la soluzione analitica, e in seguito per il gradiente e il gradiente stocastico, utilizzeremo il codice in due varianti: usando la libreria NumPy e senza il suo utilizzo, avremo bisogno di una formattazione dei dati adeguata (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 caricato i dati, verificato la correttezza del caricamento e infine formattato i dati, procediamo alla prima visualizzazione. Spesso per questo si utilizza il metodo pairplot della libreria Seaborn. Nel nostro esempio, a causa della limitatezza dei numeri, non ha senso utilizzare la libreria Seaborn. Utilizzeremo la normale libreria Matplotlib e visualizzeremo solo il grafico a dispersione.

Codice del grafico a dispersione

print 'Grafico n. 1 "Dipendenza del fatturato dal mese dell'anno"'

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

Grafico n. 1 «Dipendenza del fatturato dal mese dell'anno»

Soluzione analitica

Utilizzeremo i più semplici strumenti in python e risolveremo il 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 comune, così come i determinanti per e per , dopo di che, dividendo il determinante per per il determinante comune, 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. Tracciamo una linea nel grafico a dispersione in base ai coefficienti trovati.

Codice della retta 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 n. 2 "Risposte corrette e calcolate"

Possiamo osservare il grafico degli scostamenti per ogni mese. Nel nostro caso, non trarremo alcun valore pratico significativo da esso, ma soddisferemo la curiosità su quanto bene l'equazione della regressione lineare semplice caratterizzi la dipendenza dei ricavi dal mese dell'anno.

Codice per il 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 n. 3 "Scostamenti, %"

Non è perfetto, ma abbiamo completato il nostro compito.

Scriviamo una funzione che per determinare i coefficienti e utilizza la libreria NumPy, più precisamente – scriviamo due funzioni: una utilizzando la pseudo-inversa della matrice (sconsigliata 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 esecuzione del 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 esecuzione del calcolo dei coefficienti utilizzando la matrice pseudo-inversa:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Tempo di esecuzione del calcolo dei coefficienti utilizzando l'equazione matriciale:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Su un numero limitato di dati, la funzione «fatta in casa» emerge, trovando i coefficienti con il metodo di Cramer.

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

Discesa del gradiente

Iniziamo definendo che cos'è un gradiente. In termini semplici, un gradiente è un segmento che indica la direzione di massima crescita di una funzione. Analogamente a una salita in montagna, dove il gradiente punta, lì si trova la pendenza più ripida verso la vetta. Sviluppando l'analogia con la montagna, ricordiamo che in realtà abbiamo bisogno della discesa più ripida per raggiungere il basso il più rapidamente possibile, ovvero il minimo, il punto in cui la funzione non cresce e non decresce. In questo punto, la derivata sarà uguale a zero. Pertanto, abbiamo bisogno non di un gradiente, ma di un anti-gradiente. Per trovare l'anti-gradiente, basta moltiplicare il gradiente per -1 (meno uno).

Notiamo che una funzione può avere più minimi, e scendendo in uno di essi utilizzando l'algoritmo proposto in seguito, non saremo in grado di trovare un altro minimo che potrebbe trovarsi più in basso rispetto a quello trovato. Rilassiamoci, questo non ci minaccia! Nel nostro caso, abbiamo a che fare con un unico minimo, poiché la nostra funzione sul grafico rappresenta una normale parabola. E come tutti noi dovremmo sapere bene dal corso di matematica delle scuole superiori, per una parabola esiste un solo minimo.

Dopo aver chiarito a cosa ci servisse il gradiente e che il gradiente è un segmento, ovvero un vettore con coordinate specifiche, che sono i coefficienti stessi. e Possiamo implementare il gradiente discendente.

Prima di iniziare, ti consiglio di leggere qualche frase sull'algoritmo di discesa:

• Definiamo le coordinate dei coefficienti in modo pseudocasuale. e Nel nostro esempio, definiremo i coefficienti vicino allo zero. Questa è una prassi comune, anche se ogni caso potrebbe avere la sua metodologia.
• Dalla coordinata sottraiamo il valore della derivata parziale di primo ordine nel punto . Così, se la derivata è positiva, la funzione cresce. Pertanto, sottraendo il valore della derivata, ci muoveremo in direzione opposta alla crescita, cioè verso la discesa. Se la derivata è negativa, vuol dire che la funzione sta decrescendo in quel punto e sottraendo il valore della derivata, ci muoviamo verso la discesa.
• Eseguiamo un'operazione analoga con la coordinata : sottraiamo il valore della derivata parziale nel punto. .
• Per non saltare il minimo e non volare nel lontano spazio, è necessario impostare la dimensione del passo in direzione della discesa. In generale, si potrebbe scrivere un intero articolo su come impostare correttamente il passo e come modificarlo durante la discesa per ridurre i costi di calcolo. Ma ora ci troviamo di fronte a un compito leggermente diverso, e con il metodo scientifico del 'tentativo ed errore', o come si dice comunemente, in modo empirico, stabiliremo la dimensione del passo.
• Dopo aver sottratto dai coordinate date e i valori delle derivate, otteniamo nuove coordinate e . Effettuiamo il prossimo passo (sottrazione), già dalle coordinate calcolate. E così il ciclo ricomincia, finché non si raggiunge la convergenza richiesta.

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

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

Siamo scesi fino al fondo della fossa delle Marianne e lì abbiamo trovato gli stessi valori dei coefficienti e , come ci si aspettava.

Effettueremo un'altra immersione, solo che questa volta, l'attrezzatura del nostro veicolo sottomarino utilizzerà tecnologie diverse, ovvero la biblioteca 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

Valori dei coefficienti e invariati.

Esaminiamo come è cambiato l'errore durante la discesa del gradiente, ovvero come è variata la somma dei quadrati degli scarti a ogni passo.

Codice per il grafico delle somme dei quadrati degli scarti

print 'Grafico№4 "Somma dei quadrati degli scarti 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 degli scarti', size=16)
plt.show()

Grafico №4 «Somma dei quadrati degli scarti durante la 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.

Infine, 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 l'uso della 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 con l'uso della libreria NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Forse stiamo sbagliando qualcosa, ma di nuovo una semplice funzione 'scritta a mano' che non utilizza una libreria NumPy battete in termini di tempo di esecuzione la funzione che utilizza la libreria NumPy.

Ma non ci fermiamo e ci muoviamo verso l'esplorazione di un altro affascinante modo di risolvere l'equazione della regressione lineare semplice. Ecco a voi!

Discesa del gradiente stocastica

Per comprendere più rapidamente il principio del funzionamento della discesa del gradiente stocastica, è meglio definire le differenze rispetto alla discesa del gradiente ordinaria. Nel caso della discesa del gradiente, nel caso delle derivate dalle e abbiamo utilizzato le somme dei valori di tutte le caratteristiche e delle risposte effettive presenti nel campione (cioè le somme di tutti e ). Nella discesa del gradiente stocastica non utilizzeremo tutti i valori presenti nel campione, ma invece selezioneremo casualmente un cosiddetto indice del campione e utilizzeremo i suoi valori.

Per esempio, se l'indice è stato determinato con il numero 3 (tre), allora prendiamo i valori e , quindi inseriamo i valori nelle equazioni delle derivate e definiamo nuove coordinate. Poi, una volta definite le coordinate, determiniamo di nuovo casualmente l'indice del campionamento, inserendo i valori corrispondenti all'indice nelle equazioni delle derivate parziali per rideterminare le coordinate. e e così via fino al rinverdimento della convergenza. A prima vista, potrebbe sembrare incomprensibile come ciò possa funzionare, eppure funziona. È vero che non l'errore si riduce ad ogni passaggio, ma c'è sicuramente una tendenza.

Quali sono i vantaggi del gradiente stocastico rispetto a quello ordinario? Se abbiamo una dimensione del campione molto grande, misurata in decine di migliaia di valori, è molto più semplice elaborare, ad esempio, un migliaio casuale di essi piuttosto che l'intero campione. È in questo caso che viene avviato il gradiente stocastico. Nel nostro caso, non noteremo ovviamente una grande differenza.

Diamo un'occhiata al codice.

Codice per il gradiente 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 ritroviamo a chiederci: «Come mai?». Abbiamo ottenuto valori diversi per i coefficienti. e Forse il gradiente stocastico ha trovato parametri più ottimali per l'equazione? Purtroppo, no. Basta guardare la somma dei quadrati degli scarti per vedere che con i nuovi valori dei coefficienti, l'errore aumenta. Non ci scoraggiamo. Disegniamo il grafico della variazione dell'errore.

Codice per il grafico della somma dei quadrati degli scarti nel gradiente stocastico

print 'Grafico n. 5 "Somma dei quadrati degli scarti passo dopo passo"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Passaggi (Iterazione)', size=16)
plt.ylabel('Somma dei quadrati degli scarti', size=16)
plt.show()

Grafico n. 5 «Somma dei quadrati degli scarti nel gradiente stocastico»

Guardando il grafico, tutto diventa chiaro e ora correggeremo tutto.

Итак, что же произошло? Произошло следующее. Когда мы выбираем случайным образом месяц, то именно для выбранного месяца наш алгоритм стремится уменьшить ошибку в расчете выручки. Затем выбираем другой месяц и повторяем расчет, но ошибку уменьшаем уже для второго выбранного месяца. А теперь вспомним, что у нас первые два месяца существенно отклоняются от линии уравнения простой линейной регрессии. Это значит, что когда выбирается любой из этих двух месяцев, то уменьшая ошибку каждого из них, наш алгоритм серьезно увеличивает ошибку по всей выборке. Так что же делать? Ответ простой: надо уменьшить шаг спуска. Ведь уменьшив шаг спуска, ошибка так же перестанет «скакать» то вверх, то вниз. Вернее, ошибка «скакать» не перестанет, но будет это делать не так прытко:) Проверим.

Код для запуска SGD с меньшим шагом

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

График №6 «Сумма квадратов отклонений при стохастическом градиентном спуске (80 тыс. шагов)»

I valori dei coefficienti sono migliorati, ma non sono ancora ideali. In teoria, possiamo risolvere questo problema in questo modo. Selezioniamo, ad esempio, per le ultime 1000 iterazioni, i valori dei coefficienti per cui è stato commesso l'errore minimo. Tuttavia, per fare ciò, dovremo registrare anche i valori stessi dei coefficienti. Non lo faremo, ma piuttosto diamo un'occhiata al grafico. Sembra liscio e l'errore sembra diminuire uniformemente. In realtà non è così. Guardiamo 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 per passaggi. 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 per passaggi. 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 (ultimi 1000 passi)»

All'inizio del calo osserviamo una diminuzione dell'errore piuttosto uniforme e ripida. Negli ultimi passi vediamo che l'errore si attesta intorno al valore di 1,475 e in alcuni momenti raggiunge persino questo valore ottimale, ma poi risale... Ripeto, si possono registrare i valori dei coefficienti e , e poi scegliere quelli per i quali l'errore è minimo. Tuttavia, abbiamo riscontrato un problema più serio: abbiamo dovuto fare 80.000 passi (vedi codice), per ottenere valori vicini agli ottimali. Questo contraddice già l'idea di risparmiare tempo di calcolo rispetto al gradiente. Cosa si può correggere e migliorare? Non è difficile notare che nelle prime iterazioni scendiamo con decisione e, di conseguenza, dovremmo mantenere un grande passo nelle prime iterazioni e diminuire il passo man mano che avanziamo. Non lo faremo in questo articolo - è già abbastanza lungo. Chi lo desidera può pensarci da solo, non è difficile 🙂

Ora eseguiremo la discesa gradiente stocastica utilizzando la libreria NumPy (e non ci inciamperemo sui sassi che abbiamo identificato in precedenza)

Codice per la discesa gradiente stocastica (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 identici a quelli della discesa senza utilizzo NumPy. Tuttavia, questo ha senso.

Scopriamo quanto tempo ci hanno impiegato le discese gradiente stocastiche.

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

print ' 33[1m' + ' 33[4m' +
"Tempo di esecuzione della discesa gradiente stocastica 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 gradiente stocastica 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 si addentra nel bosco, più scure diventano le nuvole: eppure, la formula 'fatta in casa' mostra ancora una volta i migliori risultati. Tutto ciò invita a riflettere su come debbano esistere modi ancora più raffinati per utilizzare la libreria. NumPy, che accelerano davvero le operazioni di calcolo. In questo articolo non ne parleremo più. Sarà qualcosa 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 fin dei conti, queste "fatiche" con le discese, perché dobbiamo salire e scendere da una montagna (soprattutto scendere) per trovare la tanto cercata valle, se nelle nostre mani abbiamo uno strumento così potente e semplice, sotto forma di una soluzione analitica, che ci teletrasporta istantaneamente dove ci serve?

La risposta a questa domanda è evidente. Finora abbiamo esaminato un esempio molto semplice, in cui la risposta vera dipende da un singolo attributo . Nella vita si incontra raramente qualcosa del genere, quindi immaginiamo di avere 2, 30, 50 o più caratteristiche. Aggiungiamo a questo migliaia, o addirittura decine di migliaia di valori per ogni caratteristica. In questo caso, una soluzione analitica potrebbe non resistere alla prova e dare errore. A sua volta, il gradiente discendente 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 esploreremo ancora metodi 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 presentato nell'articolo aiuti i principianti 'data scientist' a comprendere come risolvere le equazioni della regressione lineare semplice (e non solo).

In secondo luogo, abbiamo esaminato diversi metodi di risoluzione delle equazioni. Ora, a seconda della situazione, possiamo scegliere quello che meglio si adatta alla risoluzione del compito assegnato.

In terzo luogo, abbiamo osservato l'importanza delle impostazioni aggiuntive, in particolare la lunghezza del passo del gradiente. Questo parametro non può essere trascurato. Come sottolineato in precedenza, per ridurre i costi di calcolo, la lunghezza del passo dovrebbe essere modificata durante il processo di discesa.

In quarto luogo, nel nostro caso, le funzioni 'fatte a mano' hanno mostrato i migliori tempi di calcolo. Probabilmente ciò è dovuto all'uso non molto professionale delle capacità della libreria. NumPy. Tuttavia, la conclusione è la seguente. Da un lato, è a volte opportuno mettere in discussione le opinioni consolidate, e dall'altro, non sempre è necessario complicare le cose: anzi, a volte un modo più semplice si dimostra più efficace per risolvere il problema. E poiché il nostro obiettivo era esaminare tre approcci alla soluzione dell'equazione di regressione lineare semplice, l'uso delle funzioni ‘fatte a mano’ è 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