Memecahkan persamaan regresi linier sederhana

Artikel ini membahas beberapa cara menentukan persamaan matematis garis regresi sederhana (berpasangan).

Semua metode penyelesaian persamaan yang dibahas di sini didasarkan pada metode kuadrat terkecil. Mari kita nyatakan metodenya sebagai berikut:

• Solusi analitis
• Penurunan Gradien
• Penurunan gradien stokastik

Untuk masing-masing metode penyelesaian persamaan garis lurus, artikel ini menyediakan berbagai fungsi, yang sebagian besar dibagi menjadi yang ditulis tanpa menggunakan perpustakaan JumlahPy dan yang digunakan untuk perhitungan JumlahPy. Hal ini diyakini bahwa penggunaan terampil JumlahPy akan mengurangi biaya komputasi.

Semua kode yang diberikan dalam artikel ditulis dalam bahasa tersebut ular piton 2.7 menggunakan Notebook Jupyter. Kode sumber dan file dengan data sampel diposting Github

Artikel ini lebih ditujukan untuk pemula dan mereka yang secara bertahap mulai menguasai studi bagian yang sangat luas dalam kecerdasan buatan - pembelajaran mesin.

Untuk mengilustrasikan materi, kami menggunakan contoh yang sangat sederhana.

Contoh kondisi

Kami memiliki lima nilai yang menjadi ciri ketergantungan Y dari X (Tabel No. 1):

Tabel No. 1 “Contoh kondisi”

Kami akan berasumsi bahwa nilainya adalah bulan dalam setahun, dan — pendapatan bulan ini. Dengan kata lain, pendapatan bergantung pada bulan dalam setahun, dan - satu-satunya tanda yang menjadi sandaran pendapatan.

Contohnya biasa saja, baik dari sudut pandang ketergantungan bersyarat pendapatan pada bulan dalam setahun, dan dari sudut pandang jumlah nilai - jumlahnya sangat sedikit. Namun, penyederhanaan seperti itu akan memungkinkan, seperti yang mereka katakan, untuk menjelaskan, tidak selalu dengan mudah, materi yang diasimilasi oleh pemula. Dan juga kesederhanaan angka-angkanya akan memungkinkan mereka yang ingin menyelesaikan contoh di atas kertas tanpa biaya tenaga kerja yang signifikan.

Mari kita asumsikan bahwa ketergantungan yang diberikan dalam contoh dapat didekati dengan baik menggunakan persamaan matematis dari garis regresi sederhana (berpasangan) berbentuk:

dimana adalah bulan di mana pendapatan diterima, — pendapatan sesuai dengan bulan tersebut, и adalah koefisien regresi dari garis estimasi.

Perhatikan bahwa koefisien sering disebut kemiringan atau gradien garis perkiraan; mewakili jumlah dimana ketika itu berubah .

Jelasnya, tugas kita dalam contoh ini adalah memilih koefisien-koefisien tersebut dalam persamaan и , yang merupakan penyimpangan nilai pendapatan yang kami hitung berdasarkan bulan dari jawaban sebenarnya, yaitu. nilai yang disajikan dalam sampel akan minimal.

Metode kuadrat terkecil

Menurut metode kuadrat terkecil, simpangan harus dihitung dengan mengkuadratkannya. Teknik ini memungkinkan Anda untuk menghindari saling pembatalan penyimpangan jika memiliki tanda yang berlawanan. Misalnya jika dalam satu kasus terjadi penyimpangan +5 (ditambah lima), dan yang lainnya -5 (dikurangi lima), maka jumlah simpangan tersebut akan saling menghilangkan dan berjumlah 0 (nol). Anda tidak perlu mengkuadratkan deviasinya, tetapi gunakan properti modulus dan semua deviasi akan menjadi positif dan terakumulasi. Kami tidak akan membahas hal ini secara rinci, tetapi hanya menunjukkan bahwa untuk kenyamanan perhitungan, biasanya deviasi dikuadratkan.

Seperti inilah rumus yang akan digunakan untuk menentukan jumlah deviasi kuadrat (kesalahan) terkecil:

dimana adalah fungsi perkiraan jawaban yang benar (yaitu, pendapatan yang kami hitung),

adalah jawaban yang benar (pendapatan yang diberikan dalam sampel),

adalah indeks sampel (jumlah bulan di mana deviasi ditentukan)

