حل معادلة انحدار خطي بسيطة

تتناول المقالة عدة طرق لتحديد المعادلة الرياضية لخط انحدار بسيط (زوج).

تعتمد جميع طرق حل المعادلة المذكورة هنا على طريقة المربعات الصغرى. نشير إلى الطرق على النحو التالي:

الحل التحليلي
نزول متدرج
الانحدار العشوائي

لكل طريقة من طرق حل معادلة الخط المستقيم ، تقدم المقالة وظائف مختلفة ، والتي تنقسم بشكل أساسي إلى تلك التي تمت كتابتها دون استخدام المكتبة نمباي وتلك التي تستخدم في الحسابات نمباي. ويعتقد أن ماهرا في الاستخدام نمباي سيقلل من تكلفة الحوسبة.

تمت كتابة كل التعليمات البرمجية في هذه المقالة بيثون 2.7 مع مفكرة Jupyter. يتوفر كود المصدر وملف البيانات النموذجي في جيثب

تركز المقالة بشكل أكبر على كل من المبتدئين وأولئك الذين بدأوا بالفعل تدريجيًا في إتقان دراسة قسم واسع جدًا في الذكاء الاصطناعي - التعلم الآلي.

دعنا نستخدم مثالًا بسيطًا جدًا لتوضيح المادة.

مثال على الشروط

لدينا خمس قيم تميز الإدمان Y من X (الجدول 1):

الجدول رقم 1 "شروط المثال"

سنفترض أن القيم هو شهر السنة ، و - أرباح هذا الشهر. بمعنى آخر ، تعتمد الإيرادات على شهر السنة ، و - العلامة الوحيدة التي تعتمد عليها الإيرادات.

هذا المثال هو كذلك ، سواء من حيث الاعتماد المشروط للإيرادات على شهر من السنة ، ومن حيث عدد القيم - هناك عدد قليل جدًا منهم. ومع ذلك ، سيسمح مثل هذا التبسيط ، كما يقولون على الأصابع ، بشرح ، ليس دائمًا بسهولة ، المواد التي استوعبها المبتدئين. وأيضًا ، ستسمح بساطة الأرقام لأولئك الذين يرغبون في حل المثال على "الورق" دون تكاليف العمالة الكبيرة.

افترض أن الاعتماد الوارد في المثال يمكن تقريبه جيدًا بالمعادلة الرياضية لخط انحدار بسيط (زوج) للنموذج:

حيث هو الشهر الذي تم فيه استلام العائدات ، - الإيرادات المقابلة للشهر ، и هي معاملات الانحدار للخط المقدر.

لاحظ أن المعامل غالبًا ما يشار إليه على أنه منحدر أو انحدار الخط المقدر ؛ هو المبلغ الذي عندما يتغير .

من الواضح أن مهمتنا في المثال هي اختيار مثل هذه المعاملات في المعادلة и ، حيث تكون انحرافات قيم الإيرادات المقدرة لدينا حسب الأشهر عن الإجابات الصحيحة ، أي القيم المقدمة في العينة ستكون ضئيلة.

طريقة المربعات الصغرى

وفقًا لطريقة المربعات الصغرى ، يجب حساب الانحراف عن طريق تربيعه. تسمح هذه التقنية بتجنب السداد المتبادل للانحرافات ، إذا كانت لها علامات معاكسة. على سبيل المثال ، إذا كان الانحراف في حالة واحدة +5 (زائد خمسة) ، والآخر -5 (ناقص خمسة) ، فسيتم إلغاء مجموع الانحرافات بشكل متبادل وسيكون 0 (صفر). من الممكن عدم تربيع الانحراف ، ولكن استخدام خاصية المقياس ، وبعد ذلك ستكون جميع الانحرافات موجبة وسوف تتراكم. لن نتطرق إلى هذه النقطة بالتفصيل ، لكننا ببساطة نشير إلى أنه من أجل تسهيل العمليات الحسابية ، من المعتاد تربيع الانحراف.

