Menyelesaikan persamaan regresi linear mudah

Artikel ini membincangkan beberapa cara untuk menentukan persamaan matematik bagi garis regresi mudah (berpasangan).

Semua kaedah untuk menyelesaikan persamaan yang dibincangkan di sini adalah berdasarkan kaedah kuasa dua terkecil. Mari kita nyatakan kaedah seperti berikut:

• Penyelesaian analitikal
• Keturunan Kecerunan
• Penurunan kecerunan stokastik

Bagi setiap kaedah menyelesaikan persamaan garis lurus, artikel itu menyediakan pelbagai fungsi, yang kebanyakannya dibahagikan kepada yang ditulis tanpa menggunakan perpustakaan numpy dan yang digunakan untuk pengiraan numpy. Adalah dipercayai bahawa penggunaan mahir numpy akan mengurangkan kos pengkomputeran.

Semua kod yang diberikan dalam artikel ditulis dalam bahasa ular sawa 2.7 menggunakan Notebook Jupyter. Kod sumber dan fail dengan data sampel disiarkan pada Github

Artikel ini lebih ditujukan kepada kedua-dua pemula dan mereka yang telah secara beransur-ansur mula menguasai kajian bahagian yang sangat luas dalam kecerdasan buatan - pembelajaran mesin.

Untuk menggambarkan bahan, kami menggunakan contoh yang sangat mudah.

Contoh syarat

Kami mempunyai lima nilai yang mencirikan pergantungan Y daripada X (Jadual No. 1):

Jadual No. 1 "Syarat contoh"

Kami akan menganggap bahawa nilai ialah bulan dalam setahun, dan — pendapatan bulan ini. Dalam erti kata lain, pendapatan bergantung pada bulan dalam setahun, dan - satu-satunya tanda di mana pendapatan bergantung.

Contohnya begitu-begitu, kedua-duanya dari sudut pergantungan bersyarat hasil pada bulan dalam setahun, dan dari sudut pandangan bilangan nilai - terdapat sangat sedikit daripada mereka. Walau bagaimanapun, penyederhanaan sedemikian akan memungkinkan, seperti yang mereka katakan, untuk menerangkan, tidak selalu dengan mudah, bahan yang diasimilasikan oleh pemula. Dan juga kesederhanaan nombor akan membolehkan mereka yang ingin menyelesaikan contoh di atas kertas tanpa kos buruh yang ketara.

Mari kita anggap bahawa pergantungan yang diberikan dalam contoh boleh dianggarkan dengan baik oleh persamaan matematik garis regresi mudah (berpasangan) dalam bentuk:

mana adalah bulan di mana pendapatan diterima, - hasil yang sepadan dengan bulan, и ialah pekali regresi bagi garis anggaran.

Perhatikan bahawa pekali sering dipanggil cerun atau kecerunan garis anggaran; mewakili jumlah yang apabila ia berubah .

Jelas sekali, tugas kita dalam contoh adalah untuk memilih pekali sedemikian dalam persamaan и , di mana sisihan nilai hasil kami dikira mengikut bulan daripada jawapan sebenar, i.e. nilai yang ditunjukkan dalam sampel akan menjadi minimum.

Kaedah kuasa dua terkecil

Mengikut kaedah kuasa dua terkecil, sisihan hendaklah dikira dengan menduakannya. Teknik ini membolehkan anda mengelakkan pembatalan bersama penyimpangan jika mereka mempunyai tanda yang bertentangan. Sebagai contoh, jika dalam satu kes, sisihan adalah +5 (tambah lima), dan yang lain -5 (tolak lima), maka jumlah sisihan akan membatalkan satu sama lain dan berjumlah 0 (sifar). Adalah mungkin untuk tidak mengkuadratkan sisihan, tetapi menggunakan sifat modulus dan kemudian semua sisihan akan menjadi positif dan akan terkumpul. Kami tidak akan membincangkan perkara ini secara terperinci, tetapi hanya menunjukkan bahawa untuk kemudahan pengiraan, adalah lazim untuk mengkuadratkan sisihan.

Inilah rupa formula yang dengannya kita akan menentukan jumlah sisihan kuasa dua terkecil (ralat):

mana ialah fungsi penghampiran jawapan benar (iaitu, hasil yang kami kira),

adalah jawapan yang benar (hasil yang disediakan dalam sampel),

ialah indeks sampel (bilangan bulan di mana sisihan ditentukan)

