Sadə xətalı reqressiya tənliklərini həll edirik

Məqalədə sadə (cüt) reqressiya xəttinin riyazi tənliyini müəyyən etmənin bir neçə yolu müzakirə olunur.

Burada müzakirə edilən bütün həll üsulları, ən kiçik kvadratlar metoduna əsaslanır. Üsullar aşağıdakı kimi adlandırılır:

• Analitik həll
• Qradient enməsi
• Stoxastik qradient enməsi

Xəttin tənliyinin həlli üçün hər bir üsul üzrə, məqalədə əsasən kitabxana istifadə etmədən yazılmış müxtəlif funksiyalar təqdim olunur NumPy və hesablama aparmaq üçün istifadə edənlər NumPy. Bacarıqlı istifadə edildikdə, bu NumPy hesablama xərclərini azaltmağa kömək edəcəkdir.

Məqalədə təqdim olunan bütün kodlar python 2.7 dilində yazılmışdır Jupyter Notebook. Mənbə kodu və nümunə məlumatları faylı Github-da

Məqalə həm başlayana, həm də artıq süni intellektin - maşın öyrənmənin geniş sahəsini öyrənməyə başlamış olanlara istiqamətlənmişdir.

Materialı təsvir etmək üçün çox sadə bir nümunə istifadə edirik.

Nümunənin şərtləri

Bizim beş dəyərimiz var, bunlar Y aylıq X (Cədvəl №1):

Cədvəl №1 "Nümunənin şərtləri"

Qoy düşünək ki, dəyərlər ildəki ayı, bu ayda gəlir. Başqa sözlə, gəlir ildəki aya bağlıdır və gəlirə təsir edən yeganə əlamətdir.

Nümunə olduqca sadədir, həm gəlirin il ayı ilə əlaqəsi, həm də dəyər sayı — çox azdır. Lakin bu sadələşdirmə, yeni başlayanların anlayıb qavramaq üçün çətinlik çəkmədən izah etməyə imkan verəcək. Həmçinin, sadə rəqəmlər istəyəni "kağızda" nümunəni həll etmək üçün ağır iş yükü tələb etmədən icra etməyə imkan verəcək.

Tutaq ki, nümunədə təqdim olunan asılılıq, sadə (cüt) reqressiya xəttinin tənliyi ilə kifayət qədər yaxşı aproksimasiya edilə bilər:

burada gəlirin əldə edildiyi ayı, ayına uyğun gəlir, və reqressiya xəttinin koeffisientləri.

Qeyd edək ki, koeffisient tez-tez qiymətləndirilmiş xəttin meyl koeffisienti olaraq adlandırılır; dəyişdikdə dəyişəcək deyil .

. Aydındır ki, nümunədə bizim vəzifəmiz tənlikdə belə koeffisientləri tapmaqdır ki, aylar üzrə hesablama dəyərlərimizin əsl cavablardan (yəni, nümunədə təqdim olunan dəyərlərdən) sapmaları minimal olsun. və Ən kiçik kvadratlar metodu

Minimum kvadratlar metodu

Ən kiçik kvadratlar metoduna uyğun olaraq, sapmanın hesablanması üçün onu kvadratlaşdırmaq lazımdır. Belə bir yanaşma, əks işarəyə malik olan sapmaların qarşılıqlı ləğv edilməsinin qarşısını alır. Məsələn, bir halda sapma +5 (beş), digər halda isə -5 (-müxalif beş) olur, beləliklə sapmaların cəmi bir-birini ləğv edir və 0 (sıfır) olur. Sapmanı kvadratlaşdırmadan istifadə etməyə də bilərik, modulun xüsusiyyətindən yararlana bilərik və bütün sapmalar pozitiv olacaq və yığılacaq. Bu nöqtədə dərinləşməyəcəyik, sadəcə bildirmək istəyirik ki, hesablamalar üçün rahatlıq məqsədilə sapmanın kvadratlaşdırılması qəbul edilmişdir.

Budur, biz sapmaların kvadratlarının ən kiçik cəmini müəyyənləşdirmək üçün istifadə edəcəyimiz formul:

burada — bu, həqiqi cavabların (yəni bizim hesabladığımız gəlir) aproksimasiya funksiyasıdır,

— bu, həqiqi cavablar (nümunədə təqdim olunan gəlir)dır,