Mari kita bedakan fungsinya, tentukan persamaan diferensial parsialnya, dan bersiap untuk beralih ke solusi analitis. Namun pertama-tama, mari kita lihat sekilas tentang apa itu diferensiasi dan mengingat arti geometri dari turunan.

Diferensiasi

Diferensiasi adalah operasi mencari turunan suatu fungsi.

Untuk apa turunannya digunakan? Turunan suatu fungsi mencirikan laju perubahan fungsi tersebut dan memberi tahu kita arahnya. Jika turunan pada suatu titik bernilai positif, maka fungsi tersebut meningkat; jika tidak, fungsinya menurun. Dan semakin besar nilai turunan absolutnya maka semakin tinggi pula laju perubahan nilai fungsinya, serta semakin curam kemiringan grafik fungsinya.

Misalnya, pada kondisi sistem koordinat kartesius, nilai turunan di titik M(0,0) sama dengan +25 artinya pada suatu titik tertentu, ketika nilainya digeser ke kanan dengan satuan konvensional, nilai meningkat sebesar 25 unit konvensional. Pada grafik terlihat kenaikan nilai yang cukup tajam dari suatu titik tertentu.

Contoh lain. Nilai turunannya sama -0,1 berarti ketika dipindahkan per satu unit konvensional, nilai berkurang hanya 0,1 unit konvensional. Pada saat yang sama, pada grafik fungsi, kita dapat mengamati kemiringan ke bawah yang hampir tidak terlihat. Jika dianalogikan dengan gunung, kita seolah-olah sedang menuruni lereng gunung yang landai dengan sangat perlahan, berbeda dengan contoh sebelumnya yang harus mendaki puncak yang sangat curam :)

Jadi, setelah membedakan fungsinya secara kebetulan и , kita mendefinisikan persamaan diferensial parsial orde 1. Setelah menentukan persamaan, kita akan menerima sistem dua persamaan, dengan menyelesaikannya kita akan dapat memilih nilai koefisien tersebut и , yang nilai turunannya pada titik tertentu berubah dengan jumlah yang sangat, sangat kecil, dan dalam kasus solusi analitis tidak berubah sama sekali. Dengan kata lain, fungsi kesalahan pada koefisien yang ditemukan akan mencapai minimum, karena nilai turunan parsial pada titik-titik tersebut akan sama dengan nol.

Jadi, menurut aturan diferensiasi, persamaan turunan parsial orde 1 terhadap koefisien akan berbentuk:

Persamaan turunan parsial orde 1 terhadap akan berbentuk:

Hasilnya, kami memperoleh sistem persamaan yang memiliki solusi analitis yang cukup sederhana:

mulai{persamaan*}
mulai{kasus}
na + bjumlahbatas_{i=1}^nx_i — jumlahbatas_{i=1}^ny_i = 0

batas jumlah_{i=1}^nx_i(a +bjumlahbatas_{i=1}^nx_i — batas jumlah_{i=1}^ny_i) = 0
akhir{kasus}
akhir{persamaan*}

Sebelum menyelesaikan persamaan, mari kita lakukan pramuat, periksa apakah pemuatan sudah benar, dan format datanya.

Memuat dan memformat data

Perlu dicatat bahwa karena fakta bahwa untuk solusi analitis, dan selanjutnya untuk penurunan gradien gradien dan stokastik, kita akan menggunakan kode dalam dua variasi: menggunakan perpustakaan JumlahPy dan tanpa menggunakannya, maka kita memerlukan format data yang sesuai (lihat kode).

Pemuatan data dan kode pemrosesan

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

Nah setelah kita pertama memuat data, kedua memeriksa kebenaran pemuatan dan terakhir memformat data, kita akan melakukan visualisasi pertama. Cara yang sering digunakan untuk ini adalah plot pasangan perpustakaan Yg keturunan dr laut. Dalam contoh kita, karena terbatasnya jumlah, tidak ada gunanya menggunakan perpustakaan Yg keturunan dr laut. Kami akan menggunakan perpustakaan biasa Matplotlib dan lihat saja diagram sebarnya.

Kode plot sebar

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

Bagan No. 1 “Ketergantungan pendapatan pada bulan dalam setahun”

Solusi analitis

Mari gunakan alat yang paling umum ular sanca dan selesaikan sistem persamaan:

mulai{persamaan*}
mulai{kasus}
na + bjumlahbatas_{i=1}^nx_i — jumlahbatas_{i=1}^ny_i = 0