هذه هي الطريقة التي تبدو بها الصيغة ، والتي من خلالها سنحدد أصغر مجموع مربعات الانحرافات (الأخطاء):

حيث هي دالة لتقريب الإجابات الصحيحة (أي الإيرادات التي نحسبها) ،

هي الإجابات الصحيحة (الإيرادات الواردة في العينة) ،

هو مؤشر العينة (رقم الشهر الذي يتم فيه تحديد الانحراف)

دعونا نفرق بين الدالة ، ونحدد المعادلات التفاضلية الجزئية ، ونستعد للانتقال إلى حل تحليلي. لكن أولاً ، لنأخذ استطراداً قصيرًا حول ماهية الاشتقاق ونتذكر المعنى الهندسي للمشتقة.

التفاضل

التفاضل هو عملية إيجاد مشتق دالة.

ما هو المشتق ل؟ يميز مشتق الدالة معدل تغير الوظيفة ويظهر لنا اتجاهها. إذا كان المشتق عند نقطة معينة موجبًا ، فإن الدالة تتزايد ؛ وإلا فإن الدالة تتناقص. وكلما زادت قيمة مشتق modulo ، زاد معدل تغير قيم الدالة ، كما زاد انحدار الرسم البياني للوظيفة.

على سبيل المثال ، في ظل ظروف نظام الإحداثيات الديكارتية ، فإن قيمة المشتق عند النقطة M (0,0،XNUMX) تساوي 25 يعني أنه عند نقطة معينة ، عندما يتم إزاحة القيمة إلى اليمين بواسطة وحدة تقليدية ، القيمة يزيد بمقدار 25 وحدة تقليدية. على الرسم البياني ، يبدو أنه زاوية ارتفاع حادة إلى حد ما في القيم من نقطة معينة.

مثال آخر. قيمة المشتق -0,1 يعني ذلك عند التحول لكل وحدة تقليدية واحدة ، القيمة ينخفض بنسبة 0,1 وحدة تقليدية فقط. في الوقت نفسه ، على الرسم البياني للدالة ، يمكننا أن نلاحظ ميلًا هبوطيًا بالكاد ملحوظًا. عند رسم تشبيه بجبل ، يبدو الأمر كما لو كنا ننزلق ببطء شديد من منحدر لطيف من جبل ، على عكس المثال السابق ، حيث كان علينا أن نأخذ قمم شديدة الانحدار :)

وهكذا ، بعد التفريق بين الوظيفة عن طريق الصعاب и ، نحدد معادلات المشتقات الجزئية من الدرجة الأولى. بعد تحديد المعادلتين ، سنحصل على نظام من معادلتين ، يمكننا حلهما من خلال اختيار هذه القيم للمعاملات и ، حيث تتغير قيم المشتقات المقابلة عند نقاط معينة بمقدار صغير جدًا جدًا ، وفي حالة الحل التحليلي لا تتغير على الإطلاق. بمعنى آخر ، ستصل دالة الخطأ عند المعامِلات الموجودة إلى الحد الأدنى ، لأن قيم المشتقات الجزئية عند هذه النقاط ستساوي صفرًا.

لذلك ، وفقًا لقواعد التفاضل ، فإن معادلة المشتق الجزئي من الدرجة الأولى فيما يتعلق بالمعامل سوف يأخذ النموذج:

معادلة مشتقة جزئية من الدرجة الأولى فيما يتعلق بـ سوف يأخذ النموذج:

نتيجة لذلك ، حصلنا على نظام معادلات يحتوي على حل تحليلي بسيط إلى حد ما:

ابدأ {المعادلة *}
تبدأ {الحالات}
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
إنهاء {الحالات}
نهاية {المعادلة *}

قبل حل المعادلة ، دعنا نقوم بالتحميل المسبق ، والتحقق من صحة التحميل وتنسيق البيانات.

تحميل وتنسيق البيانات