Mari bezakan fungsi, tentukan persamaan pembezaan separa, dan bersedia untuk beralih kepada penyelesaian analisis. Tetapi pertama-tama, mari kita bersiar-siar singkat tentang apa itu pembezaan dan ingat makna geometri terbitan itu.

Pembezaan

Pembezaan ialah operasi mencari terbitan bagi suatu fungsi.

Derivatif digunakan untuk apa? Derivatif fungsi mencirikan kadar perubahan fungsi dan memberitahu kita arahnya. Jika derivatif pada titik tertentu adalah positif, maka fungsi bertambah; jika tidak, fungsi berkurang. Dan semakin besar nilai terbitan mutlak, semakin tinggi kadar perubahan nilai fungsi, serta semakin curam kecerunan graf fungsi.

Sebagai contoh, di bawah keadaan sistem koordinat Cartesan, nilai terbitan pada titik M(0,0) adalah sama dengan +25 bermakna pada titik tertentu, apabila nilai dialihkan ke kanan oleh unit konvensional, nilai meningkat sebanyak 25 unit konvensional. Pada graf ia kelihatan seperti peningkatan yang agak curam dalam nilai dari titik tertentu.

Contoh yang lain. Nilai terbitan adalah sama -0,1 bermakna apabila disesarkan setiap satu unit konvensional, nilai berkurangan hanya 0,1 unit konvensional. Pada masa yang sama, pada graf fungsi, kita boleh melihat cerun menurun yang hampir tidak ketara. Melukis analogi dengan gunung, seolah-olah kita menuruni cerun yang lembut dari gunung, tidak seperti contoh sebelumnya, di mana kita terpaksa mendaki puncak yang sangat curam :)

Oleh itu, selepas membezakan fungsi secara kebetulan и , kami mentakrifkan persamaan pembezaan separa tertib pertama. Selepas menentukan persamaan, kami akan menerima sistem dua persamaan, dengan menyelesaikannya kami akan dapat memilih nilai pekali tersebut и , yang mana nilai derivatif yang sepadan pada titik tertentu berubah dengan jumlah yang sangat, sangat kecil, dan dalam kes penyelesaian analitik tidak berubah sama sekali. Dengan kata lain, fungsi ralat pada pekali yang ditemui akan mencapai minimum, kerana nilai derivatif separa pada titik ini akan sama dengan sifar.

Jadi, mengikut peraturan pembezaan, persamaan derivatif separa bagi tertib pertama berkenaan dengan pekali akan mengambil borang:

Persamaan terbitan separa tertib pertama berkenaan dengan akan mengambil borang:

Akibatnya, kami menerima sistem persamaan yang mempunyai penyelesaian analitik yang agak mudah:

