Risoluzione dell'equazione della regressione lineare semplice

L'articolo discute diversi modi per determinare l'equazione matematica di una semplice retta di regressione (accoppiata).

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

• Soluzione analitica
• Discesa gradiente
• Discesa stocastica del gradiente

Per ogni metodo di risoluzione dell'equazione di una retta, l'articolo fornisce varie funzioni, che si dividono principalmente in quelle scritte senza utilizzare la libreria NumPy e quelli che usano per i calcoli NumPy. Si ritiene che l'uso abile NumPy ridurrà i costi informatici.

Tutto il codice fornito nell'articolo è scritto nella lingua python 2.7 con Notebook Jupyter. Il codice sorgente e il file con i dati di esempio vengono pubblicati su Github

L'articolo è più rivolto sia ai principianti che a coloro che hanno già iniziato gradualmente a padroneggiare lo studio di una sezione molto ampia dell'intelligenza artificiale: l'apprendimento automatico.

Per illustrare il materiale, utilizziamo un esempio molto semplice.

Condizioni di esempio

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

Tabella n.1 “Condizioni di esempio”

Assumeremo che i valori è il mese dell'anno e - entrate di questo mese. In altre parole, le entrate dipendono dal mese dell'anno e - l'unico segno da cui dipendono le entrate.

L'esempio è così così, sia dal punto di vista della dipendenza condizionale delle entrate dal mese dell'anno, sia dal punto di vista del numero di valori: ce ne sono pochissimi. Tuttavia, tale semplificazione consentirà, come si suol dire, di spiegare, non sempre con facilità, il materiale assimilato dai principianti. E anche la semplicità dei numeri consentirà a chi lo desidera di risolvere l'esempio su carta senza notevoli costi di manodopera.

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

dove è il mese in cui sono state ricevute le entrate, — entrate corrispondenti al mese, и sono i coefficienti di regressione della retta stimata.

Si noti che il coefficiente spesso chiamato pendenza o gradiente della linea stimata; rappresenta l'importo con cui il quando cambia .

Ovviamente, il nostro compito nell'esempio è selezionare tali coefficienti nell'equazione и , a cui le deviazioni dei nostri valori di entrate calcolati per mese dalle risposte vere, ad es. i valori presentati nel campione saranno minimi.

Metodo dei minimi quadrati

Secondo il metodo dei minimi quadrati, la deviazione dovrebbe essere calcolata elevandola al quadrato. Questa tecnica consente di evitare la reciproca cancellazione delle deviazioni se hanno segni opposti. Ad esempio, se in un caso la deviazione è +5 (più cinque), e nell'altro -5 (meno cinque), la somma delle deviazioni si annullerà a vicenda e ammonterà a 0 (zero). È possibile non elevare al quadrato la deviazione, ma utilizzare la proprietà del modulo e quindi tutte le deviazioni saranno positive e si accumuleranno. Non ci soffermeremo su questo punto in dettaglio, ma indicheremo semplicemente che per comodità di calcolo è consuetudine elevare al quadrato la deviazione.

Ecco come appare la formula con la quale determineremo la somma minima delle deviazioni quadrate (errori):

dove è una funzione di approssimazione delle risposte vere (ovvero delle entrate che abbiamo calcolato),

sono le risposte vere (entrate fornite nel campione),

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

Differenziamo la funzione, definiamo le equazioni alle derivate parziali e prepariamoci a passare alla soluzione analitica. Ma prima facciamo una breve escursione su cos'è la differenziazione e ricordiamo il significato geometrico della derivata.

Differenziazione

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

A cosa serve il derivato? La derivata di una funzione caratterizza la velocità di variazione della funzione e ci dice la sua direzione. Se la derivata in un dato punto è positiva la funzione aumenta, altrimenti diminuisce. E maggiore è il valore della derivata assoluta, maggiore è il tasso di variazione dei valori della funzione, nonché più ripida la pendenza del grafico della funzione.

Ad esempio, nelle condizioni di un sistema di coordinate cartesiane, il valore della derivata nel punto M(0,0) è uguale a + 25 significa che in un dato punto, quando il valore viene spostato a destra da un'unità convenzionale, valore aumenta di 25 unità convenzionali. Sul grafico sembra un aumento abbastanza marcato dei valori da un dato punto.

Un altro esempio. Il valore del derivato è uguale all'0,1 ottobre significa che quando spostato per una unità convenzionale, valore diminuisce solo di 0,1 unità convenzionali. Allo stesso tempo, sul grafico della funzione, possiamo osservare una pendenza verso il basso appena percettibile. Facendo un'analogia con una montagna, è come se stessimo discendendo molto lentamente un dolce pendio da una montagna, a differenza dell'esempio precedente, dove dovevamo scalare cime molto ripide :)

