Resolució de l'equació de regressió lineal simple

L'article analitza diverses maneres de determinar l'equació matemàtica d'una recta de regressió simple (aparellada).

Tots els mètodes per resoldre l'equació que es discuteixen aquí es basen en el mètode dels mínims quadrats. Denotem els mètodes de la següent manera:

  • Solució analítica
  • Descens Gradient
  • Descens del gradient estocàstic

Per a cada mètode de resolució de l'equació d'una línia recta, l'article ofereix diverses funcions, que es divideixen principalment en les que s'escriuen sense utilitzar la biblioteca. numpy i els que s'utilitzen per a càlculs numpy. Es creu que l'ús hàbil numpy reduirà els costos informàtics.

Tot el codi que es dóna a l'article està escrit en l'idioma python 2.7 utilitzant Jupyter Llibreta. El codi font i el fitxer amb dades de mostra es publiquen a Github

L'article està més adreçat tant a principiants com a aquells que ja han començat a dominar l'estudi d'un apartat molt ampli de la intel·ligència artificial: l'aprenentatge automàtic.

Per il·lustrar el material, utilitzem un exemple molt senzill.

Condicions d'exemple

Tenim cinc valors que caracteritzen la dependència Y d' X (Taula núm. 1):

Taula núm. 1 "Condicions d'exemple"

Resolució de l'equació de regressió lineal simple

Suposarem que els valors Resolució de l'equació de regressió lineal simple és el mes de l'any, i Resolució de l'equació de regressió lineal simple - ingressos aquest mes. En altres paraules, els ingressos depenen del mes de l'any, i Resolució de l'equació de regressió lineal simple - l'únic signe del qual depenen els ingressos.

L'exemple és així, tant des del punt de vista de la dependència condicional dels ingressos del mes de l'any, com des del punt de vista del nombre de valors, n'hi ha molt pocs. Tanmateix, aquesta simplificació permetrà, com diuen, explicar, no sempre amb facilitat, el material que assimilen els principiants. I també la senzillesa dels números permetrà a qui vulgui resoldre l'exemple en paper sense costos laborals importants.

Suposem que la dependència donada a l'exemple es pot aproximar força bé mitjançant l'equació matemàtica d'una recta de regressió simple (aparellada) de la forma:

Resolució de l'equació de regressió lineal simple

on Resolució de l'equació de regressió lineal simple és el mes en què es van rebre els ingressos, Resolució de l'equació de regressió lineal simple - ingressos corresponents al mes, Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple són els coeficients de regressió de la recta estimada.

Tingueu en compte que el coeficient Resolució de l'equació de regressió lineal simple sovint anomenat pendent o gradient de la línia estimada; representa la quantitat en què el Resolució de l'equació de regressió lineal simple quan canvia Resolució de l'equació de regressió lineal simple.

Òbviament, la nostra tasca a l'exemple és seleccionar aquests coeficients a l'equació Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, en què les desviacions dels nostres valors d'ingressos calculats per mes de les respostes reals, és a dir. els valors presentats a la mostra seran mínims.

Mètode dels mínims quadrats

Segons el mètode dels mínims quadrats, la desviació s'ha de calcular al quadrat. Aquesta tècnica permet evitar la cancel·lació mútua de desviacions si tenen signes oposats. Per exemple, si en un cas, la desviació és +5 (més cinc), i en l'altre -5 (menys cinc), aleshores la suma de les desviacions s'anul·larà mútuament i ascendirà a 0 (zero). És possible no quadrar la desviació, sinó utilitzar la propietat del mòdul i llavors totes les desviacions seran positives i s'acumularan. No ens detenem en aquest punt en detall, sinó que simplement indicarem que, per a la comoditat dels càlculs, s'acostuma a quadrar la desviació.

Així és la fórmula amb la qual determinarem la mínima suma de desviacions al quadrat (errors):

Resolució de l'equació de regressió lineal simple

on Resolució de l'equació de regressió lineal simple és una funció de l'aproximació de respostes veritables (és a dir, els ingressos que hem calculat),