mulakan{persamaan*}
mulakan{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
tamat{kes}
tamat{persamaan*}

Sebelum menyelesaikan persamaan, mari kita pramuat, semak sama ada pemuatan adalah betul, dan format data.

Memuatkan dan memformat data

Perlu diingatkan bahawa disebabkan oleh fakta bahawa untuk penyelesaian analitikal, dan seterusnya untuk kecerunan dan keturunan kecerunan stokastik, kami akan menggunakan kod dalam dua variasi: menggunakan perpustakaan numpy dan tanpa menggunakannya, maka kami memerlukan pemformatan data yang sesuai (lihat kod).

Pemuatan data dan kod pemprosesan

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

Visualisasi

Sekarang, selepas kami, pertama, memuatkan data, kedua, menyemak ketepatan pemuatan dan akhirnya memformat data, kami akan menjalankan visualisasi pertama. Kaedah yang sering digunakan untuk ini ialah petak berpasangan perpustakaan Laut Laut. Dalam contoh kami, disebabkan bilangan yang terhad, tidak ada gunanya menggunakan perpustakaan Laut Laut. Kami akan menggunakan perpustakaan biasa Matplotlib dan lihat sahaja pada scatterplot.

Kod Scatterplot

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

Carta No. 1 "Pergantungan hasil pada bulan dalam setahun"

Penyelesaian analitikal

Mari gunakan alat yang paling biasa dalam ular sawa dan selesaikan sistem persamaan:

mulakan{persamaan*}
mulakan{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
tamat{kes}
tamat{persamaan*}

Mengikut peraturan Cramer kita akan mencari penentu am, serta penentu oleh dan oleh , selepas itu, membahagikan penentu dengan kepada penentu am - cari pekali , begitu juga kita dapati pekali .

Kod penyelesaian analitikal

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

Inilah yang kami dapat:

Jadi, nilai pekali telah dijumpai, jumlah sisihan kuasa dua telah ditetapkan. Mari kita lukis garis lurus pada histogram serakan mengikut pekali yang ditemui.

Kod baris regresi

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

Carta No. 2 “Jawapan yang betul dan dikira”

Anda boleh melihat graf sisihan bagi setiap bulan. Dalam kes kami, kami tidak akan memperoleh sebarang nilai praktikal yang ketara daripadanya, tetapi kami akan memenuhi rasa ingin tahu kami tentang sejauh mana persamaan regresi linear mudah mencirikan pergantungan hasil pada bulan dalam setahun.

Kod carta sisihan

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

Carta No. 3 "Sisihan, %"

Tidak sempurna, tetapi kami menyelesaikan tugas kami.

Mari kita tulis fungsi yang, untuk menentukan pekali и menggunakan perpustakaan numpy, lebih tepat lagi, kami akan menulis dua fungsi: satu menggunakan matriks pseudoinverse (tidak disyorkan dalam amalan, kerana prosesnya kompleks dan tidak stabil dari segi pengiraan), satu lagi menggunakan persamaan matriks.

Kod Penyelesaian Analitikal (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

Mari kita bandingkan masa yang digunakan untuk menentukan pekali и , selaras dengan 3 kaedah yang dibentangkan.

Kod untuk mengira masa pengiraan

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)

Dengan jumlah data yang kecil, fungsi "tulisan sendiri" muncul di hadapan, yang mencari pekali menggunakan kaedah Cramer.

Kini anda boleh beralih kepada cara lain untuk mencari pekali и .

Keturunan Kecerunan

Pertama, mari kita tentukan apa itu kecerunan. Ringkasnya, kecerunan ialah segmen yang menunjukkan arah pertumbuhan maksimum sesuatu fungsi. Secara analogi dengan mendaki gunung, di mana muka kecerunan adalah tempat pendakian paling curam ke puncak gunung itu. Membangunkan contoh dengan gunung, kita ingat bahawa sebenarnya kita memerlukan keturunan yang paling curam untuk sampai ke tanah pamah secepat mungkin, iaitu minimum - tempat di mana fungsi tidak meningkat atau berkurangan. Pada ketika ini derivatif akan sama dengan sifar. Oleh itu, kita tidak memerlukan kecerunan, tetapi antikecerunan. Untuk mencari antikecerunan anda hanya perlu mendarabkan kecerunan dengan -1 (tolak satu).

Marilah kita memberi perhatian kepada fakta bahawa fungsi boleh mempunyai beberapa minima, dan setelah turun ke salah satu daripadanya menggunakan algoritma yang dicadangkan di bawah, kita tidak akan dapat mencari minimum lain, yang mungkin lebih rendah daripada yang ditemui. Mari berehat, ini bukan ancaman kepada kami! Dalam kes kami, kami berurusan dengan minimum tunggal, sejak fungsi kami pada graf ialah parabola sekata. Dan seperti yang kita semua harus tahu dengan baik dari kursus matematik sekolah kita, parabola hanya mempunyai satu minimum.

Selepas kami mengetahui mengapa kami memerlukan kecerunan, dan juga bahawa kecerunan adalah segmen, iaitu, vektor dengan koordinat yang diberikan, yang merupakan pekali yang sama. и kita boleh melaksanakan penurunan kecerunan.

Sebelum memulakan, saya cadangkan membaca hanya beberapa ayat tentang algoritma keturunan:

• Kami menentukan secara pseudo-rawak koordinat bagi pekali и . Dalam contoh kami, kami akan menentukan pekali berhampiran sifar. Ini adalah amalan biasa, tetapi setiap kes mungkin mempunyai amalannya sendiri.
• Dari koordinat tolak nilai terbitan separa tertib pertama pada titik itu . Jadi, jika derivatifnya positif, maka fungsinya bertambah. Oleh itu, dengan menolak nilai terbitan, kita akan bergerak ke arah pertumbuhan yang bertentangan, iaitu ke arah keturunan. Jika terbitan adalah negatif, maka fungsi pada titik ini berkurangan dan dengan menolak nilai terbitan kita bergerak ke arah penurunan.
• Kami menjalankan operasi yang sama dengan koordinat : tolak nilai terbitan separa pada titik .
• Agar tidak melompat ke atas minimum dan terbang ke angkasa lepas, perlu menetapkan saiz langkah ke arah penurunan. Secara umum, anda boleh menulis keseluruhan artikel tentang cara menetapkan langkah dengan betul dan cara mengubahnya semasa proses penurunan untuk mengurangkan kos pengiraan. Tetapi sekarang kami mempunyai tugas yang sedikit berbeza di hadapan kami, dan kami akan menetapkan saiz langkah menggunakan kaedah saintifik "poke" atau, seperti yang mereka katakan dalam bahasa biasa, secara empirik.
• Sebaik sahaja kita dari koordinat yang diberikan и tolak nilai derivatif, kita mendapat koordinat baharu и . Kami mengambil langkah seterusnya (penolakan), sudah dari koordinat yang dikira. Jadi kitaran bermula lagi dan lagi, sehingga penumpuan yang diperlukan dicapai.

Semua! Kini kami bersedia untuk pergi mencari gaung paling dalam di Palung Mariana. Mari kita mulakan.

Kod untuk keturunan kecerunan

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

Kami menyelam ke bahagian paling bawah Palung Mariana dan di sana kami menemui semua nilai pekali yang sama и , iaitu apa yang dijangkakan.

Mari kita menyelam lagi, cuma kali ini kenderaan laut dalam kita akan dipenuhi dengan teknologi lain iaitu perpustakaan numpy.

Kod untuk keturunan kecerunan (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

Nilai pekali и tidak boleh diubah.

Mari kita lihat bagaimana ralat berubah semasa penurunan kecerunan, iaitu, bagaimana jumlah sisihan kuasa dua berubah dengan setiap langkah.

Kod untuk memplot jumlah sisihan kuasa dua

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

Graf No. 4 "Jumlah sisihan kuasa dua semasa penurunan kecerunan"

Pada graf kita melihat bahawa dengan setiap langkah ralat berkurangan, dan selepas beberapa lelaran tertentu kita memerhatikan garis yang hampir mendatar.

Akhir sekali, mari kita anggarkan perbezaan dalam masa pelaksanaan kod:

Kod untuk menentukan masa pengiraan penurunan kecerunan

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)

Mungkin kita melakukan sesuatu yang salah, tetapi sekali lagi ia adalah fungsi "tulisan rumah" yang mudah yang tidak menggunakan perpustakaan numpy mengatasi masa pengiraan fungsi menggunakan perpustakaan numpy.

Tetapi kita tidak berdiam diri, tetapi sedang bergerak ke arah mengkaji satu lagi cara yang menarik untuk menyelesaikan persamaan regresi linear mudah. jumpa!

Penurunan kecerunan stokastik

Untuk memahami dengan cepat prinsip operasi keturunan kecerunan stokastik, adalah lebih baik untuk menentukan perbezaannya daripada keturunan kecerunan biasa. Kami, dalam kes keturunan kecerunan, dalam persamaan derivatif bagi и menggunakan jumlah nilai semua ciri dan jawapan benar yang tersedia dalam sampel (iaitu, jumlah semua и ). Dalam keturunan kecerunan stokastik, kami tidak akan menggunakan semua nilai yang ada dalam sampel, tetapi sebaliknya, pilih secara rawak pseudo indeks sampel dan gunakan nilainya.

Sebagai contoh, jika indeks ditentukan menjadi nombor 3 (tiga), maka kita mengambil nilainya и , kemudian kita menggantikan nilai ke dalam persamaan derivatif dan menentukan koordinat baharu. Kemudian, setelah menentukan koordinat, kami sekali lagi pseudo-rawak menentukan indeks sampel, menggantikan nilai yang sepadan dengan indeks ke dalam persamaan pembezaan separa, dan menentukan koordinat dengan cara yang baru и dan lain-lain. sehingga penumpuan bertukar hijau. Pada pandangan pertama, ia mungkin tidak kelihatan seperti ini boleh berfungsi sama sekali, tetapi ia berfungsi. Memang benar bahawa perlu diperhatikan bahawa ralat tidak berkurangan dengan setiap langkah, tetapi pasti ada kecenderungan.

Apakah kelebihan penurunan kecerunan stokastik berbanding konvensional? Jika saiz sampel kami sangat besar dan diukur dalam puluhan ribu nilai, maka lebih mudah untuk memproses, katakan, seribu rawak daripadanya, berbanding keseluruhan sampel. Di sinilah penurunan kecerunan stokastik dimainkan. Dalam kes kami, sudah tentu, kami tidak akan melihat banyak perbezaan.

Mari lihat kodnya.

Kod untuk keturunan kecerunan stokastik

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

Kami melihat dengan teliti pada pekali dan mendapati diri kami bertanya soalan "Bagaimana ini boleh berlaku?" Kami mendapat nilai pekali lain и . Mungkin keturunan kecerunan stokastik telah menemui parameter yang lebih optimum untuk persamaan? Malangnya tidak. Ia cukup untuk melihat jumlah sisihan kuasa dua dan melihat bahawa dengan nilai baru pekali, ralatnya lebih besar. Kami tidak tergesa-gesa untuk berputus asa. Mari bina graf perubahan ralat.

Kod untuk memplot jumlah sisihan kuasa dua dalam keturunan kecerunan stokastik

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

Graf No. 5 "Jumlah sisihan kuasa dua semasa penurunan kecerunan stokastik"

Melihat kepada jadual, segala-galanya jatuh ke tempatnya dan kini kami akan membetulkan segala-galanya.

Jadi apa yang berlaku? Perkara berikut berlaku. Apabila kami memilih bulan secara rawak, maka untuk bulan yang dipilih itulah algoritma kami berusaha untuk mengurangkan ralat dalam mengira hasil. Kemudian kami memilih bulan lain dan ulangi pengiraan, tetapi kami mengurangkan ralat untuk bulan kedua yang dipilih. Sekarang ingat bahawa dua bulan pertama menyimpang dengan ketara daripada garis persamaan regresi linear mudah. Ini bermakna apabila mana-mana dua bulan ini dipilih, dengan mengurangkan ralat setiap satu daripadanya, algoritma kami secara serius meningkatkan ralat untuk keseluruhan sampel. Jadi apa yang perlu dilakukan? Jawapannya mudah: anda perlu mengurangkan langkah turun. Lagipun, dengan mengurangkan langkah penurunan, ralat juga akan berhenti "melompat" ke atas dan ke bawah. Atau sebaliknya, ralat "melompat" tidak akan berhenti, tetapi ia tidak akan melakukannya dengan cepat :) Mari semak.

Kod untuk menjalankan SGD dengan kenaikan yang lebih kecil

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

Graf No. 6 "Jumlah sisihan kuasa dua semasa penurunan kecerunan stokastik (80 ribu langkah)"

Pekali telah bertambah baik, tetapi masih tidak ideal. Secara hipotesis, ini boleh diperbetulkan dengan cara ini. Kami memilih, sebagai contoh, dalam 1000 lelaran terakhir nilai pekali yang ralat minimum dibuat. Benar, untuk ini kita juga perlu menulis nilai pekali itu sendiri. Kami tidak akan melakukan ini, sebaliknya memberi perhatian kepada jadual. Ia kelihatan lancar dan ralat kelihatan berkurangan sama rata. Sebenarnya ini tidak benar. Mari lihat 1000 lelaran pertama dan bandingkan dengan yang terakhir.

Kod untuk carta SGD (1000 langkah pertama)

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

Graf No. 7 “Jumlah sisihan kuasa dua SGD (1000 langkah pertama)”

Graf No. 8 “Jumlah sisihan kuasa dua SGD (1000 langkah terakhir)”

Pada awal-awal keturunan, kami memerhatikan penurunan ralat yang agak seragam dan curam. Dalam lelaran terakhir, kita melihat bahawa ralat berlaku di sekitar nilai 1,475 dan pada beberapa ketika menyamai nilai optimum ini, tetapi kemudian ia masih meningkat... Saya ulangi, anda boleh menulis nilai-nilai pekali и , dan kemudian pilih yang ralatnya adalah minimum. Walau bagaimanapun, kami mempunyai masalah yang lebih serius: kami terpaksa mengambil 80 ribu langkah (lihat kod) untuk mendapatkan nilai yang hampir optimum. Dan ini sudah bercanggah dengan idea untuk menjimatkan masa pengiraan dengan keturunan kecerunan stokastik berbanding dengan keturunan kecerunan. Apa yang boleh diperbetulkan dan diperbaiki? Tidak sukar untuk menyedari bahawa dalam lelaran pertama kita dengan yakin akan menurun dan, oleh itu, kita harus meninggalkan satu langkah besar dalam lelaran pertama dan mengurangkan langkah semasa kita bergerak ke hadapan. Kami tidak akan melakukan ini dalam artikel ini - ia sudah terlalu panjang. Mereka yang ingin boleh berfikir sendiri bagaimana untuk melakukan ini, ia tidak sukar :)