Quindi, dopo aver differenziato la funzione per probabilità и , definiamo equazioni alle derivate parziali del 1° ordine. Dopo aver determinato le equazioni, riceveremo un sistema di due equazioni, risolvendo il quale potremo selezionare tali valori dei coefficienti и , per il quale i valori delle derivate corrispondenti in determinati punti cambiano di una quantità molto, molto piccola, e nel caso di una soluzione analitica non cambiano affatto. In altre parole, la funzione di errore ai coefficienti trovati raggiungerà il minimo, poiché i valori delle derivate parziali in questi punti saranno pari a zero.

Quindi, secondo le regole della derivazione, l'equazione delle derivate parziali del 1° ordine rispetto al coefficiente assumerà la forma:

Equazione alle derivate parziali del 1° ordine rispetto a assumerà la forma:

Di conseguenza, abbiamo ricevuto un sistema di equazioni che ha una soluzione analitica abbastanza semplice:

inizio{equazione*}
iniziare{casi}
na + bsomma_limiti_{i=1}^nx_i — somma_limiti_{i=1}^ny_i = 0

somma_limiti_{i=1}^nx_i(a +b sommalimiti_{i=1}^nx_i — somma_limiti_{i=1}^ny_i) = 0
conclusione{casi}
fine{equazione*}

Prima di risolvere l'equazione precarichiamo, controlliamo che il caricamento sia corretto e formattiamo i dati.

Caricamento e formattazione dei dati

Da notare che poiché per la soluzione analitica, e successivamente per la discesa del gradiente e stocastico, utilizzeremo il codice in due varianti: utilizzando la libreria NumPy e senza utilizzarlo, avremo bisogno di un'adeguata formattazione dei dati (vedi codice).

Codice di caricamento ed elaborazione 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 prima caricato i dati, poi verificato la correttezza del caricamento e infine formattato i dati, effettueremo la prima visualizzazione. Il metodo spesso utilizzato per questo è trama a coppie Biblioteca Seaborn. Nel nostro esempio, a causa del numero limitato, non ha senso utilizzare la biblioteca Seaborn. Utilizzeremo la biblioteca normale matplotlib e basta guardare il grafico a dispersione.

Codice del grafico a dispersione

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

Grafico n. 1 “Dipendenza delle entrate dal mese dell'anno”

Soluzione analitica

Usiamo gli strumenti più comuni in python e risolvi il sistema di equazioni:

inizio{equazione*}
iniziare{casi}
na + bsomma_limiti_{i=1}^nx_i — somma_limiti_{i=1}^ny_i = 0

somma_limiti_{i=1}^nx_i(a +b sommalimiti_{i=1}^nx_i — somma_limiti_{i=1}^ny_i) = 0
conclusione{casi}
fine{equazione*}

Secondo la regola di Cramer troveremo il determinante generale, così come i determinanti di e , dopodiché, dividendo il determinante per al determinante generale: trova il coefficiente , allo stesso modo troviamo 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, sono stati trovati i valori dei coefficienti, è stata stabilita la somma delle deviazioni al quadrato. Disegniamo una linea retta sull'istogramma di dispersione 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 n.2 “Risposte corrette e calcolate”

Puoi guardare il grafico della deviazione per ogni mese. Nel nostro caso, non ne trarremo alcun valore pratico significativo, ma soddisferemo la nostra curiosità su quanto bene la semplice equazione di regressione lineare caratterizzi la dipendenza delle entrate dal mese dell'anno.

Codice del grafico delle deviazioni

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

Non perfetto, ma abbiamo completato il nostro compito.

Scriviamo una funzione che, per determinare i coefficienti и utilizza la libreria NumPy, più precisamente, scriveremo due funzioni: una utilizzando una matrice pseudoinversa (sconsigliata nella pratica, poiché il processo è computazionalmente complesso e instabile), l'altra utilizzando un'equazione di matrice.

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 и , secondo i 3 metodi presentati.

Codice per il calcolo del tempo di calcolo

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)

Con una piccola quantità di dati, viene fuori una funzione “autoscritta”, che trova i coefficienti utilizzando il metodo di Cramer.

Ora puoi passare ad altri modi per trovare i coefficienti и .

Discesa gradiente