— bu, nümunə indeksi (sapmanın müəyyən edildiyi ayın nömrəsi)dir.

Funksiyanı diferensiallaşdıracağıq, qismən törəmələrin tənliklərini müəyyənləşdirəcəyik və analitik həllə keçməyə hazır olacağıq. Lakin əvvəlcə diferensiallaşdırma nədir və törəmənin geometrik mənasını xatırlamaq üçün kiçik bir ekskursiya keçirək.

Diferensiallaşdırma

Diferensiallaşdırma, funksiyanın törəməsini tapma əməliyyatıdır.

Törəmə nə üçündür? Funksiyanın törəməsi onun dəyişmə sürətini xarakterizə edir və bizə istiqamətini göstərir. Əgər verilmiş nöqtədə törəmə müsbət idisə, funksiya artır; əks halda — funksiya azalır. Törəmənin qiyməti nə qədər çoxdursa, funksiya dəyərlərinin dəyişmə sürəti də o qədər yüksəkdir, həmçinin funksiyanın qrafikinin meyli də daha dikdir.

Məsələn, Dekart koordinat sistemində, M(0,0) nöqtəsində törəmənin dəyəri +25 savaşı göstərir ki, verilmiş nöqtədə, dəyəri şərti vahid sağa dəyişdikdə, dəyəri 25 şərti vahid artır. Qrafikdə bu, verilmiş nöqtədən kifayət qədər dik yüksəlmə bucağı kimi görünür. Başqa bir nümunə. Törəmənin dəyəri

bizə göstərir ki, dəyəri -0,1 bir şərti vahid dəyişdikdə, dəyəri yalnız 0.1 şərti vahid azalır. Bu halda, funksiyanın qrafikində, biz azca görünən aşağı meyli görə bilərik. Dağ ilə analoji quraraq, görünür ki, biz çox yavaşca enirik, düz bir meyl altında, əks tərəfdən isə əvvəlki nümunədə olduğu kimi, çox dik zirvələr götürməliydik :) Beləliklə, funksiyanı diferensiallaşdırdıqdan sonra

koeffisientlər üzrə. koeffisientlər üzrə və , biz 1-ci dərəcəli qismən diferensial tənlikləri müəyyən edəcəyik. Tənlikləri müəyyən etdikdən sonra iki tənlikdən ibarət sistem əldə edəcəyik, bunu həll etməklə biz müvafiq koeffisiyentlərin dəyərlərini tapacaq. və , bu dəyərlər müvafiq diferensialların təyin olunmuş nöqtələrdə çox, çox kiçik dəyişikliyə məruz qalır və analitik həll ilə heç bir dəyişikliyə məruz qalmır. Başqa sözlə, tapılmış koeffisiyentlərdə xətanın funksiyası minimuma çatır, çünki bu nöqtələrdəki qismən diferensialların dəyərləri sıfıra bərabər olacaq.

Beləliklə, diferensiallaşdırma qaydalarına əsasən, koeffisiyentə görə 1-ci dərəcəli qismən diferensial tənliyi şəklini alacaq:

1-ci dərəcəli qismən diferensial tənliyi şəklini alacaq:

Nəticədə sadə analitik həlli olan tənliklər sistemini əldə etdik:

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

Tənliyi həll etməzdən əvvəl, yükləməyi, yükləmənin düzgünlüyünü yoxlamağı və məlumatları formatlaşdırmağı unutmayaq.

Məlumatların yüklənməsi və formatlaşdırılması

Qeyd etmək lazımdır ki, analitik həll üçün və daha sonra gradien və stoxastik gradien enməsi üçün kodu iki variantda tətbiq edəcəyimiz üçün: birincisi, kitabxana istifadə etməklə NumPy və ikincisi, kitabxana olmadan, bizə müvafiq məlumatların formatlaşdırılması tələb olunacaq (baxın: kod).

Məlumatların yüklənməsi və işlənməsi kodu

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

Vizualizasiya

İndi, birincisi, məlumatları yüklədikdən, ikincisi, yükləmənin düzgüyn olduğunu yoxladıqdan və nəhayət məlumatları formatlaşdırdıqdan sonra, ilk vizualizasiyanı keçirəcəyik. Bunun üçün tez-tez pairplot metodundan istifadə olunur Seaborn. Bizim nümunəmizdə, ədədlərin məhdudiyyətinə görə kitabxananı tətbiq etməyin bir mənası yoxdur. SeabornSadə kitabxana Matplotlib istifadə edəcəyik və sadəcə dispersiya diagramına baxacağıq.