batas jumlah_{i=1}^nx_i(a +bjumlahbatas_{i=1}^nx_i — batas jumlah_{i=1}^ny_i) = 0
akhir{kasus}
akhir{persamaan*}

Menurut aturan Cramer kita akan menemukan determinan umum, serta determinan dengan dan , setelah itu membagi determinannya dengan ke determinan umum - temukan koefisiennya , demikian pula kita menemukan koefisiennya .

Kode solusi analitis

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

Jadi, nilai koefisien telah ditemukan, jumlah simpangan kuadrat telah ditentukan. Mari kita menggambar garis lurus pada histogram hamburan sesuai dengan koefisien yang ditemukan.

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

Bagan No. 2 “Jawaban yang benar dan penuh perhitungan”

Anda dapat melihat grafik deviasi setiap bulannya. Dalam kasus kami, kami tidak akan memperoleh nilai praktis apa pun darinya, namun kami akan memuaskan rasa ingin tahu kami tentang seberapa baik persamaan regresi linier sederhana mencirikan ketergantungan pendapatan pada bulan dalam setahun.

Kode grafik deviasi

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

Bagan No. 3 “Penyimpangan, %”

Tidak sempurna, tapi kami menyelesaikan tugas kami.

Mari kita tuliskan suatu fungsi untuk menentukan koefisien и menggunakan perpustakaan JumlahPy, lebih tepatnya, kita akan menulis dua fungsi: satu menggunakan matriks pseudoinvers (tidak disarankan dalam praktiknya, karena prosesnya rumit secara komputasi dan tidak stabil), yang lain menggunakan persamaan matriks.