تجدر الإشارة إلى أنه نظرًا لحقيقة أنه بالنسبة للحل التحليلي ، وفي المستقبل لنسب التدرج اللوني العشوائي ، سنستخدم الكود في شكلين مختلفين: استخدام المكتبة نمباي وبدون استخدامه ، نحتاج إلى تنسيق البيانات المناسب (انظر الكود).

كود تحميل ومعالجة البيانات

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

تصور

الآن ، بعد أن قمنا أولاً بتحميل البيانات ، وثانياً ، تحققنا من صحة التحميل وقمنا أخيرًا بتنسيق البيانات ، سنقوم بتنفيذ التصور الأول. غالبًا ما يتم استخدام هذه الطريقة زوج مكتبة سيبورن. في مثالنا ، نظرًا لمحدودية الأرقام ، لا معنى لاستخدام المكتبة سيبورن. سوف نستخدم المكتبة المعتادة ماتبلوتليب وانظر فقط إلى مخطط التشتت.

كود مبعثر

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

الرسم البياني رقم 1 "اعتماد الإيرادات على شهر من السنة"

الحل التحليلي

نحن نستخدم الأدوات الأكثر شيوعًا في الثعبان وحل نظام المعادلات:

ابدأ {المعادلة *}
تبدأ {الحالات}
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
إنهاء {الحالات}
نهاية {المعادلة *}

حسب قاعدة كرامر العثور على المحدد المشترك ، وكذلك المحددات بواسطة و وبعد ذلك قسمة المحدد على لمحدد مشترك - أوجد المعامل ، وبالمثل نجد المعامل .

كود الحل التحليلي

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

هذا ما حصلنا عليه:

لذلك ، تم العثور على قيم المعاملات ، يتم تعيين مجموع الانحرافات التربيعية. لنرسم خطًا مستقيمًا على الرسم البياني للتشتت وفقًا للمعاملات الموجودة.

كود خط الانحدار

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

مخطط رقم 2 "الإجابات الصحيحة والمحسوبة"

يمكنك إلقاء نظرة على الرسم البياني للانحرافات لكل شهر. في حالتنا ، لن نستمد أي قيمة عملية كبيرة منها ، لكننا سنرضي الفضول في مدى جودة معادلة الانحدار الخطي البسيطة التي تميز اعتماد الإيرادات على شهر من العام.

رمز مخطط الانحراف

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

مخطط رقم 3 "الانحرافات ،٪"

ليست مثالية ، لكننا قمنا بعملنا.

لنكتب دالة لتحديد المعاملات и يستخدم المكتبة نمباي، بشكل أكثر دقة ، سنكتب وظيفتين: الأولى تستخدم مصفوفة عكسية زائفة (غير موصى بها في الممارسة العملية ، لأن العملية معقدة حسابيًا وغير مستقرة) ، والأخرى تستخدم معادلة مصفوفة.

كود الحل التحليلي (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

قارن الوقت المستغرق لتحديد المعاملات и ، وفقًا للطرق الثلاثة المقدمة.

رمز لحساب وقت الحساب

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)

على كمية صغيرة من البيانات ، تظهر وظيفة "مكتوبة ذاتيًا" ، والتي تجد المعاملات باستخدام طريقة كرامر.

يمكنك الآن الانتقال إلى طرق أخرى لإيجاد المعاملات и .

نزول متدرج

أولاً ، دعنا نحدد ما هو التدرج اللوني. بطريقة بسيطة ، التدرج اللوني عبارة عن مقطع يشير إلى اتجاه أقصى نمو للدالة. عن طريق القياس مع التسلق صعودًا ، حيث يبدو التدرج ، يوجد أكبر تسلق إلى قمة الجبل. عند تطوير مثال الجبل ، تذكر أننا في الواقع نحتاج إلى أقصى نزول من أجل الوصول إلى أدنى مستوى بأسرع ما يمكن ، أي الحد الأدنى - المكان الذي لا تزيد فيه الوظيفة ولا تنقص. عند هذه النقطة ، سيكون المشتق صفرًا. لذلك ، لا نحتاج إلى التدرج اللوني ، بل نحتاج إلى التدرج المضاد. للعثور على مضاد التدرج ، تحتاج فقط إلى ضرب التدرج في -1 (ناقص واحد).

