Giải phương trình hồi quy tuyến tính đơn giản

Bài viết thảo luận về một số cách để xác định phương trình toán học của đường hồi quy đơn giản (theo cặp).

Tất cả các phương pháp giải phương trình được thảo luận ở đây đều dựa trên phương pháp bình phương tối thiểu. Hãy biểu thị các phương pháp như sau:

• Giải pháp phân tích
• Xuống dốc
• Giảm dần độ dốc ngẫu nhiên

Với mỗi phương pháp giải phương trình đường thẳng, bài viết cung cấp các hàm khác nhau, chủ yếu chia thành các hàm viết không sử dụng thư viện numpy và những thứ dùng để tính toán numpy. Người ta tin rằng việc sử dụng khéo léo numpy sẽ giảm chi phí tính toán.

Tất cả các mã đưa ra trong bài viết được viết bằng ngôn ngữ trăn 2.7 sử dụng Máy tính xách tay Jupyter. Mã nguồn và file dữ liệu mẫu được đăng tải trên Github

Bài viết hướng đến cả những người mới bắt đầu và những người đã dần dần nắm vững việc nghiên cứu một phần rất rộng trong trí tuệ nhân tạo - học máy.

Để minh họa tài liệu, chúng tôi sử dụng một ví dụ rất đơn giản.

Điều kiện ví dụ

Chúng tôi có năm giá trị đặc trưng cho sự phụ thuộc Y từ X (Bảng số 1):

Bảng số 1 “Điều kiện ví dụ”

Chúng ta sẽ giả định rằng các giá trị là tháng trong năm và - doanh thu trong tháng này. Nói cách khác, doanh thu phụ thuộc vào tháng trong năm và - dấu hiệu duy nhất mà doanh thu phụ thuộc vào.

Ví dụ này là như vậy, cả từ quan điểm về sự phụ thuộc có điều kiện của doanh thu vào tháng trong năm, lẫn từ quan điểm về số lượng giá trị - có rất ít trong số đó. Tuy nhiên, sự đơn giản hóa như vậy sẽ giúp bạn có thể giải thích, không phải lúc nào cũng dễ dàng, những tài liệu mà người mới bắt đầu tiếp thu. Và sự đơn giản của các con số sẽ cho phép những người muốn giải ví dụ trên giấy mà không tốn nhiều chi phí lao động.

Chúng ta hãy giả sử rằng sự phụ thuộc được đưa ra trong ví dụ có thể được xấp xỉ khá tốt bằng phương trình toán học của đường hồi quy đơn giản (theo cặp) có dạng:

đâu là tháng nhận được doanh thu, - doanh thu tương ứng với tháng, и là các hệ số hồi quy của đường ước lượng.

Lưu ý rằng hệ số thường được gọi là độ dốc hoặc độ dốc của đường ước lượng; đại diện cho số tiền mà khi nó thay đổi .

Rõ ràng, nhiệm vụ của chúng ta trong ví dụ là chọn các hệ số như vậy trong phương trình и , tại đó độ lệch của giá trị doanh thu được tính toán của chúng tôi theo tháng so với các câu trả lời đúng, tức là. các giá trị được trình bày trong mẫu sẽ ở mức tối thiểu.

Phương pháp bình phương tối thiểu

Theo phương pháp bình phương tối thiểu, độ lệch phải được tính bằng cách bình phương nó. Kỹ thuật này cho phép bạn tránh loại bỏ các sai lệch lẫn nhau nếu chúng có dấu hiệu trái ngược nhau. Ví dụ, nếu trong một trường hợp, độ lệch là +5 (cộng năm), và mặt khác -5 (trừ năm), khi đó tổng các độ lệch sẽ triệt tiêu lẫn nhau và bằng 0 (không). Có thể không bình phương độ lệch mà sử dụng tính chất của mô đun và khi đó tất cả các độ lệch sẽ dương và sẽ tích lũy. Chúng tôi sẽ không đi sâu vào điểm này một cách chi tiết mà chỉ chỉ ra rằng để thuận tiện cho việc tính toán, người ta thường bình phương độ lệch.