Kode Solusi Analitik (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 waktu yang dihabiskan untuk menentukan koefisien и , sesuai dengan 3 metode yang disajikan.

Kode untuk menghitung waktu perhitungan

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 sedikit, fungsi “yang ditulis sendiri” akan muncul, yang mencari koefisien menggunakan metode Cramer.

Sekarang Anda dapat beralih ke cara lain untuk mencari koefisien и .

Penurunan Gradien

Pertama, mari kita definisikan apa itu gradien. Sederhananya, gradien adalah segmen yang menunjukkan arah pertumbuhan maksimum suatu fungsi. Analoginya dengan mendaki gunung, dimana kemiringan mukanya merupakan tempat pendakian paling curam menuju puncak gunung. Mengembangkan contoh gunung, kita ingat bahwa sebenarnya kita membutuhkan turunan yang paling curam agar bisa mencapai dataran rendah secepat mungkin, yaitu minimum - tempat yang fungsinya tidak bertambah atau berkurang. Pada titik ini turunannya akan sama dengan nol. Oleh karena itu, kita tidak memerlukan gradien, melainkan antigradien. Untuk mencari antigradien, Anda hanya perlu mengalikan gradien dengan -1 (kurang satu).

Mari kita perhatikan fakta bahwa suatu fungsi dapat memiliki beberapa nilai minimum, dan dengan turun ke salah satunya menggunakan algoritma yang diusulkan di bawah ini, kita tidak akan dapat menemukan nilai minimum lain, yang mungkin lebih rendah dari nilai minimum yang ditemukan. Tenang saja, ini bukan ancaman bagi kami! Dalam kasus kita, kita berhadapan dengan satu minimum, karena fungsi kita pada grafik adalah parabola beraturan. Dan seperti yang kita semua ketahui dengan baik dari pelajaran matematika sekolah, parabola hanya memiliki satu minimum.

Setelah kita mengetahui mengapa kita memerlukan gradien, dan juga bahwa gradien adalah sebuah segmen, yaitu vektor dengan koordinat tertentu, yang koefisiennya persis sama и kita dapat menerapkan penurunan gradien.

Sebelum memulai, saya sarankan membaca beberapa kalimat tentang algoritma keturunan:

• Kami menentukan koordinat koefisien secara pseudo-acak и . Dalam contoh kita, kita akan mendefinisikan koefisien mendekati nol. Ini adalah praktik umum, namun setiap kasus mungkin memiliki praktiknya sendiri.
• Dari koordinat kurangi nilai turunan parsial orde 1 di titik tersebut . Jadi, jika turunannya positif, maka fungsinya bertambah. Oleh karena itu, dengan mengurangkan nilai turunannya, kita akan bergerak berlawanan arah dengan pertumbuhannya, yaitu ke arah turun. Jika turunannya negatif, maka fungsi pada titik ini berkurang dan dengan mengurangkan nilai turunannya kita bergerak ke arah turun.
• Kami melakukan operasi serupa dengan koordinat : mengurangkan nilai turunan parsial di titik tersebut .
• Agar tidak melompati batas minimum dan terbang ke luar angkasa, perlu mengatur ukuran langkah ke arah turun. Secara umum, Anda dapat menulis artikel lengkap tentang cara menyetel langkah dengan benar dan cara mengubahnya selama proses penurunan untuk mengurangi biaya komputasi. Tapi sekarang kita memiliki tugas yang sedikit berbeda di depan kita, dan kita akan menentukan ukuran langkah menggunakan metode ilmiah “poke” atau, seperti yang mereka katakan dalam bahasa umum, secara empiris.
• Setelah kita berada dari koordinat yang diberikan и kurangi nilai turunannya, kita mendapatkan koordinat baru и . Kita ambil langkah selanjutnya (pengurangan), sudah dari koordinat yang dihitung. Dan siklus tersebut dimulai lagi dan lagi, hingga konvergensi yang diperlukan tercapai.

Semua! Sekarang kami siap berangkat mencari ngarai terdalam di Palung Mariana. Mari kita mulai.

Kode untuk penurunan gradien

# напишем функцию градиентного спуска без использования библиотеки 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 dasar Palung Mariana dan di sana kami menemukan semua nilai koefisien yang sama и , dan itulah yang diharapkan.

Mari kita menyelam lagi, hanya saja kali ini wahana laut dalam kita akan diisi dengan teknologi lain yaitu perpustakaan JumlahPy.

Kode untuk penurunan gradien (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 koefisien и tidak dapat diubah.

Mari kita lihat bagaimana kesalahan berubah selama penurunan gradien, yaitu bagaimana jumlah deviasi kuadrat berubah pada setiap langkah.

Kode untuk memplot jumlah deviasi kuadrat

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

Grafik No. 4 “Jumlah simpangan kuadrat selama penurunan gradien”

Pada grafik kita melihat bahwa dengan setiap langkah kesalahannya berkurang, dan setelah beberapa iterasi tertentu kita mengamati garis yang hampir horizontal.

Terakhir, mari kita perkirakan perbedaan waktu eksekusi kode:

Kode untuk menentukan waktu perhitungan penurunan gradien

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, tapi sekali lagi ini adalah fungsi sederhana yang “ditulis di rumah” yang tidak menggunakan perpustakaan JumlahPy mengungguli waktu perhitungan suatu fungsi menggunakan perpustakaan JumlahPy.

Namun kami tidak tinggal diam, tetapi kami sedang mempelajari cara menarik lainnya untuk menyelesaikan persamaan regresi linier sederhana. Bertemu!

Penurunan gradien stokastik

Untuk memahami dengan cepat prinsip pengoperasian penurunan gradien stokastik, lebih baik menentukan perbedaannya dari penurunan gradien biasa. Kita, dalam kasus penurunan gradien, dalam persamaan turunan dari и menggunakan jumlah nilai semua fitur dan jawaban benar yang tersedia dalam sampel (yaitu, jumlah semuanya и ). Dalam penurunan gradien stokastik, kita tidak akan menggunakan semua nilai yang ada dalam sampel, melainkan memilih secara acak semu apa yang disebut indeks sampel dan menggunakan nilainya.

Misalnya indeksnya bernomor 3 (tiga), maka kita ambil nilainya и , lalu kita substitusikan nilainya ke dalam persamaan turunan dan tentukan koordinat baru. Kemudian, setelah menentukan koordinat, kita kembali menentukan indeks sampel secara acak semu, mensubstitusikan nilai-nilai yang sesuai dengan indeks tersebut ke dalam persamaan diferensial parsial, dan menentukan koordinat dengan cara baru. и dll. sampai konvergensi berubah menjadi hijau. Pada pandangan pertama, ini mungkin tampak tidak berhasil sama sekali, tetapi ternyata berhasil. Memang benar bahwa perlu dicatat bahwa kesalahan tidak berkurang pada setiap langkah, tetapi pasti ada kecenderungannya.

Apa kelebihan penurunan gradien stokastik dibandingkan konvensional? Jika ukuran sampel kita sangat besar dan diukur dalam puluhan ribu nilai, maka akan jauh lebih mudah untuk memproses, katakanlah, ribuan nilai secara acak, dibandingkan keseluruhan sampel. Di sinilah penurunan gradien stokastik berperan. Dalam kasus kami, tentu saja, kami tidak akan melihat banyak perbedaan.

Mari kita lihat kodenya.

Kode untuk penurunan gradien 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])