Sekarang mari kita lakukan penurunan kecerunan stokastik menggunakan perpustakaan numpy (dan jangan kita tersandung pada batu yang kita kenal pasti sebelum ini)

Kod untuk Keturunan Kecerunan Stokastik (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

Nilainya ternyata hampir sama seperti ketika menurun tanpa menggunakan numpy. Walau bagaimanapun, ini adalah logik.

Mari kita ketahui berapa lama penurunan kecerunan stokastik mengambil masa kita.

Kod untuk menentukan masa pengiraan SGD (80 ribu langkah)

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)

Semakin jauh ke dalam hutan, semakin gelap awan: sekali lagi, formula "tulisan sendiri" menunjukkan hasil terbaik. Semua ini menunjukkan bahawa mesti ada cara yang lebih halus untuk menggunakan perpustakaan numpy, yang benar-benar mempercepatkan operasi pengiraan. Dalam artikel ini kita tidak akan belajar tentang mereka. Akan ada sesuatu untuk difikirkan pada masa lapang anda :)

Meringkaskan

Sebelum meringkaskan, saya ingin menjawab soalan yang kemungkinan besar timbul daripada pembaca yang dikasihi. Mengapa, sebenarnya, "penyiksaan" seperti itu dengan keturunan, mengapa kita perlu berjalan naik dan turun gunung (kebanyakannya turun) untuk mencari tanah pamah yang berharga, jika kita ada di tangan kita alat yang kuat dan mudah, di bentuk penyelesaian analitikal, yang membawa kita serta-merta ke tempat Betul?