دعنا ننتبه إلى حقيقة أن الوظيفة يمكن أن يكون لها العديد من الحدود الدنيا ، وبعد أن تنزل إلى واحدة منها وفقًا للخوارزمية المقترحة أدناه ، لن نتمكن من إيجاد حد أدنى آخر ، قد يكون أقل من الموجود. استرخي ، نحن لسنا في خطر! في حالتنا ، نحن نتعامل مع حد أدنى واحد ، منذ وظيفتنا على الرسم البياني هو القطع المكافئ العادي. وكما يجب أن نعلم جميعًا جيدًا من دورة الرياضيات المدرسية ، فإن القطع المكافئ له حد أدنى واحد فقط.

بعد أن اكتشفنا سبب حاجتنا إلى التدرج اللوني ، وكذلك أن التدرج هو مقطع ، أي متجه بإحداثيات معينة ، وهي نفس المعاملات فقط и يمكننا تنفيذ نزول متدرج.

قبل البدء ، أقترح قراءة بضع جمل فقط حول خوارزمية النسب:

نحدد بشكل عشوائي شبه عشوائي إحداثيات المعاملات и . في مثالنا ، سنحدد المعاملات بالقرب من الصفر. هذه ممارسة شائعة ، ولكن قد يكون لكل حالة ممارسة خاصة بها.
من تنسيق اطرح قيمة المشتق الجزئي من الدرجة الأولى عند النقطة . لذلك ، إذا كان المشتق موجبًا ، فإن الدالة تتزايد. لذلك ، بطرح قيمة المشتق ، سوف نتحرك في الاتجاه المعاكس للنمو ، أي في اتجاه الهبوط. إذا كان المشتق سالبًا ، فإن الوظيفة عند هذه النقطة تتناقص وتطرح قيمة المشتق ، نتحرك نحو الهبوط.
نقوم بعملية مماثلة مع الإحداثيات : اطرح قيمة المشتق الجزئي عند النقطة .
لكي لا تقفز فوق الحد الأدنى ولا تطير بعيدًا في الفضاء السحيق ، من الضروري ضبط حجم الخطوة في اتجاه الهبوط. بشكل عام ، يمكن للمرء أن يكتب مقالة كاملة حول كيفية ضبط الخطوة بشكل صحيح وكيفية تغييرها أثناء الهبوط لتقليل تكلفة العمليات الحسابية. لكن الآن لدينا مهمة مختلفة قليلاً ، وسوف نحدد حجم الخطوة بالطريقة العلمية "الوخز" أو ، كما يقولون في عامة الناس ، تجريبيًا.
بمجرد خروجنا من الإحداثيات المعطاة и اطرح قيم المشتقات ، نحصل على إحداثيات جديدة и . نأخذ الخطوة التالية (الطرح) ، بالفعل من الإحداثيات المحسوبة. وهكذا تبدأ الدورة مرارًا وتكرارًا ، حتى الوصول إلى نقطة التقارب المطلوبة.

الجميع! نحن الآن على استعداد للذهاب للبحث عن أعمق ممر خندق ماريانا. هيا بنا نبدأ.

رمز نزول التدرج

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

لقد غوصنا في قاع خندق ماريانا وهناك وجدنا جميع القيم نفسها للمعاملات и وهو أمر متوقع بالفعل.

دعنا نغطس مرة أخرى ، هذه المرة فقط ، حشو مركبتنا في أعماق البحار ستكون تقنيات أخرى ، وهي المكتبة نمباي.