Innanzitutto, definiamo cos'è un gradiente. In poche parole, il gradiente è un segmento che indica la direzione di massima crescita di una funzione. Per analogia con la scalata di una montagna, dove le pareti in pendenza sono dove si trova la salita più ripida verso la cima della montagna. Sviluppando l'esempio con la montagna, ricordiamo che in realtà abbiamo bisogno della discesa più ripida per raggiungere il più velocemente possibile la pianura, cioè il minimo, il luogo dove la funzione non aumenta né diminuisce. A questo punto la derivata sarà pari a zero. Pertanto non abbiamo bisogno di un gradiente, ma di un antigradiente. Per trovare l'antigradiente devi solo moltiplicare il gradiente per -1 (meno uno).

Prestiamo attenzione al fatto che una funzione può avere diversi minimi, e scendendo in uno di essi utilizzando l'algoritmo proposto di seguito, non saremo in grado di trovare un altro minimo, che potrebbe essere inferiore a quello trovato. Rilassatevi, questa non è una minaccia per noi! Nel nostro caso abbiamo a che fare con un unico minimo, essendo la nostra funzione sul grafico c'è una parabola regolare. E come tutti dovremmo sapere molto bene dal nostro corso di matematica scolastica, una parabola ha un solo minimo.

Dopo aver scoperto perché avevamo bisogno di un gradiente e anche che il gradiente è un segmento, cioè un vettore con determinate coordinate, che sono esattamente gli stessi coefficienti и possiamo implementare la discesa del gradiente.

Prima di iniziare, suggerisco di leggere solo alcune frasi sull’algoritmo di discesa:

• Determiniamo in modo pseudo-casuale le coordinate dei coefficienti и . Nel nostro esempio, determineremo coefficienti prossimi allo zero. Questa è una pratica comune, ma ogni caso può avere la propria pratica.
• Dalle coordinate sottrarre il valore della derivata parziale di 1° ordine nel punto . Quindi, se la derivata è positiva, la funzione aumenta. Pertanto, sottraendo il valore del derivato, ci sposteremo nella direzione opposta alla crescita, cioè nella direzione della discesa. Se la derivata è negativa allora la funzione a questo punto diminuisce e sottraendo il valore della derivata ci muoviamo nella direzione di discesa.
• Eseguiamo un'operazione simile con la coordinata : sottrarre il valore della derivata parziale nel punto .
• Per non saltare oltre il minimo e volare nello spazio profondo, è necessario impostare la dimensione del gradino nella direzione di discesa. In generale si potrebbe scrivere un articolo intero su come impostare correttamente il passo e su come modificarlo durante il processo di discesa in modo da ridurre i costi computazionali. Ma ora abbiamo davanti a noi un compito leggermente diverso e stabiliremo la dimensione del passo utilizzando il metodo scientifico del “poke” o, come si dice nel linguaggio comune, empiricamente.
• Una volta che siamo dalle coordinate indicate и sottraiamo i valori delle derivate, otteniamo nuove coordinate и . Facciamo il passo successivo (sottrazione), già dalle coordinate calcolate. E così il ciclo ricomincia ancora e ancora, finché non viene raggiunta la convergenza richiesta.

Tutto! Ora siamo pronti per andare alla ricerca della gola 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 tuffati fino al fondo della Fossa delle Marianne e lì abbiamo trovato tutti gli stessi valori di coefficiente и , che è esattamente quello che c'era da aspettarselo.

Facciamo un'altra immersione, solo che questa volta il nostro veicolo sottomarino sarà riempito con altre tecnologie, vale a dire una 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 и immutabile.

Diamo un'occhiata a come è cambiato l'errore durante la discesa del gradiente, ovvero come è cambiata la somma delle deviazioni al quadrato ad ogni passaggio.

Codice per tracciare somme di deviazioni quadrate

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

Grafico n. 4 “Somma dei quadrati delle deviazioni durante la discesa del gradiente”

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

Infine, stimiamo la differenza nel tempo di esecuzione del codice:

Codice per determinare il tempo di calcolo della discesa del gradiente

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)

Forse stiamo facendo qualcosa di sbagliato, ma anche in questo caso si tratta di una semplice funzione “scritta in casa” che non utilizza la libreria NumPy supera il tempo di calcolo di una funzione utilizzando la libreria NumPy.

Ma non stiamo fermi, ma ci stiamo muovendo verso lo studio di un altro entusiasmante modo per risolvere la semplice equazione di regressione lineare. Incontrare!

Discesa stocastica del gradiente