Kita memperhatikan koefisiennya dengan cermat dan mendapati diri kita bertanya pada diri sendiri, “Bagaimana hal ini bisa terjadi?” Kami mendapat nilai koefisien lainnya и . Mungkin penurunan gradien stokastik telah menemukan parameter yang lebih optimal untuk persamaan tersebut? Sayangnya tidak ada. Cukup dengan melihat jumlah deviasi kuadrat dan melihat bahwa dengan nilai koefisien yang baru, kesalahannya lebih besar. Kami tidak terburu-buru untuk putus asa. Mari kita buat grafik perubahan kesalahan.

Kode untuk memplot jumlah deviasi kuadrat dalam penurunan gradien 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()

Grafik No. 5 “Jumlah deviasi kuadrat selama penurunan gradien stokastik”

Melihat jadwalnya, semuanya berjalan sesuai rencana dan sekarang kami akan memperbaiki semuanya.

Jadi apa yang terjadi? Hal berikut ini terjadi. Saat kami memilih suatu bulan secara acak, maka pada bulan yang dipilih itulah algoritme kami berupaya mengurangi kesalahan dalam menghitung pendapatan. Kemudian kita pilih bulan lain dan ulangi perhitungannya, tapi kita kurangi errornya untuk bulan kedua yang dipilih. Sekarang mari kita ingat bahwa dua bulan pertama menyimpang secara signifikan dari garis persamaan regresi linier sederhana. Ini berarti bahwa ketika salah satu dari dua bulan ini dipilih, dengan mengurangi kesalahan masing-masing bulan, algoritme kami secara serius meningkatkan kesalahan untuk keseluruhan sampel. Jadi apa yang harus dilakukan? Jawabannya sederhana: Anda perlu mengurangi langkah turun. Toh, dengan mengurangi langkah turun, error tersebut juga akan berhenti “melompat” ke atas dan ke bawah. Atau lebih tepatnya, kesalahan “melompat” tidak akan berhenti, tetapi tidak akan terjadi begitu cepat :) Mari kita periksa.

Kode untuk menjalankan SGD dengan penambahan 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()

Grafik No. 6 “Jumlah deviasi kuadrat selama penurunan gradien stokastik (80 ribu langkah)”

Koefisiennya sudah membaik, namun masih belum ideal. Secara hipotesis, hal ini dapat diperbaiki dengan cara ini. Kami memilih, misalnya, dalam 1000 iterasi terakhir nilai koefisien yang menghasilkan kesalahan minimum. Benar, untuk ini kita juga harus menuliskan nilai koefisiennya sendiri. Kami tidak akan melakukan ini, melainkan memperhatikan jadwalnya. Kelihatannya mulus dan errornya sepertinya berkurang secara merata. Sebenarnya, hal ini tidak benar. Mari kita lihat 1000 iterasi pertama dan bandingkan dengan yang terakhir.

Kode untuk grafik 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()

Grafik No. 7 “Jumlah deviasi kuadrat SGD (1000 langkah pertama)”

Grafik No.8 “Jumlah deviasi kuadrat SGD (1000 langkah terakhir)”

Pada awal penurunan, kami mengamati penurunan kesalahan yang cukup seragam dan tajam. Pada iterasi terakhir, kita melihat errornya terus menerus di sekitar nilai 1,475 dan di beberapa saat bahkan sama dengan nilai optimal tersebut, namun kemudian tetap naik... Saya ulangi, Anda bisa menuliskan nilai-nilainya. koefisien и , lalu pilih yang kesalahannya minimal. Namun, kami menghadapi masalah yang lebih serius: kami harus mengambil 80 ribu langkah (lihat kode) untuk mendapatkan nilai mendekati optimal. Dan ini sudah bertentangan dengan gagasan menghemat waktu komputasi dengan penurunan gradien stokastik dibandingkan dengan penurunan gradien. Apa yang bisa diperbaiki dan ditingkatkan? Tidak sulit untuk menyadari bahwa pada iterasi pertama kita dengan percaya diri turun dan, oleh karena itu, kita harus meninggalkan satu langkah besar pada iterasi pertama dan mengurangi langkah tersebut saat kita bergerak maju. Kami tidak akan melakukan ini di artikel ini - ini sudah terlalu panjang. Yang mau bisa mikir sendiri gimana caranya, nggak susah :)