Resolució de l'equació de regressió lineal simple són les respostes reals (ingressos proporcionats a la mostra),

Resolució de l'equació de regressió lineal simple és l'índex mostral (número del mes en què es determina la desviació)

Diferenciarem la funció, definim les equacions diferencials parcials i estiguem preparats per passar a la solució analítica. Però primer, fem una petita excursió sobre què és la diferenciació i recordem el significat geomètric de la derivada.

Diferenciació

La diferenciació és l'operació de trobar la derivada d'una funció.

Per a què serveix la derivada? La derivada d'una funció caracteritza la velocitat de canvi de la funció i ens indica la seva direcció. Si la derivada en un punt donat és positiva, aleshores la funció augmenta; en cas contrari, la funció disminueix. I com més gran sigui el valor de la derivada absoluta, més gran serà la taxa de canvi dels valors de la funció, així com més pendent és la gràfica de la funció.

Per exemple, en les condicions d'un sistema de coordenades cartesianes, el valor de la derivada en el punt M(0,0) és igual a + 25 significa que en un punt donat, quan es desplaça el valor Resolució de l'equació de regressió lineal simple a la dreta per una unitat convencional, valor Resolució de l'equació de regressió lineal simple augmenta en 25 unitats convencionals. Al gràfic sembla un augment força pronunciat dels valors Resolució de l'equació de regressió lineal simple des d'un punt donat.

Un altre exemple. El valor de la derivada és igual -0,1 vol dir que quan es desplaça Resolució de l'equació de regressió lineal simple per unitat convencional, valor Resolució de l'equació de regressió lineal simple disminueix només en 0,1 unitat convencional. Al mateix temps, a la gràfica de la funció, podem observar un pendent descendent amb prou feines perceptible. Fent una analogia amb una muntanya, és com si baixéssim molt a poc a poc un pendent suau des d'una muntanya, a diferència de l'exemple anterior, on havíem de pujar cims molt pronunciats :)

Així, després de diferenciar la funció Resolució de l'equació de regressió lineal simple per probabilitats Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, definim equacions en derivades parcials de 1r ordre. Després de determinar les equacions, rebrem un sistema de dues equacions, resolent les quals podrem seleccionar aquests valors dels coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, per als quals els valors de les derivades corresponents en determinats punts canvien en una quantitat molt, molt petita, i en el cas d'una solució analítica no canvien en absolut. És a dir, la funció d'error als coeficients trobats arribarà a un mínim, ja que els valors de les derivades parcials en aquests punts seran iguals a zero.

Així, segons les regles de diferenciació, l'equació de derivada parcial de 1r ordre respecte al coeficient Resolució de l'equació de regressió lineal simple prendrà la forma:

Resolució de l'equació de regressió lineal simple

Equació de derivada parcial de 1r ordre respecte a Resolució de l'equació de regressió lineal simple prendrà la forma:

Resolució de l'equació de regressió lineal simple

Com a resultat, vam rebre un sistema d'equacions que té una solució analítica bastant senzilla:

començar{equació*}
començar{casos}
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
acabar {casos}
final{equació*}

Abans de resoldre l'equació, precarreguem, comprovem que la càrrega sigui correcta i formem les dades.

Càrrega i format de dades

Cal tenir en compte que a causa del fet que per a la solució analítica, i posteriorment per a la baixada de gradient i gradient estocàstic, utilitzarem el codi en dues variants: utilitzant la biblioteca numpy i sense utilitzar-lo, llavors necessitarem un format de dades adequat (vegeu codi).

Codi de càrrega i processament de dades

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

Visualització

Ara, després d'haver, en primer lloc, carregat les dades, en segon lloc, comprovat la correcció de la càrrega i finalment format les dades, realitzarem la primera visualització. El mètode utilitzat sovint per a això és trama de parella biblioteca Nascut al mar. En el nostre exemple, a causa del nombre limitat, no té sentit utilitzar la biblioteca Nascut al mar. Farem servir la biblioteca habitual matplotlib i només mira el diagrama de dispersió.

Codi de diagrama de dispersió

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

Gràfic núm. 1 “Dependència dels ingressos del mes de l'any”