Per comprendere rapidamente il principio di funzionamento della discesa del gradiente stocastica, è meglio determinarne le differenze rispetto alla discesa del gradiente ordinaria. Noi, nel caso della discesa del gradiente, nelle equazioni delle derivate di и utilizzato le somme dei valori di tutte le caratteristiche e le risposte vere disponibili nel campione (ovvero le somme di tutti и ). Nella discesa del gradiente stocastico, non utilizzeremo tutti i valori presenti nel campione, ma selezioneremo invece in modo pseudo-casuale il cosiddetto indice del campione e utilizzeremo i suoi valori.

Ad esempio, se si determina che l'indice è il numero 3 (tre), prendiamo i valori и , quindi sostituiamo i valori nelle equazioni derivate e determiniamo nuove coordinate. Quindi, dopo aver determinato le coordinate, determiniamo nuovamente in modo pseudo-casuale l'indice del campione, sostituiamo i valori corrispondenti all'indice nelle equazioni alle derivate parziali e determiniamo le coordinate in un modo nuovo и eccetera. finché la convergenza diventa verde. A prima vista potrebbe non sembrare che possa funzionare affatto, ma è così. È vero che vale la pena notare che l’errore non diminuisce ad ogni passo, ma sicuramente c’è una tendenza.

Quali sono i vantaggi della discesa del gradiente stocastico rispetto a quella convenzionale? Se la dimensione del nostro campione è molto ampia e misurata in decine di migliaia di valori, allora è molto più semplice elaborarne, ad esempio, un migliaio a caso, piuttosto che l’intero campione. È qui che entra in gioco la discesa stocastica del gradiente. Nel nostro caso, ovviamente, non noteremo molta differenza.

Diamo un'occhiata al codice.

Codice per la discesa del 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 sorprendiamo a porci la domanda "Come può essere?" Abbiamo altri valori di coefficiente и . Forse la discesa del gradiente stocastico ha trovato parametri più ottimali per l'equazione? Sfortunatamente no. Basta guardare la somma dei quadrati delle deviazioni e vedere che con nuovi valori dei coefficienti l'errore è maggiore. Non abbiamo fretta di disperare. Costruiamo un grafico della modifica dell'errore.

Codice per tracciare la somma delle deviazioni quadrate nella discesa del gradiente stocastico

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

Grafico n. 5 “Somma dei quadrati delle deviazioni durante la discesa del gradiente stocastico”

Guardando il programma, tutto va a posto e ora sistemeremo tutto.

Allora, cos'è successo? È successo quanto segue. Quando selezioniamo un mese in modo casuale, è per il mese selezionato che il nostro algoritmo cerca di ridurre l'errore nel calcolo delle entrate. Quindi selezioniamo un altro mese e ripetiamo il calcolo, ma riduciamo l'errore per il secondo mese selezionato. Ora ricorda che i primi due mesi si discostano in modo significativo dalla linea della semplice equazione di regressione lineare. Ciò significa che quando viene selezionato uno qualsiasi di questi due mesi, riducendo l'errore di ciascuno di essi, il nostro algoritmo aumenta notevolmente l'errore per l'intero campione. Quindi che si fa? La risposta è semplice: bisogna ridurre il gradino di discesa. Dopotutto, riducendo il passo di discesa, anche l'errore smetterà di “saltare” su e giù. O meglio, l'errore di “salto” non si fermerà, ma non lo farà così velocemente :) Controlliamo.

Codice per eseguire SGD con incrementi più piccoli

# запустим функцию, уменьшив шаг в 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 durante la discesa del gradiente stocastico (80mila passi)”

I coefficienti sono migliorati, ma non sono ancora ideali. Ipoteticamente, questo può essere corretto in questo modo. Selezioniamo, ad esempio, nelle ultime 1000 iterazioni i valori dei coefficienti con cui è stato commesso l'errore minimo. È vero, per questo dovremo annotare anche i valori dei coefficienti stessi. Non lo faremo, ma presteremo attenzione al programma. Sembra liscio e l'errore sembra diminuire in modo uniforme. In realtà, questo non è vero. Diamo un'occhiata alle prime 1000 iterazioni e confrontiamole con le ultime.

Codice per grafico SGD (primi 1000 passi)

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

Grafico n. 7 “Somma dei quadrati delle deviazioni SGD (primi 1000 passi)”

Grafico n. 8 “Somma dei quadrati delle deviazioni SGD (ultimi 1000 passi)”