Jawapan kepada soalan ini terletak di permukaan. Sekarang kita telah melihat contoh yang sangat mudah, di mana jawapan yang benar adalah bergantung pada satu tanda . Anda tidak sering melihat ini dalam hidup, jadi mari bayangkan bahawa kita mempunyai 2, 30, 50 atau lebih tanda. Mari tambahkan pada ini beribu-ribu, malah berpuluh-puluh ribu nilai untuk setiap atribut. Dalam kes ini, penyelesaian analitik mungkin tidak dapat menahan ujian dan gagal. Sebaliknya, keturunan kecerunan dan variasinya akan secara perlahan tetapi pasti membawa kita lebih dekat kepada matlamat - fungsi minimum. Dan jangan risau tentang kelajuan - kami mungkin akan melihat cara yang membolehkan kami menetapkan dan mengawal panjang langkah (iaitu kelajuan).

Dan sekarang ringkasan ringkas sebenar.

Pertama, saya berharap bahan yang dibentangkan dalam artikel itu akan membantu memulakan "saintis data" dalam memahami cara menyelesaikan persamaan regresi linear yang mudah (dan bukan sahaja).

Kedua, kami melihat beberapa cara untuk menyelesaikan persamaan. Sekarang, bergantung kepada keadaan, kita boleh memilih yang paling sesuai untuk menyelesaikan masalah.

