פתרון משוואת רגרסיה לינארית פשוטה

המאמר דן במספר דרכים לקביעת המשוואה המתמטית של קו רגרסיה פשוט (מזוג).

כל השיטות לפתרון המשוואה הנדונה כאן מבוססות על שיטת הריבועים הקטנים ביותר. נסמן את השיטות באופן הבא:

• פתרון אנליטי
• ירידה בשיפוע
• ירידה בשיפוע סטוכסטי

עבור כל שיטה לפתרון משוואת קו ישר, המאמר מספק פונקציות שונות, המחולקות בעיקר לאלו שנכתבות ללא שימוש בספריה רדום ואלה שמשתמשים לחישובים רדום. הוא האמין כי שימוש מיומן רדום יקטין את עלויות המחשוב.

כל הקוד שניתן במאמר כתוב בשפה פיתון 2.7 עם מחברת צדק. קוד המקור והקובץ עם נתונים לדוגמה מפורסמים ב Github

המאמר מכוון יותר הן למתחילים והן לאלו שכבר החלו להשתלט בהדרגה בלימוד קטע רחב מאוד בבינה מלאכותית - למידת מכונה.

כדי להמחיש את החומר, אנו משתמשים בדוגמה פשוטה מאוד.

תנאים לדוגמה

יש לנו חמישה ערכים המאפיינים תלות Y מ X (טבלה מס' 1):

טבלה מס' 1 "תנאים לדוגמה"

נניח שהערכים הוא חודש השנה, ו - הכנסות החודש. במילים אחרות, ההכנסה תלויה בחודש בשנה, וכן - הסימן היחיד שבו תלויות ההכנסות.

הדוגמה היא כך וכך, הן מנקודת המבט של התלות המותנית של הפדיון בחודש השנה, והן מבחינת מספר הערכים - יש מעט מאוד מהם. עם זאת, פישוט כזה יאפשר, כמו שאומרים, להסביר, לא תמיד בקלות, את החומר שמטמיעים מתחילים. וגם פשטות המספרים תאפשר למי שרוצה לפתור את הדוגמה על הנייר ללא עלויות עבודה משמעותיות.

הבה נניח שניתן להעריך היטב את התלות הניתנת בדוגמה על ידי המשוואה המתמטית של קו רגרסיה פשוט (מזווג) של הצורה:

איפה הוא החודש שבו התקבלה ההכנסה, - הכנסה התואמת לחודש, и הם מקדמי הרגרסיה של הקו המשוער.

שימו לב שהמקדם נקרא לעתים קרובות השיפוע או השיפוע של הקו המשוער; מייצג את הסכום שבו כאשר זה משתנה .

ברור שהמשימה שלנו בדוגמה היא לבחור מקדמים כאלה במשוואה и , שבהן הסטיות של ערכי ההכנסה המחושבים שלנו לפי חודש מהתשובות האמיתיות, כלומר. הערכים המוצגים במדגם יהיו מינימליים.

השיטה הכי פחות ריבועית

לפי שיטת הריבועים הקטנים, יש לחשב את הסטייה בריבוע. טכניקה זו מאפשרת לך להימנע מביטול הדדי של סטיות אם יש להן סימנים הפוכים. לדוגמה, אם במקרה אחד, הסטייה היא +5 (פלוס חמש), ובשני -5 (מינוס חמש), אז סכום הסטיות יבטל זו את זו ויסתכם ב-0 (אפס). אתה לא צריך לריבוע את הסטייה, אלא להשתמש במאפיין המודולוס ואז כל הסטיות יהיו חיוביות ויצטברו. לא נתעכב על נקודה זו בהרחבה, אלא נציין רק שלצורך נוחות החישובים נהוג לריבוע את הסטייה.

כך נראית הנוסחה שבעזרתה נקבע את הסכום הקטן ביותר של סטיות בריבוע (שגיאות):

איפה הוא פונקציה של קירוב של תשובות אמיתיות (כלומר, ההכנסה שחישבנו),

הן התשובות האמיתיות (ההכנסות מסופקות במדגם),

הוא מדד המדגם (מספר החודש שבו נקבעת החריגה)

בואו נבדיל את הפונקציה, נגדיר את משוואות הדיפרנציאליות החלקיות, ונהיה מוכנים לעבור לפתרון האנליטי. אבל ראשית, בואו נצא לסיור קצר על מהי בידול ונזכור את המשמעות הגיאומטרית של הנגזרת.

בידול

דיפרנציאציה היא פעולת מציאת הנגזרת של פונקציה.

למה משמשת הנגזרת? הנגזרת של פונקציה מאפיינת את קצב השינוי של הפונקציה ואומרת לנו את הכיוון שלה. אם הנגזרת בנקודה נתונה חיובית, הפונקציה גדלה; אחרת, הפונקציה יורדת. וככל שהערך של הנגזרת המוחלטת גדול יותר, כך קצב השינוי של ערכי הפונקציה גבוה יותר, וכך גם השיפוע של גרף הפונקציות תלול יותר.

לדוגמה, בתנאים של מערכת קואורדינטות קרטזית, ערך הנגזרת בנקודה M(0,0) שווה ל +25 פירושו שבנקודה נתונה, כאשר הערך מוזז ימינה לפי יחידה קונבנציונלית, ערך גדל ב-25 יחידות קונבנציונליות. בגרף זה נראה כמו עלייה די תלולה בערכים מנקודה נתונה.

דוגמה אחרת. הערך הנגזר שווה -0,1 אומר שכאשר נעקור ליחידה קונבנציונלית אחת, ערך יורד ב-0,1 יחידה קונבנציונלית בלבד. יחד עם זאת, בגרף של הפונקציה, אנו יכולים לראות שיפוע בקושי מורגש כלפי מטה. אם מציירים אנלוגיה להר, זה כאילו אנחנו יורדים לאט מאוד במדרון מתון מההר, שלא כמו הדוגמה הקודמת, שבה היינו צריכים לטפס על פסגות תלולות מאוד :)