Đây là công thức mà chúng ta sẽ xác định tổng độ lệch bình phương (lỗi) nhỏ nhất:

đâu là hàm gần đúng của các câu trả lời đúng (nghĩa là doanh thu mà chúng tôi đã tính toán),

là những câu trả lời đúng (doanh thu được cung cấp trong mẫu),

là chỉ số mẫu (số tháng xác định độ lệch)

Hãy vi phân hàm số, xác định các phương trình vi phân từng phần và sẵn sàng chuyển sang giải pháp phân tích. Nhưng trước tiên, chúng ta hãy cùng tìm hiểu một chút về vi phân là gì và ghi nhớ ý nghĩa hình học của đạo hàm.

Sự khác biệt

Vi phân là phép toán tìm đạo hàm của một hàm số.

Đạo hàm dùng để làm gì? Đạo hàm của một hàm đặc trưng cho tốc độ thay đổi của hàm và cho chúng ta biết hướng của nó. Nếu đạo hàm tại một điểm cho trước là dương thì hàm số tăng; ngược lại hàm số sẽ giảm. Và giá trị đạo hàm tuyệt đối càng lớn thì tốc độ thay đổi của các giá trị hàm càng cao, cũng như độ dốc của đồ thị hàm số càng dốc.

Ví dụ, trong điều kiện của hệ tọa độ Descartes, giá trị đạo hàm tại điểm M(0,0) bằng +25 có nghĩa là tại một thời điểm nhất định, khi giá trị được dịch chuyển ở bên phải bởi đơn vị quy ước, giá trị tăng thêm 25 đơn vị thông thường. Trên biểu đồ, có vẻ như giá trị tăng khá nhanh từ một điểm nhất định.

Một vi dụ khac. Giá trị phái sinh bằng nhau -0,1 có nghĩa là khi di dời trên một đơn vị thông thường, giá trị chỉ giảm 0,1 đơn vị thông thường. Đồng thời, trên biểu đồ của hàm số, chúng ta có thể quan sát thấy độ dốc đi xuống hầu như không đáng chú ý. Hãy so sánh với một ngọn núi, nó giống như thể chúng ta đang đi xuống một con dốc thoai thoải từ một ngọn núi rất chậm, không giống như ví dụ trước, nơi chúng ta phải leo lên những đỉnh núi rất dốc :)

Như vậy, sau khi lấy đạo hàm theo tỷ lệ cược и , chúng tôi xác định phương trình vi phân từng phần bậc 1. Sau khi xác định các phương trình, chúng ta sẽ nhận được một hệ gồm hai phương trình, bằng cách giải chúng ta sẽ chọn được các giá trị như vậy của các hệ số и , tại đó giá trị của đạo hàm tương ứng tại các điểm đã cho thay đổi một lượng rất nhỏ và trong trường hợp nghiệm giải tích thì không thay đổi chút nào. Nói cách khác, hàm sai số tại các hệ số tìm được sẽ đạt mức tối thiểu, vì các giá trị của đạo hàm riêng tại các điểm này sẽ bằng 0.

Vì vậy, theo quy tắc đạo hàm riêng, phương trình đạo hàm riêng bậc 1 đối với hệ số sẽ có dạng:

phương trình đạo hàm riêng bậc 1 đối với sẽ có dạng:

Kết quả ta thu được hệ phương trình có nghiệm giải tích khá đơn giản:

bắt đầu{phương trình*}
bắt đầu{trường hợp}
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
kết thúc{trường hợp}
kết thúc{phương trình*}

Trước khi giải phương trình, hãy tải trước, kiểm tra xem việc tải có đúng không và định dạng dữ liệu.

Đang tải và định dạng dữ liệu

Cần lưu ý rằng do thực tế là đối với giải pháp phân tích và sau đó là độ dốc và độ dốc giảm độ dốc ngẫu nhiên, chúng tôi sẽ sử dụng mã theo hai biến thể: sử dụng thư viện numpy và không sử dụng nó thì chúng ta sẽ cần định dạng dữ liệu phù hợp (xem mã).

Mã tải và xử lý dữ liệu

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

Hình dung