Resolució de l'equació de regressió lineal simple

Solució analítica

Utilitzem les eines més habituals pitó i resol el sistema d'equacions:

començar{equació*}
començar{casos}
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
acabar {casos}
final{equació*}

Segons la regla de Cramer trobarem el determinant general, així com els determinants per Resolució de l'equació de regressió lineal simple i per Resolució de l'equació de regressió lineal simple, després de la qual cosa, dividint el determinant per Resolució de l'equació de regressió lineal simple al determinant general: troba el coeficient Resolució de l'equació de regressió lineal simple, de la mateixa manera trobem el coeficient Resolució de l'equació de regressió lineal simple.

Codi de solució analítica

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

Això és el que tenim:

Resolució de l'equació de regressió lineal simple

Així doncs, s'han trobat els valors dels coeficients, s'ha establert la suma de les desviacions al quadrat. Tracem una línia recta a l'histograma de dispersió d'acord amb els coeficients trobats.

Codi de línia de regressió

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

Gràfic núm. 2 “Respostes correctes i calculades”

Resolució de l'equació de regressió lineal simple

Podeu consultar el gràfic de desviació de cada mes. En el nostre cas, no en treurem cap valor pràctic significatiu, però satisfarem la nostra curiositat sobre com caracteritza l'equació de regressió lineal simple la dependència dels ingressos del mes de l'any.

Codi gràfic de desviació

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

Gràfic núm. 3 "Desviacions, %"

Resolució de l'equació de regressió lineal simple

No és perfecte, però hem completat la nostra tasca.

Escrivim una funció que, per determinar els coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple utilitza la biblioteca numpy, més precisament, escriurem dues funcions: una utilitzant una matriu pseudoinversa (no recomanable a la pràctica, ja que el procés és computacionalment complex i inestable), l'altra utilitzant una equació matricial.

Codi de solució analítica (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

Comparem el temps dedicat a la determinació dels coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, d'acord amb els 3 mètodes presentats.

Codi per calcular el temps de càlcul

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)

Resolució de l'equació de regressió lineal simple

Amb una petita quantitat de dades, apareix una funció "autoescrita", que troba els coeficients mitjançant el mètode de Cramer.

Ara podeu passar a altres maneres de trobar coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple.

Descens Gradient

Primer, anem a definir què és un gradient. En poques paraules, el gradient és un segment que indica la direcció del creixement màxim d'una funció. Per analogia amb escalar una muntanya, on el desnivell s'enfronta és on hi ha la pujada més pronunciada fins al cim de la muntanya. Desenvolupant l'exemple amb la muntanya, recordem que, de fet, necessitem la baixada més pronunciada per arribar el més ràpid possible a la terra baixa, és a dir, el mínim, el lloc on la funció no augmenta ni disminueix. En aquest punt la derivada serà igual a zero. Per tant, no necessitem un gradient, sinó un antigradient. Per trobar l'antigradient només heu de multiplicar el gradient per -1 (menys un).

Fixem-nos en el fet que una funció pot tenir diversos mínims, i havent baixat a un d'ells mitjançant l'algorisme proposat a continuació, no podrem trobar un altre mínim, que pot ser inferior al trobat. Relaxem-nos, això no és una amenaça per a nosaltres! En el nostre cas estem davant d'un únic mínim, ja que la nostra funció Resolució de l'equació de regressió lineal simple a la gràfica hi ha una paràbola regular. I com tots hauríem de saber molt bé pel curs de matemàtiques de l'escola, una paràbola només té un mínim.

Després vam descobrir per què necessitàvem un gradient, i també que el gradient és un segment, és a dir, un vector amb coordenades donades, que són precisament els mateixos coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple podem implementar el descens de gradients.