Dispensiya diaqramı kodu

print 'Qrafik №1 "İlin aylarına görə gəlirin asılılığı"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Aylar$', size=16)
plt.ylabel('$Gəlirlər$', size=16)
plt.show()

Qrafik №1 «İlin aylarına görə gəlirin asılılığı»

Analitik həll

Ən sadə alətlərdən istifadə edərək python sistemini həll edəcəyik:

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

Kramerin qaydasına görə ümumi determinantı, həmçinin və üzrə determinantları tapacağıq, sonra determinantı ümumi determinantla bölərək koeffisiyentini tapacağıq , eynilə koeffisiyentini .

tapacağıq.

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

Analitik həllin kodu

Beləliklə, koeffisiyentlərin dəyərləri tapılıb, anomaliya cəmləri müəyyən edilib. Tapılmış koeffisiyentlərə uyğun olarak, dispersiya histogramında düz xətt çəkəcəyik.

Rejresiya xəttinin kodu

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

Qrafik №2 «Düzgün və hesablanmış cavablar»

Hər ay üçün dəyişiklikləri göstərən qrafiki nəzərdən keçirə bilərsiniz. Bizim halda, bunun praktiki əhəmiyyəti yoxdur, lakin ilin ayına görə gəlirin asılılığını sadə xətti reqressiya tənliyi ilə nə dərəcədə yaxşı xarakterizə etdiyimizi bilmək üçün marağımızı ödəyəcək.

Dəyişikliklərin qrafik kodu

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

Qrafik №3 "Dəyişikliklər, %"

Tam mükəmməl deyil, amma vəzifəmizi yerinə yetirdik.

Koeffisientlərin müəyyənləşdirilməsi üçün bir funksiya yazacağıq və bibilotekanı istifadə etməklə NumPy, dəqiq desək — iki funksiya yazacağıq: biri pseudoinvers matrisdən istifadə edən (praktikada tövsiyyə edilmir, çünki hesablama prosesi çətin və qeyri-sabitdir), digəri isə matris tənliyi istifadə edən.