Bây giờ, trước tiên, sau khi chúng ta tải dữ liệu, thứ hai, kiểm tra tính chính xác của quá trình tải và cuối cùng là định dạng dữ liệu, chúng ta sẽ tiến hành hiển thị đầu tiên. Phương pháp thường được sử dụng cho việc này là sự ghép đôi thư viện sơ sinh. Trong ví dụ của chúng tôi, do số lượng có hạn nên không có ích gì khi sử dụng thư viện sơ sinh. Chúng tôi sẽ sử dụng thư viện thông thường Matplotlib và chỉ cần nhìn vào biểu đồ phân tán.

Mã biểu đồ phân tán

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

Biểu đồ số 1 “Sự phụ thuộc của doanh thu vào các tháng trong năm”

Giải pháp phân tích

Hãy sử dụng những công cụ phổ biến nhất trong mãng xà và giải hệ phương trình:

bắt đầu{phương trình*}
bắt đầu{trường hợp}
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
kết thúc{trường hợp}
kết thúc{phương trình*}

Theo quy tắc Cramer chúng ta sẽ tìm định thức tổng quát, cũng như các định thức bằng cách và , sau đó chia định thức cho đến định thức tổng quát - tìm hệ số , tương tự ta tìm hệ số .

Mã giải pháp phân tích

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

Đây là những gì chúng tôi có:

Như vậy, các giá trị của các hệ số đã tìm được, tổng các bình phương đã được xác lập. Hãy vẽ một đường thẳng trên biểu đồ tán xạ theo các hệ số tìm được.

Mã dòng hồi quy

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

Biểu đồ số 2 “Câu trả lời đúng và có tính toán”

Bạn có thể nhìn vào biểu đồ độ lệch cho mỗi tháng. Trong trường hợp của chúng tôi, chúng tôi sẽ không thu được bất kỳ giá trị thực tế đáng kể nào từ nó, nhưng chúng tôi sẽ thỏa mãn sự tò mò của mình về phương trình hồi quy tuyến tính đơn giản mô tả sự phụ thuộc của doanh thu vào tháng trong năm như thế nào.

Mã biểu đồ độ lệch

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

Biểu đồ số 3 “Độ lệch, %”

Không hoàn hảo nhưng chúng tôi đã hoàn thành nhiệm vụ của mình.

Hãy viết hàm xác định các hệ số и sử dụng thư viện numpy, chính xác hơn, chúng ta sẽ viết hai hàm: một hàm sử dụng ma trận giả nghịch đảo (không được khuyến nghị trong thực tế, vì quá trình tính toán phức tạp và không ổn định), hàm kia sử dụng phương trình ma trận.

Mã giải pháp phân tích (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

Hãy so sánh thời gian dành cho việc xác định các hệ số и , theo 3 phương pháp đã trình bày.

Code tính thời gian tính toán

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)

Với một lượng nhỏ dữ liệu, một hàm “tự viết” sẽ xuất hiện, tìm ra các hệ số bằng phương pháp của Cramer.

Bây giờ bạn có thể chuyển sang các cách khác để tìm hệ số и .

Xuống dốc

Đầu tiên, hãy định nghĩa gradient là gì. Nói một cách đơn giản, gradient là một đoạn biểu thị hướng tăng trưởng tối đa của một hàm. Tương tự như việc leo núi, nơi các mặt dốc là nơi có độ dốc cao nhất để leo lên đỉnh núi. Phát triển ví dụ với ngọn núi, chúng ta nhớ rằng trên thực tế, chúng ta cần độ dốc cao nhất để đến vùng đất thấp càng nhanh càng tốt, tức là mức tối thiểu - nơi mà hàm không tăng hoặc giảm. Lúc này đạo hàm sẽ bằng 0. Do đó, chúng ta không cần gradient mà cần antigradient. Để tìm phản gradient bạn chỉ cần nhân gradient với -1 (trừ đi một).