Abans de començar, suggereixo llegir només unes quantes frases sobre l'algoritme de descens:

  • Determinem de manera pseudoaleatòria les coordenades dels coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple. En el nostre exemple, determinarem coeficients propers a zero. Aquesta és una pràctica habitual, però cada cas pot tenir la seva pròpia pràctica.
  • Des de la coordenada Resolució de l'equació de regressió lineal simple resta el valor de la derivada parcial de 1r ordre en el punt Resolució de l'equació de regressió lineal simple. Per tant, si la derivada és positiva, la funció augmenta. Per tant, restant el valor de la derivada, ens mourem en el sentit contrari del creixement, és a dir, en el sentit del descens. Si la derivada és negativa, aleshores la funció en aquest punt disminueix i en restar el valor de la derivada ens movem en la direcció del descens.
  • Realitzem una operació similar amb la coordenada Resolució de l'equació de regressió lineal simple: resta el valor de la derivada parcial en el punt Resolució de l'equació de regressió lineal simple.
  • Per no saltar per sobre del mínim i volar a l'espai profund, cal establir la mida del pas en la direcció del descens. En general, podríeu escriure un article sencer sobre com configurar el pas correctament i com canviar-lo durant el procés de descens per tal de reduir els costos computacionals. Però ara tenim una tasca una mica diferent per davant, i establirem la mida del pas mitjançant el mètode científic de "poke" o, com diuen en el llenguatge comú, empíricament.
  • Un cop som de les coordenades donades Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple restem els valors de les derivades, obtenim noves coordenades Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple. Fem el següent pas (resta), ja a partir de les coordenades calculades. I així el cicle comença una i altra vegada, fins que s'aconsegueix la convergència requerida.

Tot! Ara estem preparats per anar a la recerca del congost més profund de la Fossa de les Mariannes. Comencem.

Codi per a la baixada de desnivells

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

Resolució de l'equació de regressió lineal simple

Vam submergir-nos fins al fons de la fossa de les Mariannes i allà vam trobar tots els mateixos valors de coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, que és exactament el que s'esperava.

Fem una altra immersió, només que aquesta vegada, el nostre vehicle de fons s'omplirà d'altres tecnologies, és a dir, una biblioteca numpy.

Codi per a la baixada del gradient (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

Resolució de l'equació de regressió lineal simple
Valors dels coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple immutable.

Vegem com va canviar l'error durant el descens del gradient, és a dir, com va canviar la suma de les desviacions al quadrat amb cada pas.

Codi per representar sumes de desviacions al quadrat

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

Gràfic núm. 4 “Suma de desviacions al quadrat durant el descens del gradient”

Resolució de l'equació de regressió lineal simple

A la gràfica veiem que a cada pas l'error disminueix, i després d'un cert nombre d'iteracions observem una línia gairebé horitzontal.

Finalment, estimem la diferència en el temps d'execució del codi:

Codi per determinar el temps de càlcul de descens del gradient

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)

Resolució de l'equació de regressió lineal simple

Potser estem fent alguna cosa malament, però de nou és una funció simple "escrita a casa" que no utilitza la biblioteca numpy supera el temps de càlcul d'una funció utilitzant la biblioteca numpy.

Però no ens quedem quiets, sinó que estem avançant cap a estudiar una altra manera interessant de resoldre l'equació de regressió lineal simple. Trobar-se!

Descens del gradient estocàstic

Per entendre ràpidament el principi de funcionament del descens del gradient estocàstic, és millor determinar les seves diferències amb el descens del gradient normal. Nosaltres, en el cas del descens del gradient, en les equacions de derivades de Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple va utilitzar les sumes dels valors de totes les característiques i les respostes veritables disponibles a la mostra (és a dir, les sumes de totes les Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple). En el descens de gradient estocàstic, no utilitzarem tots els valors presents a la mostra, sinó que seleccionem pseudoaleatòriament l'anomenat índex de mostra i utilitzarem els seus valors.

Per exemple, si es determina que l'índex és el número 3 (tres), prenem els valors Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple, llavors substituïm els valors en les equacions derivades i determinem noves coordenades. Aleshores, havent determinat les coordenades, tornem a determinar pseudoaleatòriament l'índex de mostra, substituïm els valors corresponents a l'índex per les equacions en derivades parcials i determinem les coordenades d'una manera nova. Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple etc. fins que la convergència es torna verda. A primera vista, pot semblar que això no podria funcionar en absolut, però sí. És cert que val la pena assenyalar que l'error no disminueix a cada pas, però sens dubte hi ha una tendència.