Analitik həll kodu (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

Koeffisientlərin müəyyənləşdirilməsi üçün sərf olunan vaxtı müqayisə edəcəyik və , təqdim olunan 3 metodla.

Hesablama vaxtını hesablamaq üçün kod

print ' 33[1m' + ' 33[4m' + "Koeffisientlərin NumPy kitabxanası istifadə edilmədən hesablanma vaxtı:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Koeffisientlərin pseudoinvers matrisdən istifadə etməklə hesablanma vaxtı:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Koeffisientlərin matris tənliyindən istifadə etməklə hesablanma vaxtı:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Kiçik verilənlərdə, "özümüz yazdığımız" funksiya, Kramerin metodu ilə koeffisientləri tapır.

İndi koeffisientləri tapmanın digər üsullarına keçə bilərik və .

Qradient enməsi

Öncəliklə, gradientin nə olduğunu müəyyənləşdirək. Sadəcə desək, gradient - funksiyanın maksimum artım istiqamətini göstərən bir kəsikdir. Dağa qalxma ilə analoji, gradientin baxdığı yerdə, dağın zirvəsinə ən dik yuxarıdır. Dağı misal götürsək, əslində, ən sürətli şəkildə basamağa yüklənmək üçün ən dik aşağıdadır, yəni minimum - funksiyanın nə artmadığı, nə də azalmadığı yer. Bu yerdə törəmə sıfıra bərabər olacaq. Beləliklə, gradientə deyil, antibrebentə ehtiyacımız var. Antigandienti tapmaq üçün sadəcə gradienti -1 (minus bir) ilə vurmaq lazımdır.

Unutmayın ki, funksiyanın bir neçə minimumu ola bilər və aşağıdakı algoritm əsasında birinə düşsək, tapdığımızdan daha aşağı başqa bir minimumu tapa bilməyəcəyik. Rahatlaşaq, bu bizi narahat etmir! Bizim halımızda tək minimumla üzləşirik, çünki funksiyamız qrafikdə adi parabolu təmsil edir. Və hər birimiz məktəb matematikası kursundan yaxşı bilirik ki, parabolun yalnız bir minimumu olur.

İlk olaraq, gradientin nə üçün lazım olduğunu anladıqdan sonra, gradientin bir parça, yəni müəyyən koordinatları olan bir vektor olduğunu da başa düşdük. Bu koordinatlar, həmin əmsalların özüdür. və İndi gradient enişini həyata keçirə bilərik.

Başlamadan əvvəl, eniş alqoritmi haqqında bir neçə cümlə oxumağı təklif edirəm:

• Əmsalların koordinatlarını təxminən təsadüfi müəyyən edirik. və Məsələn, biz sıfıra yaxın əmsalları müəyyən edəcəyik. Bu, geniş yayılmış bir praktika olsa da, hər bir halda öz praktikası ola bilər.
• Koordinatdan birinci dərəcəli törəmə dəyərini çıxırıq. Beləliklə, əgər törəmə müsbət olarsa, funksiyanın artdığını qeyd edirik. Nəticədə, törəmə dəyərini çıxaraq artım istiqamətinin əksinə, yəni eniş istiqamətinə doğru irəliləyəcəyik. Əgər törəmə mənfi olarsa, deməli, funksiya həmin nöqtədə azalır və törəmə dəyərini çıxaraq eniş istiqamətinə doğru irəliləyəcəyik.
• Eyni əməliyyatı koordinatla icra edirik: həmin nöqtədəki törəmə dəyərini çıxırıq. .
• Minimuma keçməmək və uzaq kosmosa uçmamaq üçün eniş istiqamətində addımın ölçüsünü müəyyən etmək lazımdır. Ümumilikdə, addımı düzgün təyin etməklə necə dəyişdirmək barədə bir məqalə yazmaq olar ki, bu da hesablama xərclərini azaltsın. Lakin bu an əlimizdə bir qədər fərqli bir məsələ var və elmi üsullarla «sınaq» və ya sadə dillə desək, empirik metodla addımın ölçüsünü müəyyən edəcəyik.
• Verilmiş koordinatlardan və törəmə dəyərlərini çıxardıqdan sonra yeni koordinatlar əldə edirik. və Növbəti addımı (çıxma) artıq hesablanmış koordinatlardan həyata keçiririk. Və bu dövr davamlı olaraq yenidən baş verir, tələb olunan yaxınlaşmaya çatana qədər.

Hamısı! İndi Mariya çuxurunun ən dərin çuxurunu axtarmağa hazırıq. Başlayaq.

Gradient enişi üçün kod

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

Mariya çuxurunun dibinə enmişik və orada eyni əmsal dəyərlərini aşkar etdik və , bu da gözlənilən bir haldır.

Başqa bir eniş edəcəyik, yalnız bu dəfə dərinlik aparatımızın içini fərqli texnologiyalar — yəni NumPy.

Gradient enişi üçün kod (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

Əmsal dəyərləri və dəyişməzdir.

Gradient enişi zamanı xətanın necə dəyişdiyinə, yəni hər addımda kvadrat fərqlərin cəminə baxaq.

Kvadrat fərqlərin cəmini göstərmək üçün kod

print 'Qrafik№4 "Hər addımda kvadrat fərqlərin cəmi"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Addımlar (Təkrarlama)', size=16)
plt.ylabel('Kvadrat fərqlərin cəmi', size=16)
plt.show()

Qrafik №4 «Gradient enişi zamanı kvadrat fərqlərin cəmi»

Qrafikdə hər addımla səhvin azaldığını və bir neçə iterasiyadan sonra demək olar ki, horizontal bir xətt gördüyümüzü görürük.

Son olaraq kodun icra vaxtındakı fərqi qiymətləndirək:

Gradient enişinin icra vaxtını müəyyənləşdirmək üçün kod

print ' 33[1m' + ' 33[4m' + "NumPy kitabxanası istifadə etmədən gradient enişinin icra vaxtı:" + ' 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 kitabxanası istifadə edərək gradient enişinin icra vaxtı:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Bəlkə də biz bir şeyləri düzgün etmirik, amma yenə də kitabxana istifadə etmədən yazılmış sadə "öz yazdığımız" funksiyanın kitabxana istifadə edən funksiyadan daha sürətli işlədiyini görürük. NumPy amma hesablamaların icrasındakı vaxtı keçməyi bacarır. NumPy.

Amma biz yerimizdə qalmırıq, sadə xətti tənliklər problemini həll etməyin digər bir maraqlı üsulunu öyrənməyə doğru irəliləyirik. Qarşınızda!

Stoxastik qradient enməsi

Stokastik gradient enişinin iş prinsipini daha yaxşı başa düşmək üçün, əvvəlcə onun adətən gradient enişindən fərqlərini müəyyənləşdirmək lazımdır. Biz gradient enişi vəziyyətində, törəmə tənliklərdə və nümunədəki bütün xüsusiyyətlərin və real cavabların cəmlərini (yəni, bütün cəmləri və ). Stokastik gradient enişində, nümunədəki bütün qiymətləri istifadə etməyəcəyik, əvəzində təsadüfi olaraq "nümunə indeksi" seçib onun dəyərlərini istifadə edəcəyik.

Məsələn, əgər indeks 3 (üç) olaraq seçilərsə, biz və , daha sonra dəyərləri törəmə tənliklərinə daxil edib yeni koordinatları müəyyən edirik. Sonra, koordinatları müəyyən etdikdən sonra, yenidən təsadüfi olaraq nümunə indeksini müəyyənləşdirib, indiki dəyərləri qismən törəmə tənliklərinə daxil edirik, yeni koordinatları müəyyən edirik və və s. aşkarlanma ilə tamamlanana qədər. İlk baxışda, bu necə işləyə bilər deyə düşünə bilərsiniz, amma işləyir. Ancaq qeyd etmək lazımdır ki, hər addımda səhv azalır, amma meyl, əslində, var.

Stokastik gradient enişinin adi gradient enişinə nisbətən üstünlükləri nələrdir? Əgər bizim nümunə ölçümüz çox böyükdürsə və on minlərlə dəyərlə ölçülürsə, məsələn, təsadüfi bir min dəyəri emal etmək, bütün nümunəni emal etməkdən daha asandır. Belə bir halda stokastik gradient enişi işə düşür. Bizim vəziyyətimizdə isə, əlbəttə ki, böyük fərq görməyəcəyik.

Koda baxırıq.

Stokastik gradient enişi üçün kod

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

Koeffisiyentlərə diqqətlə baxırıq və özümüzü "Necə olur ki, belədir?" sualı ilə tuturuq. Bizdə fərqli koeffisiyent dəyərləri alındı. və Bəlkə stoxastik gradient enməsi tənliyin daha optimal parametrlərini tapdı? Təəssüf ki, yox. Kvadrat sapmaların cəminə baxmaq kifayətdir və yeni koeffisient dəyərləri ilə səhv daha çox olduğunu görmək olar. Ruhsuzluğa qapılmağa tələsməyin. Səhv dəyişiminin qrafikini çəkək.

Stoxastik gradient enməsi üçün kvadrat sapmaların cəmini qrafikləşdirmək üçün kod

print 'Qrafik №5 "Kvadrat sapmaların cəmi üzrə addım-addım"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Addımlar (İterasiya)', size=16)
plt.ylabel('Kvadrat sapmaların cəmi', size=16)
plt.show()