Ketiga, kami melihat kuasa tetapan tambahan, iaitu panjang langkah penurunan kecerunan. Parameter ini tidak boleh diabaikan. Seperti yang dinyatakan di atas, untuk mengurangkan kos pengiraan, panjang langkah perlu diubah semasa penurunan.

Keempat, dalam kes kami, fungsi "tulisan rumah" menunjukkan hasil masa terbaik untuk pengiraan. Ini mungkin disebabkan oleh penggunaan keupayaan perpustakaan yang paling profesional numpy. Tetapi walau bagaimanapun, kesimpulan berikut menunjukkan dirinya sendiri. Di satu pihak, kadangkala patut mempersoalkan pendapat yang telah ditetapkan, dan sebaliknya, ia tidak selalu bernilai merumitkan segala-galanya - sebaliknya, kadang-kadang cara yang lebih mudah untuk menyelesaikan masalah adalah lebih berkesan. Dan kerana matlamat kami adalah untuk menganalisis tiga pendekatan untuk menyelesaikan persamaan regresi linear yang mudah, penggunaan fungsi "tulisan sendiri" adalah cukup untuk kami.

Sastera (atau sesuatu seperti itu)

1. Regresi linear

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

2. Kaedah kuasa dua terkecil

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Terbitan

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Kecerunan

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Keturunan kecerunan

habr.com/ms/post/471458

habr.com/ms/post/307312

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

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

Sumber: www.habr.com