כך, לאחר בידול הפונקציה בסיכויים и , אנו מגדירים משוואות דיפרנציאליות חלקיות מסדר 1. לאחר קביעת המשוואות, נקבל מערכת של שתי משוואות, שעל ידי פתרונה נוכל לבחור ערכים כאלה של המקדמים и , שעבורו הערכים של הנגזרות המתאימות בנקודות נתונות משתנים בכמות מאוד מאוד קטנה, ובמקרה של פתרון אנליטי לא משתנים כלל. במילים אחרות, פונקציית השגיאה במקדמים שנמצאו תגיע למינימום, שכן ערכי הנגזרות החלקיות בנקודות אלו יהיו שווים לאפס.

אז לפי כללי ההבחנה, המשוואה הנגזרת החלקית של הסדר הראשון ביחס למקדם יקבל את הטופס:

משוואת נגזרת חלקית מסדר 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

נשווה את הזמן המושקע בקביעת המקדמים и , בהתאם ל-3 השיטות שהוצגו.

קוד לחישוב זמן חישוב

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)

עם כמות קטנה של נתונים, יוצאת קדימה פונקציה "נכתבת בעצמה", המוצאת את המקדמים בשיטת Cramer.

עכשיו אתה יכול לעבור לדרכים אחרות למצוא מקדמים и .

ירידה בשיפוע

ראשית, הבה נגדיר מהו שיפוע. במילים פשוטות, השיפוע הוא קטע המציין את כיוון הצמיחה המקסימלית של פונקציה. באנלוגיה לטיפוס על הר, כאשר פני השיפוע הם המקום בו נמצא הטיפוס התלול ביותר לראש ההר. מפתחים את הדוגמה עם ההר, אנו זוכרים שלמעשה אנו זקוקים לירידה התלולה ביותר על מנת להגיע לשפלה כמה שיותר מהר, כלומר המינימום - המקום בו הפונקציה אינה עולה או יורדת. בשלב זה הנגזרת תהיה שווה לאפס. לכן, אנחנו לא צריכים שיפוע, אלא אנטי שיפוע. כדי למצוא את אנטי-הדרגה אתה רק צריך להכפיל את השיפוע ב -1 (פחות אחד).