Qrafik №5 «Stoxastik gradient enməsi zamanı kvadrat sapmaların cəmi»

Qrafikə baxanda hər şey yerinə oturur və indi hər şeyi düzəldəcəyik.

Bəs, nə baş verdi? Belə oldu. Təsadüfi bir ay seçəndə, bizim alqoritm həmin seçilmiş ay üçün gəlir hesablamasındakı səhvi azaltmağa çalışır. Sonra başqa bir ay seçirik və hesablamanı təkrarlayırıq, amma artıq ikinci seçilmiş ay üçün səhvi azaldırıq. İndi xatırlayaq ki, ilk iki ay sadə xətti reqressiya tənliyindən əhəmiyyətli dərəcədə sapır. Bu, deməkdir ki, bu iki aylardan hər hansı biri seçildikdə, hər birinin səhvini azaldarkən, alqoritmimiz bütün nümunə üzrə səhvi ciddi şəkildə artırır. Bəs nə etməliyik? Cavab sadədir: enmə addımını azaltmalıyıq. Çünki enmə addımını azaltmaqla, səhv də artıq yuxarıya, aşağıya “sıçrayacaq”. Düzdür, səhv “sıçrayış” etməyə davam edəcək, amma bu, o qədər sürətlə olmayacaq :) Gəlin yoxlayaq.