Sekarang mari kita lakukan penurunan gradien stokastik menggunakan perpustakaan JumlahPy (dan jangan tersandung batu-batu yang kita identifikasi sebelumnya)

Kode untuk Penurunan Gradien Stochastic (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 dengan saat turun tanpa menggunakan JumlahPy. Namun, ini logis.

Mari kita cari tahu berapa lama waktu yang dibutuhkan untuk penurunan gradien stokastik.

Kode penentuan waktu perhitungan 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 awannya: sekali lagi, rumus yang “ditulis sendiri” menunjukkan hasil terbaik. Semua ini menunjukkan bahwa harus ada cara yang lebih halus dalam menggunakan perpustakaan JumlahPy, yang sangat mempercepat operasi komputasi. Pada artikel ini kita tidak akan mempelajarinya. Akan ada sesuatu untuk dipikirkan di waktu luang Anda :)

Ringkaslah

Sebelum menyimpulkan, saya ingin menjawab pertanyaan yang kemungkinan besar muncul dari pembaca yang budiman. Sebenarnya, mengapa “penyiksaan” seperti itu dilakukan dengan menuruni bukit, mengapa kita perlu berjalan naik turun gunung (kebanyakan turun) untuk menemukan dataran rendah yang berharga, jika kita memiliki perangkat yang begitu kuat dan sederhana di tangan kita? bentuk solusi analitis, yang langsung memindahkan kita ke tempat yang benar?

Jawaban atas pertanyaan ini ada di permukaan. Sekarang kita telah melihat contoh yang sangat sederhana, yang mana jawaban sebenarnya adalah tergantung pada satu tanda . Anda tidak sering melihat hal ini dalam hidup, jadi bayangkan kita memiliki 2, 30, 50 atau lebih tanda-tanda. Mari kita tambahkan ribuan, atau bahkan puluhan ribu nilai untuk setiap atribut. Dalam hal ini, solusi analitik mungkin tidak tahan terhadap pengujian dan gagal. Pada gilirannya, penurunan gradien dan variasinya perlahan tapi pasti akan membawa kita lebih dekat ke tujuan - fungsi minimum. Dan jangan khawatir tentang kecepatan - kita mungkin akan mencari cara yang memungkinkan kita mengatur dan mengatur panjang langkah (yaitu, kecepatan).

Dan sekarang ringkasan singkat sebenarnya.

Pertama, saya berharap materi yang disajikan dalam artikel ini dapat membantu “ilmuwan data” pemula dalam memahami cara menyelesaikan persamaan regresi linier sederhana (dan tidak hanya).

Kedua, kita melihat beberapa cara untuk menyelesaikan persamaan tersebut. Sekarang, tergantung situasinya, kita dapat memilih salah satu yang paling sesuai untuk menyelesaikan masalah.

Ketiga, kami melihat kekuatan pengaturan tambahan, yaitu panjang langkah penurunan gradien. Parameter ini tidak dapat diabaikan. Seperti disebutkan di atas, untuk mengurangi biaya perhitungan, panjang langkah harus diubah saat turun.

Keempat, dalam kasus kami, fungsi “yang ditulis sendiri” menunjukkan hasil waktu terbaik untuk perhitungan. Hal ini mungkin disebabkan oleh kurangnya penggunaan kemampuan perpustakaan secara profesional JumlahPy. Namun demikian, kesimpulan berikut ini menunjukkan dirinya sendiri. Di satu sisi, terkadang ada baiknya mempertanyakan pendapat yang sudah ada, dan di sisi lain, tidak selalu ada gunanya memperumit segalanya - sebaliknya, terkadang cara yang lebih sederhana untuk menyelesaikan suatu masalah lebih efektif. Dan karena tujuan kami adalah menganalisis tiga pendekatan untuk menyelesaikan persamaan regresi linier sederhana, penggunaan fungsi yang “ditulis sendiri” sudah cukup bagi kami.

Sastra (atau semacamnya)

1. Regresi linier

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

2. Metode kuadrat terkecil

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Turunan

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradien

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Penurunan gradien

habr.com/en/post/471458

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