Quins són els avantatges del descens del gradient estocàstic respecte al convencional? Si la mida de la nostra mostra és molt gran i es mesura en desenes de milers de valors, llavors és molt més fàcil processar-ne, per exemple, un miler aleatori, en lloc de la mostra sencera. Aquí és on entra en joc el descens del gradient estocàstic. En el nostre cas, és clar, no notarem gaire diferència.

Mirem el codi.

Codi per al descens de gradient estocàstic

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

Resolució de l'equació de regressió lineal simple

Mirem atentament els coeficients i ens sorprenem fent la pregunta "Com pot ser això?" Tenim altres valors de coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simple. Potser el descens del gradient estocàstic ha trobat paràmetres més òptims per a l'equació? Lamentablement no. N'hi ha prou amb mirar la suma de les desviacions al quadrat i veure que amb nous valors dels coeficients, l'error és més gran. No tenim pressa per desesperar-nos. Construïm un gràfic del canvi d'error.

Codi per representar la suma de desviacions al quadrat en descens del gradient estocàstic

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

Gràfic núm. 5 “Suma de desviacions al quadrat durant el descens del gradient estocàstic”

Resolució de l'equació de regressió lineal simple

Mirant l'horari, tot s'ajusta i ara ho arreglarem tot.

Llavors què va passar? Va passar el següent. Quan seleccionem un mes aleatòriament, és per al mes seleccionat que el nostre algorisme pretén reduir l'error en el càlcul dels ingressos. Aleshores seleccionem un altre mes i repetim el càlcul, però reduïm l'error del segon mes seleccionat. Ara recordeu que els dos primers mesos es desvien significativament de la línia de l'equació de regressió lineal simple. Això vol dir que quan es selecciona qualsevol d'aquests dos mesos, en reduir l'error de cadascun d'ells, el nostre algorisme augmenta seriosament l'error per a tota la mostra. Llavors, què fer? La resposta és senzilla: cal reduir el pas de baixada. Després de tot, en reduir el pas de baixada, l'error també deixarà de "saltar" amunt i avall. O millor dit, l'error de "salt" no s'aturarà, però no ho farà tan ràpid :) Comprovem.

Codi per executar SGD amb increments més petits

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

Resolució de l'equació de regressió lineal simple

Gràfic núm. 6 “Suma de desviacions al quadrat durant el descens del gradient estocàstic (80 mil passos)”

Resolució de l'equació de regressió lineal simple

Els coeficients han millorat, però encara no són ideals. Hipotèticament, això es pot corregir d'aquesta manera. Seleccionem, per exemple, en les últimes 1000 iteracions els valors dels coeficients amb els quals es va cometre l'error mínim. És cert que per a això també haurem d'anotar els valors dels propis coeficients. No ho farem, sinó que pararem atenció a l'horari. Sembla suau i l'error sembla que disminueix uniformement. De fet, això no és cert. Mirem les primeres 1000 iteracions i comparem-les amb les darreres.

Codi per al gràfic SGD (primers 1000 passos)

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

Gràfic núm. 7 "Suma de desviacions al quadrat SGD (primers 1000 passos)"

Resolució de l'equació de regressió lineal simple

Gràfic núm. 8 "Suma de desviacions al quadrat SGD (darrers 1000 passos)"

Resolució de l'equació de regressió lineal simple

Al principi de la baixada, observem una disminució de l'error força uniforme i pronunciada. En les últimes iteracions, veiem que l'error gira al voltant del valor d'1,475 i en alguns moments fins i tot iguala aquest valor òptim, però després encara puja... Repeteixo, podeu anotar els valors de la coeficients Resolució de l'equació de regressió lineal simple и Resolució de l'equació de regressió lineal simplei, a continuació, seleccioneu aquells per als quals l'error és mínim. Tanmateix, teníem un problema més greu: vam haver de fer 80 mil passos (vegeu codi) per aconseguir valors propers a l'òptim. I això ja contradiu la idea d'estalviar temps de càlcul amb el descens del gradient estocàstic en relació amb el descens del gradient. Què es pot corregir i millorar? No és difícil notar que en les primeres iteracions anem baixant amb confiança i, per tant, hauríem de deixar un gran pas en les primeres iteracions i reduir el pas a mesura que avancem. No ho farem en aquest article, ja és massa llarg. Qui ho desitgi pot pensar per si mateix com fer-ho, no és difícil :)