Chúng ta hãy chú ý đến thực tế là một hàm có thể có một số cực tiểu và khi đi xuống một trong số chúng bằng thuật toán được đề xuất bên dưới, chúng ta sẽ không thể tìm thấy mức tối thiểu khác, có thể thấp hơn mức tối thiểu được tìm thấy. Hãy thư giãn, đây không phải là mối đe dọa đối với chúng tôi! Trong trường hợp của chúng ta, chúng ta đang xử lý một mức tối thiểu duy nhất, vì hàm của chúng ta trên đồ thị là một parabol đều. Và như tất cả chúng ta đều biết rất rõ từ khóa học toán ở trường, một parabol chỉ có một mức tối thiểu.

Sau khi chúng tôi tìm ra lý do tại sao chúng tôi cần một gradient và gradient đó cũng là một đoạn, nghĩa là một vectơ có tọa độ cho trước, có cùng hệ số и chúng ta có thể thực hiện giảm độ dốc.

Trước khi bắt đầu, tôi khuyên bạn chỉ nên đọc một vài câu về thuật toán giảm dần:

• Chúng tôi xác định một cách giả ngẫu nhiên tọa độ của các hệ số и . Trong ví dụ của chúng tôi, chúng tôi sẽ xác định các hệ số gần bằng 0. Đây là một thực tế phổ biến, nhưng mỗi trường hợp có thể có cách thực hành riêng.
• Từ tọa độ trừ giá trị đạo hàm riêng cấp 1 tại điểm . Vì vậy, nếu đạo hàm dương thì hàm số tăng. Do đó, bằng cách trừ đi giá trị của đạo hàm, chúng ta sẽ di chuyển theo hướng tăng trưởng ngược lại, tức là theo hướng đi xuống. Nếu đạo hàm âm thì hàm số tại thời điểm này sẽ giảm và bằng cách trừ đi giá trị của đạo hàm, chúng ta sẽ di chuyển theo hướng giảm dần.
• Chúng tôi thực hiện một hoạt động tương tự với tọa độ : trừ giá trị đạo hàm riêng tại điểm .
• Để không nhảy quá mức tối thiểu và bay vào không gian sâu, cần đặt kích thước bước theo hướng đi xuống. Nói chung, bạn có thể viết cả một bài viết về cách thiết lập bước chính xác và cách thay đổi bước đó trong quá trình giảm dần để giảm chi phí tính toán. Nhưng bây giờ chúng ta có một nhiệm vụ hơi khác trước mắt và chúng ta sẽ thiết lập kích thước bước bằng phương pháp khoa học là “chọc” hoặc, như người ta nói theo cách nói thông thường, theo kinh nghiệm.
• Khi chúng ta đến từ tọa độ đã cho и trừ đi các giá trị đạo hàm ta được tọa độ mới и . Chúng tôi thực hiện bước tiếp theo (trừ), từ tọa độ được tính toán. Và cứ thế chu trình bắt đầu lại nhiều lần cho đến khi đạt được sự hội tụ cần thiết.

Tất cả! Bây giờ chúng ta đã sẵn sàng lên đường tìm kiếm hẻm núi sâu nhất của rãnh Mariana. Bắt đầu nào.

Mã giảm độ dốc

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

Chúng tôi đã lặn xuống đáy rãnh Mariana và ở đó chúng tôi tìm thấy tất cả các giá trị hệ số giống nhau и , đó chính xác là những gì được mong đợi.

Chúng ta hãy lặn lần nữa, chỉ lần này, phương tiện dưới biển sâu của chúng ta sẽ chứa đầy những công nghệ khác, cụ thể là thư viện numpy.

Mã để giảm độ dốc (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

Giá trị hệ số и không thể thay đổi.

Chúng ta hãy xem lỗi đã thay đổi như thế nào trong quá trình giảm độ dốc, tức là tổng độ lệch bình phương thay đổi như thế nào sau mỗi bước.

Mã để vẽ tổng các độ lệch bình phương

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

Đồ thị số 4 “Tổng bình phương độ lệch khi giảm độ dốc”

Trên biểu đồ, chúng ta thấy rằng với mỗi bước, sai số sẽ giảm dần và sau một số lần lặp nhất định, chúng ta quan sát thấy một đường gần như nằm ngang.

Cuối cùng, hãy ước tính sự khác biệt về thời gian thực thi mã:

Mã xác định thời gian tính toán độ dốc giảm dần

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)