Kiçik addımla SGD-ni başlatmaq üçün kod

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

Qrafik №6 «Stoxastik gradient enməsi zamanı kvadrat sapmaların cəmi (80 min addım)»

Koeffisient dəyərləri yaxşılaşdı, amma yenə də ideal deyil. Hipotetik olaraq bunu belə düzəltmək olar. Məsələn, son 1000 iterasiyada minimal səhv olan koeffisient dəyərlərini seçirik. Lakin bunun üçün koeffisientlərin öz dəyərlərini də qeyd etməliyik. Bunu etməyəcəyik, bunun əvəzinə qrafikə baxaq. O, hamar görünür və səhv bərabər şəkildə azalmağa bənzəyir. Əslində isə belə deyil. İlk 1000 iterasiyanı izləyək və onları sonlarla müqayisə edək.

SGD qrafiki (ilk 1000 addım) üçün kod

print 'Qrafik №7 "Kvadrat sapmaların cəmi üzrə addım-addım. İlk 1000 iterasiya"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Addımlar (İterasiya)', size=16)
plt.ylabel('Kvadrat sapmaların cəmi', size=16)
plt.show()

print 'Qrafik №7 "Kvadrat sapmaların cəmi üzrə addım-addım. Sonuncu 1000 iterasiya"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Addımlar (İterasiya)', size=16)
plt.ylabel('Kvadrat sapmaların cəmi', size=16)
plt.show()

Cədvəl №7 "SGD ilə İtiraf Kvadratlarının Cəmi (ilk 1000 addım)"

Cədvəl №8 "SGD ilə İtiraf Kvadratlarının Cəmi (son 1000 addım)"

En başda eniş etməkdə olduğu müddətdə xətanın qabarıqlığı nisbətən sabit və dik bir şəkildə azalmaqdadır. Son iterasiyalarda xətanın 1,475 ətrafında döndüyünü və bəzən bu optimal dəyərə bərabər olduğunu, amma sonra yenidən yüksəldiyini görürük... Bir daha təkrar edirəm, biz koeffisiyentlərin dəyərlərini yazıb, daha sonra xətanın minimum olduğu vəziyyətləri seçə bilərik. Lakin, daha ciddi bir problemin ilə qarşılaşdıq: optimal dəyərlərə yaxın olan 80 min addım atmalıyıq (kodla baxın). Və bu, stoxastik gradient enişi ilə klassik gradient enişi arasında zaman qənaəti ideyasına ziddir. Nəyi düzəldə və inkişaf etdirə bilərik? İlk iterasiyalarda aşağıya doğru sabit bir irəliləyiş olduğunu görmək çətin deyil, buna görə ilk iterasiyalarda böyük addım atmağı, irəlilədikcə addımı azaldmağı düşünmək daha məqsədəuyğundur. Bu məqalədə bunu etməyəcəyik — o, artıq uzanmışdır. İstəyənlər özləri bunu necə edəcəyini fikirləşə bilərlər, çox çətin deyil 🙂 və İndi, kitabxana istifadə edərək stoxastik gradient enişini yerinə yetirəcəyik

(və daha əvvəl aşkar etdiyimiz daşlara dözməyəcəyik) NumPy Stoxastik gradient enişi üçün kod (NumPy)

Alınan dəyərlər, kitabxana istifadə etməyən enişlə demək olar ki, eynidir

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

. Bunun məntiqli olduğunu düşünmək olar. NumPyStoxastik gradient enişinin nə qədər zaman apardığını öyrənək.

SGD üçün hesablama müddətini müəyyən edən kod (80 min addım)

print ' 33[1m' + ' 33[4m' + "NumPy kitabxanası istifadə etmədən stoxastik gradient enişinin icra müddəti:" + ' 33[0m' %timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000) print '***************************************' printprint ' 33[1m' + ' 33[4m' + "NumPy kitabxanası istifadə edərək stoxastik gradient enişinin icra müddəti:" + ' 33[0m' %timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