رمز نزول التدرج (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

قيم المعامل и لم تتغير.

لنلقِ نظرة على كيفية تغير الخطأ أثناء نزول التدرج ، أي كيف تغير مجموع الانحرافات التربيعية مع كل خطوة.

رمز مخطط انحراف المجموع التربيعي

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

الرسم البياني رقم 4 "مجموع الانحرافات التربيعية في الانحدار المتدرج"

على الرسم البياني ، نرى أن الخطأ يتناقص مع كل خطوة ، وبعد عدد معين من التكرارات ، نلاحظ خطًا أفقيًا تقريبًا.

أخيرًا ، دعنا نقيم الفرق في وقت تنفيذ الكود:

كود لحساب نزول التدرج الزمني

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)

ربما نفعل شيئًا خاطئًا ، ولكن مرة أخرى وظيفة بسيطة "مكتوبة ذاتيًا" لا تستخدم المكتبة نمباي قبل وقت حساب الوظيفة باستخدام المكتبة نمباي.

لكننا لا نقف مكتوفي الأيدي ، لكننا نتجه نحو تعلم طريقة أخرى مثيرة لحل معادلة الانحدار الخطي البسيطة. يقابل!

الانحدار العشوائي

من أجل فهم كيفية عمل النسب المتدرج العشوائي بسرعة ، من الأفضل تحديد اختلافاته عن نزول التدرج التقليدي. نحن ، في حالة الانحدار ، في معادلات مشتقات и استخدم مجموع قيم جميع الميزات والإجابات الصحيحة المتوفرة في العينة (أي مجموع الكل и ). في نزول التدرج العشوائي ، لن نستخدم جميع القيم الموجودة في العينة ، ولكن بدلاً من ذلك ، سنختار بشكل شبه عشوائي ما يسمى بمؤشر العينة ونستخدم قيمه.

على سبيل المثال ، إذا تم تعريف الفهرس على أنه رقم 3 (ثلاثة) ، فإننا نأخذ القيم и ، ثم نعوض بالقيم في معادلات المشتقات ونحدد إحداثيات جديدة. بعد ذلك ، بعد تحديد الإحداثيات ، نحدد مرة أخرى مؤشر العينة بطريقة شبه عشوائية ، ونستبدل القيم المقابلة للفهرس في المعادلات التفاضلية الجزئية ، ونحدد الإحداثيات بطريقة جديدة и إلخ. حتى التقارب الأخضر. للوهلة الأولى ، قد يبدو أنه يمكن أن يعمل على الإطلاق ، لكنه يعمل. صحيح ، تجدر الإشارة إلى أن الخطأ لا ينقص مع كل خطوة ، ولكن هناك بالتأكيد اتجاه.

ما هي مزايا الانحدار العشوائي على الانحدار التقليدي؟ إذا كان حجم العينة كبيرًا جدًا ويتم قياسه بعشرات الآلاف من القيم ، فسيكون من الأسهل بكثير معالجة ، على سبيل المثال ، ألف منهم بشكل عشوائي ، من العينة بأكملها. هذا هو المكان الذي يبدأ فيه هبوط التدرج العشوائي. في حالتنا ، بالطبع ، لن نلاحظ فرقًا كبيرًا.

لنلق نظرة على الكود.

رمز هبوط التدرج العشوائي

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

نحن ننظر بتمعن في المعاملات وندرك أنفسنا في السؤال "كيف ذلك؟". حصلنا على قيم أخرى للمعاملات и . ربما وجد نزول التدرج العشوائي المزيد من المعلمات المثلى للمعادلة؟ للاسف لا. يكفي أن ننظر إلى مجموع الانحرافات التربيعية ونرى أنه مع القيم الجديدة للمعاملات ، يكون الخطأ أكبر. نحن لسنا في عجلة من أمرنا لليأس. لنقم ببناء رسم بياني لتغيير الخطأ.

رمز لتخطيط مجموع الانحرافات التربيعية في أصل التدرج العشوائي

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

الرسم البياني رقم 5 "مجموع الانحرافات التربيعية في هبوط التدرج العشوائي"