Có lẽ chúng ta đang làm sai điều gì đó, nhưng một lần nữa, đó lại là một hàm “tự viết” đơn giản không sử dụng thư viện numpy vượt trội hơn thời gian tính toán của hàm bằng thư viện numpy.

Nhưng chúng tôi không đứng yên mà đang hướng tới việc nghiên cứu một cách thú vị khác để giải phương trình hồi quy tuyến tính đơn giản. Gặp!

Giảm dần độ dốc ngẫu nhiên

Để hiểu nhanh nguyên lý hoạt động của quá trình giảm độ dốc ngẫu nhiên, tốt hơn là xác định sự khác biệt của nó so với việc giảm độ dốc thông thường. Chúng ta, trong trường hợp giảm độ dốc, trong các phương trình đạo hàm của и đã sử dụng tổng giá trị của tất cả các đặc điểm và câu trả lời đúng có sẵn trong mẫu (nghĩa là tổng của tất cả и ). Trong quá trình giảm độ dốc ngẫu nhiên, chúng tôi sẽ không sử dụng tất cả các giá trị có trong mẫu mà thay vào đó, chọn giả ngẫu nhiên cái gọi là chỉ mục mẫu và sử dụng các giá trị của nó.

Ví dụ: nếu chỉ số được đánh số 3 (ba) thì ta lấy các giá trị и , sau đó ta thay các giá trị vào phương trình đạo hàm và xác định tọa độ mới. Sau đó, khi xác định được tọa độ, chúng ta lại xác định giả ngẫu nhiên chỉ số mẫu, thay thế các giá trị tương ứng với chỉ số vào các phương trình vi phân từng phần và xác định tọa độ theo cách mới и vân vân. cho đến khi hội tụ chuyển sang màu xanh. Thoạt nhìn, có vẻ như điều này không thể thực hiện được, nhưng thực tế là có. Đúng là điều đáng chú ý là sai số không giảm theo từng bước nhưng chắc chắn có xu hướng xảy ra.

Ưu điểm của phương pháp giảm độ dốc ngẫu nhiên so với phương pháp thông thường là gì? Nếu kích thước mẫu của chúng tôi rất lớn và được đo bằng hàng chục nghìn giá trị, thì việc xử lý một nghìn giá trị ngẫu nhiên trong số đó sẽ dễ dàng hơn nhiều thay vì toàn bộ mẫu. Đây là lúc độ dốc giảm dần ngẫu nhiên phát huy tác dụng. Tất nhiên, trong trường hợp của chúng tôi, chúng tôi sẽ không nhận thấy nhiều sự khác biệt.

Hãy nhìn vào mã.

Mã giảm dần độ dốc ngẫu nhiên

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

Chúng ta xem xét cẩn thận các hệ số và tự đặt ra câu hỏi “Làm sao điều này có thể xảy ra được?” Chúng tôi có các giá trị hệ số khác и . Có lẽ việc giảm độ dốc ngẫu nhiên đã tìm thấy các tham số tối ưu hơn cho phương trình? Tiếc là không có. Chỉ cần nhìn vào tổng độ lệch bình phương là đủ và thấy rằng với các giá trị mới của các hệ số, sai số sẽ lớn hơn. Chúng tôi không vội tuyệt vọng. Hãy xây dựng một biểu đồ về sự thay đổi lỗi.

Mã để vẽ tổng các bình phương độ lệch khi giảm độ dốc ngẫu nhiên

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

Đồ thị số 5 “Tổng bình phương độ lệch khi giảm độ dốc ngẫu nhiên”

Nhìn vào lịch trình, mọi thứ đã ổn định và bây giờ chúng tôi sẽ sửa chữa mọi thứ.