הבה נשים לב לעובדה שלפונקציה יכולות להיות כמה מינימות, ולאחר שירדנו לאחת מהן באמצעות האלגוריתם המוצע להלן, לא נוכל למצוא מינימום נוסף, שעשוי להיות נמוך מזה שנמצא. בואו נרגע, זה לא איום עלינו! במקרה שלנו עסקינן במינימום בודד, שכן תפקידנו על הגרף פרבולה רגילה. וכפי שכולנו צריכים לדעת היטב מהקורס שלנו במתמטיקה בבית הספר, לפרבולה יש רק מינימום אחד.

אחרי שגילינו למה אנחנו צריכים שיפוע, וגם שהשיפוע הוא קטע, כלומר וקטור עם קואורדינטות נתונות, שהם בדיוק אותם מקדמים и אנחנו יכולים ליישם ירידה בשיפוע.

לפני שמתחיל, אני מציע לקרוא רק כמה משפטים על אלגוריתם הירידה:

• אנו קובעים בצורה פסאודו אקראית את הקואורדינטות של המקדמים и . בדוגמה שלנו, נגדיר מקדמים ליד אפס. זהו נוהג נפוץ, אך לכל מקרה עשוי להיות נוהג משלו.
• מתוך קואורדינטה להפחית את הערך של הנגזרת החלקית מסדר 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 "סכום הסטיות בריבוע SGD (1000 הצעדים הראשונים)"

גרף מס' 8 "סכום הסטיות בריבוע SGD (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, 50 או יותר סימנים. בואו נוסיף לזה אלפי, או אפילו עשרות אלפי ערכים עבור כל תכונה. במקרה זה, ייתכן שהפתרון האנליטי לא יעמוד במבחן ויכשל. בתורו, ירידה בשיפוע והווריאציות שלה יקרבו אותנו לאט אבל בטוח אל המטרה - המינימום של הפונקציה. ואל תדאג לגבי המהירות - כנראה שנבדוק דרכים שיאפשרו לנו להגדיר ולווסת את אורך הצעד (כלומר, מהירות).

ועכשיו הסיכום הקצר ממש.

ראשית, אני מקווה שהחומר המוצג במאמר יעזור ל"מדעני נתונים" מתחילים בהבנה כיצד לפתור משוואות רגרסיה ליניאריות פשוטות (ולא רק).

שנית, בדקנו כמה דרכים לפתור את המשוואה. כעת, בהתאם למצב, נוכל לבחור את המתאים ביותר לפתרון הבעיה.

שלישית, ראינו את הכוח של הגדרות נוספות, כלומר אורך צעד הירידה בשיפוע. אי אפשר להזניח את הפרמטר הזה. כפי שצוין לעיל, על מנת להוזיל את עלות החישובים, יש לשנות את אורך הצעד במהלך הירידה.

רביעית, במקרה שלנו, פונקציות "כתובות בבית" הראו את תוצאות הזמן הטובות ביותר לחישובים. זה כנראה נובע מהשימוש הלא מקצועי ביותר ביכולות הספרייה רדום. אבל כך או כך, המסקנה הבאה מעידה על עצמה. מצד אחד, לפעמים כדאי להטיל ספק בדעות מבוססות, ומצד שני, לא תמיד כדאי לסבך הכל - להיפך, לפעמים דרך פשוטה יותר לפתרון בעיה יעילה יותר. ומכיוון שהמטרה שלנו הייתה לנתח שלוש גישות לפתרון משוואת רגרסיה לינארית פשוטה, השימוש בפונקציות "נכתבות בעצמן" הספיק לנו בהחלט.

ספרות (או משהו כזה)

1. רגרסיה לינארית

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

2. שיטת הריבועים הקטנים ביותר

mathprofi.ru/metod_naimenshih_kvadratov.html

3. נגזרת

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. מעבר צבע

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. ירידה בשיפוע

habr.com/en/post/471458

habr.com/en/post/307312

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

6. ספריית 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

מקור: www.habr.com