بعد النظر في الجدول الزمني ، كل شيء يقع في مكانه ، والآن سنصلح كل شيء.

اذا ماذا حصل؟ حدث ما يلي. عندما نختار شهرًا بشكل عشوائي ، فهذا هو الشهر المحدد الذي تسعى خوارزميتنا إلى تقليل الخطأ في حساب الإيرادات. ثم نختار شهرًا آخر ونكرر الحساب ، لكننا نقلل الخطأ للشهر الثاني المحدد. والآن لنتذكر أن أول شهرين ينحرفان بشكل ملحوظ عن خط معادلة الانحدار الخطي البسيطة. هذا يعني أنه عند اختيار أي من هذين الشهرين ، عن طريق تقليل خطأ كل منهما ، تزيد الخوارزمية بشكل خطير من الخطأ في العينة بأكملها. اذا مالعمل؟ الجواب بسيط: أنت بحاجة إلى تقليل درجة النزول. في الواقع ، من خلال تقليل خطوة الهبوط ، سيتوقف الخطأ أيضًا عن "القفز" إما لأعلى أو لأسفل. أو بالأحرى ، فإن خطأ "القفز" لن يتوقف ، لكنه لن يفعل ذلك بهذه السرعة :) دعنا نتحقق.

كود لتشغيل SGD بخطوة أصغر

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

الرسم البياني رقم 6 "مجموع الانحرافات التربيعية لنسب التدرج العشوائي (80 ألف خطوة)"

لقد تحسنت الاحتمالات ، لكنها ما زالت غير مثالية. افتراضيًا ، يمكن تصحيح هذا بهذه الطريقة. نختار ، على سبيل المثال ، في آخر 1000 تكرار قيم المعاملات التي حدث بها الحد الأدنى من الخطأ. صحيح ، لهذا سيتعين علينا كتابة قيم المعاملات نفسها. لن نفعل هذا ، بل ننتبه إلى الجدول الزمني. يبدو سلسًا ويبدو أن الخطأ ينخفض بشكل متساوٍ. في الواقع ليس كذلك. دعنا نلقي نظرة على أول 1000 تكرار ونقارنها مع التكرارات الأخيرة.

رمز مخطط SGD (أول 1000 خطوة)

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

الرسم البياني رقم 7 "مجموع الانحرافات التربيعية للدولار السنغافوري (أول 1000 خطوة)"

الرسم البياني رقم 8 "مجموع الانحرافات التربيعية للدولار السنغافوري (آخر 1000 خطوة)"

في بداية الهبوط ، نلاحظ انخفاضًا حادًا ومنتظمًا في الخطأ. في التكرارات الأخيرة ، نرى أن الخطأ يدور حول قيمة 1,475 وفي بعض اللحظات تساوي هذه القيمة المثلى ، لكنها لا تزال ترتفع ... أكرر ، يمكنك تدوين قيم معاملات и ، ثم اختر تلك التي يكون الخطأ فيها ضئيلًا. ومع ذلك ، كانت لدينا مشكلة أكبر: كان علينا أن نتخذ 80 ألف خطوة (انظر الكود) لنقترب من القيم المثلى. وهذا يتناقض بالفعل مع فكرة توفير وقت الحساب في نزول التدرج العشوائي فيما يتعلق بالنسب المتدرج. ما الذي يمكن تصحيحه وتحسينه؟ ليس من الصعب أن نرى أننا في التكرارات الأولى نسير بشكل مطرد ، وبالتالي ، يجب أن نترك خطوة كبيرة في التكرارات الأولى وننقص الخطوة بينما نمضي قدمًا. لن نفعل هذا في هذه المقالة - لقد استمر بالفعل. أولئك الذين يرغبون يمكنهم أن يفكروا بأنفسهم كيف يفعلون ذلك ، فهذا ليس بالأمر الصعب 🙂

سنقوم الآن بتنفيذ نزول التدرج العشوائي باستخدام المكتبة نمباي (ودعنا لا نتعثر فوق الصخور التي حددناها سابقًا)