Vậy chuyện gì đã xảy ra? Điều sau đây đã xảy ra. Khi chúng tôi chọn ngẫu nhiên một tháng, thì thuật toán của chúng tôi sẽ tìm cách giảm lỗi trong tháng được chọn trong tháng đã chọn. Sau đó, chúng tôi chọn một tháng khác và lặp lại phép tính, nhưng chúng tôi sẽ giảm sai số cho tháng được chọn thứ hai. Bây giờ chúng ta hãy nhớ rằng hai tháng đầu tiên có sự khác biệt đáng kể so với phương trình hồi quy tuyến tính đơn giản. Điều này có nghĩa là khi bất kỳ tháng nào trong hai tháng này được chọn, bằng cách giảm sai số của từng tháng, thuật toán của chúng tôi sẽ làm tăng sai số cho toàn bộ mẫu một cách nghiêm trọng. Vậy lam gi? Câu trả lời rất đơn giản: bạn cần giảm bước đi xuống. Rốt cuộc, bằng cách giảm bước đi xuống, lỗi cũng sẽ ngừng “nhảy” lên xuống. Hay đúng hơn là lỗi “nhảy” sẽ không dừng lại, nhưng nó sẽ không xảy ra nhanh như vậy :) Hãy kiểm tra.

Mã để chạy SGD với số gia nhỏ hơn

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

Biểu đồ số 6 “Tổng bình phương độ lệch trong quá trình giảm độ dốc ngẫu nhiên (80 nghìn bước)”

Các hệ số đã được cải thiện nhưng vẫn chưa lý tưởng. Theo giả thuyết, điều này có thể được sửa chữa theo cách này. Ví dụ, chúng tôi chọn trong 1000 lần lặp cuối cùng các giá trị của các hệ số có sai số tối thiểu. Đúng, để làm được điều này, chúng ta cũng sẽ phải viết ra giá trị của các hệ số. Chúng tôi sẽ không làm điều này mà chỉ chú ý đến lịch trình. Nó trông mượt mà và lỗi dường như giảm đều. Trên thực tế, điều này là không đúng sự thật. Hãy nhìn vào 1000 lần lặp đầu tiên và so sánh chúng với lần lặp trước.

Mã cho biểu đồ SGD (1000 bước đầu tiên)

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

Đồ thị số 7 “Tổng bình phương độ lệch SGD (1000 bước đầu tiên)”

Đồ thị số 8 “Tổng bình phương độ lệch SGD (1000 bước cuối)”

Khi bắt đầu đi xuống, chúng tôi quan sát thấy sai số giảm khá đồng đều và mạnh. Trong các lần lặp cuối cùng, chúng ta thấy rằng lỗi xoay quanh giá trị 1,475 và tại một số thời điểm thậm chí bằng giá trị tối ưu này, nhưng sau đó nó vẫn tăng lên... Tôi nhắc lại, bạn có thể viết ra các giá trị của hệ số и , rồi chọn những lỗi có lỗi tối thiểu. Tuy nhiên, chúng tôi gặp phải một vấn đề nghiêm trọng hơn: chúng tôi phải thực hiện 80 nghìn bước (xem mã) để đạt được giá trị gần mức tối ưu. Và điều này đã mâu thuẫn với ý tưởng tiết kiệm thời gian tính toán bằng cách giảm độ dốc ngẫu nhiên so với giảm độ dốc. Những gì có thể được sửa chữa và cải thiện? Không khó để nhận thấy rằng trong những lần lặp đầu tiên, chúng ta đang tự tin đi xuống và do đó, chúng ta nên để lại một bước tiến lớn trong những lần lặp đầu tiên và giảm bước khi tiến về phía trước. Chúng tôi sẽ không làm điều này trong bài viết này - nó đã quá dài rồi. Ai có nhu cầu có thể tự suy nghĩ cách làm, không khó đâu :)

Bây giờ hãy thực hiện giảm độ dốc ngẫu nhiên bằng thư viện numpy (và đừng vấp phải những viên đá mà chúng ta đã xác định trước đó)

Mã cho Giảm dần độ dốc ngẫu nhiên (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

Các giá trị hóa ra gần giống như khi giảm dần mà không sử dụng numpy. Tuy nhiên, điều này là hợp lý.

Hãy cùng tìm hiểu xem chúng ta mất bao lâu để giảm độ dốc ngẫu nhiên.

Mã xác định thời gian tính SGD (80 nghìn bước)

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)