All'inizio della discesa osserviamo una diminuzione dell'errore abbastanza uniforme e marcata. Nelle ultime iterazioni vediamo che l'errore gira intorno al valore di 1,475 e in alcuni momenti arriva anche a questo valore ottimale, ma poi sale ancora... Ripeto, potete annotare i valori del coefficienti и , quindi seleziona quelli per i quali l'errore è minimo. Abbiamo però avuto un problema più serio: abbiamo dovuto fare 80mila passi (vedi codice) per ottenere valori prossimi a quelli ottimali. E questo già contraddice l’idea di risparmiare tempo di calcolo con la discesa stocastica del gradiente rispetto alla discesa del gradiente. Cosa si può correggere e migliorare? Non è difficile notare che nelle prime iterazioni stiamo scendendo con sicurezza e, quindi, dovremmo lasciare un ampio passo nelle prime iterazioni e ridurre il passo man mano che andiamo avanti. Non lo faremo in questo articolo: è già troppo lungo. Chi lo desidera può pensare da solo a come farlo, non è difficile :)

Ora eseguiamo la discesa del gradiente stocastico utilizzando la libreria NumPy (e non inciampiamo nelle pietre che abbiamo individuato prima)

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 si sono rivelati quasi gli stessi di quando si scendeva senza utilizzare NumPy. Tuttavia, questo è logico.

Scopriamo quanto tempo ci hanno impiegato le discese a gradiente stocastico.

Codice per determinare il tempo di calcolo SGD (80mila passaggi)

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)

Più ci si addentra nella foresta, più scure sono le nuvole: ancora una volta, la formula “autoscritta” mostra il risultato migliore. Tutto ciò suggerisce che devono esserci modi ancora più sottili per utilizzare la biblioteca NumPy, che velocizzano davvero le operazioni di calcolo. In questo articolo non ne parleremo. Ci sarà qualcosa a cui pensare nel tuo tempo libero :)

Riassumiamo

Prima di riassumere, vorrei rispondere ad una domanda che molto probabilmente è nata dal nostro caro lettore. Perché, infatti, tale “tortura” con le discese, perché dobbiamo camminare su e giù per la montagna (soprattutto giù) per trovare la preziosa pianura, se abbiamo tra le mani un dispositivo così potente e semplice, nel forma di una soluzione analitica, che ci teletrasporta istantaneamente nel posto giusto?

La risposta a questa domanda si trova in superficie. Ora abbiamo esaminato un esempio molto semplice, in cui la risposta vera è dipende da un segno . Non lo vedi spesso nella vita, quindi immaginiamo di avere 2, 30, 50 o più segni. Aggiungiamo a questo migliaia, o anche decine di migliaia di valori per ciascun attributo. In questo caso, la soluzione analitica potrebbe non resistere al test e fallire. A sua volta, la discesa del gradiente e le sue variazioni ci avvicineranno lentamente ma inesorabilmente all'obiettivo: il minimo della funzione. E non preoccuparti della velocità: probabilmente esamineremo i modi che ci consentiranno di impostare e regolare la lunghezza del passo (ovvero la velocità).

E ora il vero e proprio breve riassunto.

In primo luogo, spero che il materiale presentato nell'articolo aiuti i “data scientist” principianti a capire come risolvere semplici (e non solo) equazioni di regressione lineare.

In secondo luogo, abbiamo esaminato diversi modi per risolvere l’equazione. Ora, a seconda della situazione, possiamo scegliere quello più adatto a risolvere il problema.

In terzo luogo, abbiamo visto la potenza delle impostazioni aggiuntive, ovvero la lunghezza del gradino di discesa del gradiente. Questo parametro non può essere trascurato. Come notato sopra, per ridurre il costo dei calcoli, la lunghezza del passo dovrebbe essere modificata durante la discesa.

In quarto luogo, nel nostro caso, le funzioni “scritte in casa” hanno mostrato i migliori risultati temporali per i calcoli. Ciò è probabilmente dovuto all'uso non proprio professionale delle capacità della biblioteca NumPy. Comunque sia, la seguente conclusione suggerisce da sola. Da un lato, a volte vale la pena mettere in discussione le opinioni consolidate e, dall'altro, non sempre vale la pena complicare tutto: al contrario, a volte un modo più semplice per risolvere un problema è più efficace. E poiché il nostro obiettivo era analizzare tre approcci per risolvere una semplice equazione di regressione lineare, l'uso di funzioni "autoprodotte" per noi era abbastanza.

Letteratura (o qualcosa del genere)

1. Regressione lineare

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

2. Metodo dei minimi quadrati

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Derivato

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. gradiente

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Discesa del gradiente

habr.com/en/post/471458

habr.com/en/post/307312

artemarakcheev.com//2017-12-31/linear_regression

6. Libreria 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