رمز النسب المتدرج العشوائي (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

تبين أن القيم هي نفسها تقريبًا عند التنازلي دون استخدام نمباي. ومع ذلك ، هذا منطقي.

دعنا نكتشف مقدار الوقت الذي استغرقه هبوط التدرج العشوائي.

رمز لتحديد وقت حساب SGD (80 ألف خطوة)

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)

كلما توغلنا في الغابة ، كانت الغيوم أغمق: مرة أخرى ، تظهر الصيغة "المكتوبة ذاتيًا" أفضل نتيجة. كل هذا يشير إلى أنه يجب أن تكون هناك طرق أكثر دقة لاستخدام المكتبة. نمباي، والتي تعمل على تسريع العمليات الحسابية حقًا. في هذه المقالة ، لن نتعرف عليهم. شيء لتفكر فيه في وقت فراغك :)

نحن نلخص

قبل التلخيص ، أود أن أجيب على سؤال قد نشأ على الأرجح من عزيزنا القارئ. لماذا ، في الواقع ، مثل هذه "العذاب" مع النسل ، لماذا نحتاج إلى السير صعودًا وهبوطًا في الجبل (غالبًا لأسفل) للعثور على الأرض المنخفضة العزيزة ، إذا كان لدينا مثل هذا الجهاز القوي والبسيط في أيدينا ، في شكل حل تحليلي ينقلنا فوريًا إلى المكان الصحيح؟

الجواب على هذا السؤال يكمن في السطح. الآن قمنا بتحليل مثال بسيط للغاية حيث الإجابة الصحيحة هي يعتمد على ميزة واحدة . أنت لا ترى هذا كثيرًا في الحياة ، لذا فلنتخيل أن لدينا علامتين أو 2 أو 30 أو أكثر. أضف إلى هذا الآلاف وحتى عشرات الآلاف من القيم لكل معلم. في هذه الحالة ، قد يفشل الحل التحليلي في الاختبار ويفشل. في المقابل ، فإن الانحدار المتدرج وتنوعاته ستقربنا ببطء ولكن بثبات من الهدف - الحد الأدنى من الوظيفة. ولا تقلق بشأن السرعة - ربما لا نزال نحلل الطرق التي ستسمح لنا بتعيين وتعديل طول الخطوة (أي السرعة).

والآن للحصول على ملخص قصير.

أولاً ، آمل أن تساعد المادة المقدمة في المقالة "علماء البيانات" المبتدئين في فهم كيفية حل معادلات الانحدار الخطي البسيطة (وليس فقط).

ثانيًا ، درسنا عدة طرق لحل المعادلة. الآن ، بناءً على الموقف ، يمكننا اختيار الوضع الأنسب للمهمة قيد البحث.

ثالثًا ، رأينا قوة الإعدادات الإضافية ، وهي طول خطوة نزول التدرج. يجب عدم إهمال هذه المعلمة. كما هو مذكور أعلاه ، لتقليل تكلفة الحوسبة ، يجب تغيير طول الخطوة أثناء الهبوط.

رابعًا ، في حالتنا ، أظهرت الدوال "المكتوبة ذاتيًا" أفضل نتيجة زمنية للحسابات. ربما يرجع ذلك إلى عدم الاستخدام الأكثر احترافًا لقدرات المكتبة. نمباي. ولكن مهما كان الأمر ، فإن الاستنتاج يقترح نفسه على النحو التالي. من ناحية ، يجدر أحيانًا التشكيك في الآراء الراسخة ، ومن ناحية أخرى ، لا يستحق الأمر دائمًا تعقيد كل شيء - بل على العكس ، في بعض الأحيان تكون الطريقة الأبسط لحل المشكلة أكثر فعالية. وبما أن هدفنا كان تحليل ثلاث طرق لحل معادلة الانحدار الخطي البسيطة ، فإن استخدام الدوال "المكتوبة ذاتيًا" كان كافياً بالنسبة لنا.