Càng đi sâu vào rừng, mây càng đen: một lần nữa, công thức “tự viết” cho kết quả tốt nhất. Tất cả điều này gợi ý rằng phải có những cách tinh tế hơn nữa để sử dụng thư viện numpy, điều này thực sự tăng tốc các hoạt động tính toán. Trong bài viết này chúng ta sẽ không tìm hiểu về chúng. Sẽ có điều gì đó để bạn suy nghĩ lúc rảnh rỗi :)

Tóm tắt

Trước khi tóm tắt, tôi muốn trả lời một câu hỏi rất có thể được đặt ra từ độc giả thân yêu của chúng tôi. Trên thực tế, tại sao lại phải “tra tấn” bằng những đoạn đường dốc như vậy, tại sao chúng ta phải đi bộ lên xuống núi (chủ yếu là đi xuống) để tìm được vùng đất trũng quý giá, nếu chúng ta có trong tay một thiết bị mạnh mẽ và đơn giản như vậy, trong một giải pháp phân tích sẽ đưa chúng ta ngay lập tức đến Đúng nơi?

Câu trả lời cho câu hỏi này nằm trên bề mặt. Bây giờ chúng ta đã xem xét một ví dụ rất đơn giản, trong đó câu trả lời đúng là phụ thuộc vào một dấu hiệu . Bạn không thường thấy điều này trong cuộc sống, vì vậy hãy tưởng tượng rằng chúng ta có 2, 30, 50 hoặc nhiều dấu hiệu hơn. Hãy thêm vào hàng nghìn, thậm chí hàng chục nghìn giá trị cho mỗi thuộc tính. Trong trường hợp này, dung dịch phân tích có thể không chịu được thử nghiệm và thất bại. Đổi lại, việc giảm độ dốc và các biến thể của nó sẽ chậm rãi nhưng chắc chắn sẽ đưa chúng ta đến gần hơn với mục tiêu - mức tối thiểu của hàm. Và đừng lo lắng về tốc độ - có thể chúng ta sẽ xem xét các cách cho phép chúng ta thiết lập và điều chỉnh độ dài bước (tức là tốc độ).

Và bây giờ là bản tóm tắt ngắn gọn thực tế.

Đầu tiên, tôi hy vọng rằng tài liệu được trình bày trong bài viết sẽ giúp những người mới bắt đầu làm “nhà khoa học dữ liệu” hiểu cách giải các phương trình hồi quy tuyến tính đơn giản (và không chỉ).

Thứ hai, chúng tôi đã xem xét một số cách để giải phương trình. Bây giờ, tùy theo tình huống mà chúng ta có thể chọn cách phù hợp nhất để giải quyết vấn đề.

Thứ ba, chúng tôi thấy được sức mạnh của các cài đặt bổ sung, cụ thể là độ dài bước giảm độ dốc. Thông số này không thể bỏ qua. Như đã lưu ý ở trên, để giảm chi phí tính toán, độ dài bước cần được thay đổi trong quá trình giảm dần.

Thứ tư, trong trường hợp của chúng tôi, các hàm “viết tại nhà” cho kết quả tính toán theo thời gian tốt nhất. Điều này có lẽ là do việc sử dụng các khả năng của thư viện chưa chuyên nghiệp nhất numpy. Nhưng dù có thể như vậy thì kết luận sau đây cũng tự nó gợi ý. Một mặt, đôi khi cần đặt câu hỏi về những quan điểm đã có sẵn, mặt khác, không phải lúc nào cũng đáng để phức tạp hóa mọi thứ - ngược lại, đôi khi cách giải quyết vấn đề đơn giản hơn sẽ hiệu quả hơn. Và vì mục tiêu của chúng tôi là phân tích ba cách tiếp cận để giải một phương trình hồi quy tuyến tính đơn giản nên việc sử dụng các hàm “tự viết” là khá đủ đối với chúng tôi.

Văn học (hoặc một cái gì đó như thế)

1. Hồi quy tuyến tính

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

2. Phương pháp bình phương tối thiểu

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Phái sinh

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Độ dốc

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Độ dốc gốc

habr.com/vi/post/471458

habr.com/vi/post/307312

artemarakcheev.com//2017-12-31/Tuyến_regression

6. Thư viện 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

Nguồn: www.habr.com