Ara realitzem un descens de gradient estocàstic mitjançant la biblioteca numpy (i no ensopeguem amb les pedres que hem identificat abans)

Codi per al descens del gradient estocàstic (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

Resolució de l'equació de regressió lineal simple

Els valors van resultar ser gairebé els mateixos que en baixar sense utilitzar numpy. Tanmateix, això és lògic.

Descobrim quant de temps ens van trigar els descensos en gradient estocàstic.

Codi per determinar el temps de càlcul de SGD (80 mil passos)

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)

Resolució de l'equació de regressió lineal simple

Com més endins al bosc, més foscos són els núvols: de nou, la fórmula "autoescrita" mostra el millor resultat. Tot això fa pensar que hi ha d'haver maneres encara més subtils d'utilitzar la biblioteca numpy, que realment acceleren les operacions de càlcul. En aquest article no aprendrem sobre ells. Hi haurà alguna cosa en què pensar en el teu temps lliure :)

Resumim

Abans de resumir, m'agradaria respondre una pregunta que probablement va sorgir del nostre estimat lector. Per què, de fet, tanta "tortura" amb baixades, per què necessitem caminar amunt i avall de la muntanya (principalment avall) per trobar la terra baixa atresorada, si tenim a les nostres mans un aparell tan potent i senzill, en el forma d'una solució analítica, que ens teletransporta instantàniament al lloc correcte?

La resposta a aquesta pregunta es troba a la superfície. Ara hem mirat un exemple molt senzill, en el qual la resposta real és Resolució de l'equació de regressió lineal simple depèn d'un signe Resolució de l'equació de regressió lineal simple. Això no ho veus sovint a la vida, així que imaginem que tenim 2, 30, 50 o més signes. Afegim a això milers, o fins i tot desenes de milers de valors per a cada atribut. En aquest cas, la solució analítica pot no suportar la prova i fallar. Al seu torn, el descens de gradients i les seves variacions ens acostaran a poc a poc a l'objectiu: el mínim de la funció. I no us preocupeu per la velocitat: probablement mirarem maneres que ens permetin establir i regular la longitud del pas (és a dir, la velocitat).

I ara el breu resum real.

En primer lloc, espero que el material presentat a l'article ajudi els "científics de dades" que inicien a entendre com resoldre equacions de regressió lineal senzilles (i no només).

En segon lloc, vam mirar diverses maneres de resoldre l'equació. Ara, segons la situació, podem escollir la que més s'adapti per resoldre el problema.

En tercer lloc, vam veure el poder de la configuració addicional, és a dir, la longitud del pas de descens del gradient. Aquest paràmetre no es pot descuidar. Com s'ha indicat anteriorment, per reduir el cost dels càlculs, la longitud del pas s'ha de canviar durant el descens.

En quart lloc, en el nostre cas, les funcions "autoescrites" van mostrar els millors resultats de temps per als càlculs. Això probablement es deu a l'ús no més professional de les capacitats de la biblioteca numpy. Però sigui com sigui, es suggereix la següent conclusió. D'una banda, de vegades val la pena qüestionar les opinions consolidades i, de l'altra, no sempre val la pena complicar-ho tot, al contrari, de vegades és més eficaç una manera més senzilla de resoldre un problema. I com que el nostre objectiu era analitzar tres enfocaments per resoldre una equació de regressió lineal senzilla, l'ús de funcions "autoescrites" ens va ser suficient.

Literatura (o alguna cosa així)

1. Regressió lineal

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

2. Mètode dels mínims quadrats

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Derivada

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4 Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Descens en desnivell

habr.com/en/post/471458

habr.com/en/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

Font: www.habr.com

Afegeix comentari