Meşə içində daha irəlilədikcə tuman daha tündləşir: yenə "öz-özünə yazılmış" formul ən yaxşı nəticəni göstərir. Bütün bunlar, kitabxananın daha incə istifadə yollarının olmalı olduğunu düşündürür

, hansı ki, həqiqətən də hesablamaların sürətini artırır. Bu məqalədə bunlardan artıq öyrənməyəcəyik. Boş vaxtda düşünmək üçün bir şey olacaq :) NumPyXülasə edək

Xülasə etməzdən əvvəl, sevgili oxucumuzun şübhəsiz ki, düşündüyü bir sualı cavablandırmaq istəyirəm. Niyə, əslində, bu cür "sıxıntılar" ilə enişləri edirik, niyə dağda yuxarı və aşağı getməliyik (əsasən aşağı), arzuladığımız aşağılığa çatmaq üçün, əlimizdə dərhal bizi lazımi yerə teleport edən güclü və sadə bir vasitə - analitik həll var?

Xülasə etmədən əvvəl, əziz oxucumuzun ağıllarında yaranmış ola biləcək sualı cavablandırmaq istərdik. Niyə belə "əziyyət" çəkərək dağları dırmaşır, yuxarı və aşağı (əsasən aşağı) gedirik, əgər əlimizdə bizə lazım olan yerdə dərhal teleportasiyanı həyata keçirən, analitik həll kimi güclü və sadə bir cihaz varsa?

Bu suala verilən cavab açıq-aşkardadır. İndi biz çox sadə bir nümunəni müzakirə etdik, burada doğru cavab bir əlamətə bağlıdır . Həyatda bu cür nadir rast gəlinir, buna görə də fərz edək ki, bizdə 2, 30, 50 və ya daha çox əlamət var. Bunun üstündə hər bir əlamət üçün minlər, hətta on minlərlə qiymət əlavə edək. Bu halda analitik həll çətinlik çəkə bilər və çökə bilər. Öz növbəsində, gradient enişi və onun varyasiyaları bizi məqsədimizə - funksiyanın minimumuna yavaş-yavaş yaxınlaşdıracaq. Və sürət məsələsində narahat olmayın - biz mütləq bir daha müzakirə edəcəyik ki, biz addımın uzunluğunu (yəni sürətini) necə tənzimləyirik.

İndi isə əsas xülasəyə keçək.

Birincisi, ümid edirəm ki, məqalədə təqdim olunan material, yeni başlayan "data scientist"lərə sadə (və təkcə deyil) xətti regresiya tənliklərini necə həll etməyin başa düşülməsində kömək edəcək.

İkincisi, biz tənliyin bir neçə həll üsulunu müzakirə etdik. İndi vəziyyətdən asılı olaraq, biz ən uyğun olanını seçə bilərik.

Üçüncüsü, addım uzunluğunun gradient enişinin əlavə tənzimlemələrinin gücünü gördük. Bu parametrə laqeyd yanaşmaq olmaz. Yuxarıda qeyd edildiyi kimi, hesablama xərclərini azaltmaq məqsədilə, addım uzunluğunu enişin gedişində dəyişdirmək lazımdır.

Dördüncüsü, bizim halımızda "öz yazdığımız" funksiyalar ən yaxşı vaxt nəticəsini göstərdi. Bu, kitabxananın imkanlarından professional istifadə edilməməsi ilə bağlı ola bilər. NumPyAmma necə olursa-olsun, nəticə belədir. Bir tərəfdən, bəzən qəbul edilmiş fikirləri şübhə altında qoymaq vacibdir, digər tərəfdən isə hər şeyi çətinləşdirmək lazım deyil - əksinə, bəzən daha sadə bir həll yolu daha effektiv ola bilər. Və məqsədimiz sadə xətti regresiya tənliyi üçün üç yanaşmanı müzakirə etmək olduğundan, "öz yazdığımız" funksiyalarımız bizə tamamilə kifayət etdi.

Ədəbiyyat (ya da buna bənzər bir şey)

1. Xətti regresiya

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

2. Ən kiçik kvadratlar metodu

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Törəmə

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Gradient enişi

habr.com/ru/post/471458

habr.com/ru/post/307312

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

6. NumPy kitabxanası

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

Mənbə: habr.com