رگراسیون لجستیک در پایتون!

سلام! (= توی این مطلب میخوایم رگراسیون لجستیک (یا به قول بعضیا لاااجستیک!) رو به صورت وکتورایزد (vectorized) یا برداری شده پیاده کنیم! و البته! بدون استفاده از لایبرری sklearn! همش با خودمون! منظورم از وکتورایزد هم اینه که یه کاری کنیم که به جای اجرا شدن با لوپ و این نوع محاسبات وقتگیر، از محاسبات برداری و ماتریسی جبر خطی استفاده بشه که باعث میشه سرعت اجرا فوق العاده بالا بره. قبل از شروع، اگه با رگراسیون لجستیک به صورت تئوری آشنا نیستین، حتما این مطلبم درباره همین موضوع رو بخونین! خب! شروع کنیم!

واسه ی شروع، از کتابخونه های نامپای (Numpy)، پانداز (Pandas) و مت پلات لیب (Matplotlib) استفاده میکنیم که به ترتیب واسه ی کار با ماتریس ها و بردارها، کار با دیتاست و نمایش نمودار استفاده میشن. این سه تا کتابخونه رو ایمپورت میکنیم:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

این لینک دیتافریمی هست که قراره ازش استفاده کنیم! که تشکیل شده از یه تعداد عدد به عنوان x و یه تعداد 0 و 1 به عنوان y. ینی این شکلی:

هدف اینه که بتونیم یه تابع سیگموید (sigmoid) رو طوری روی این نمودار جا بدیم که به ازای مقادیر بزرگتر یا مساوی 25 (تقریبا!) مقدار تابع از 0.5 بزرگتر یا مساویش بشه (و در نتیجه جوابش رو 1 درنظر بگیریم)؛ و همینطور برای مقادیر کوچیکتر از اون مرز، مقدار تابع هم از 0.5 کوچیکتر باشه. خب. یه تابع initializer تعریف میکنیم که دیتا رو برامون لود میکنه:

def initializer():
    df = pd.read_csv('dataframe.csv')
    x = df.x.values.reshape(-1,1).T
    y = df.y.values.reshape(-1,1).T
    w = np.random.random((1,x.shape[0]))
    b = 0
    x = np.divide(x, 100)
    print('x shape:',x.shape,', w shape:',w.shape,', y shape:',y.shape)
    return x, y, w, b

خب یه سری نکته هست که درباره ی این چند خط باید بگم. توی خط اول که تابع رو تعریف کردیم، توی خط دوم به کمک کتابخونه پانداز، دیتافریم رو لود کردیم و ریختیم توی یه متغیر؛ و اما دو خط بعدی. از اونجایی که تصمیم گرفتیم به صورت وکتورایزد عمل کنیم، از این روش واسه ی تعریف x و y استفاده کردیم. واسه ی این که دقیقتر متوجه بشین توی اون دو خط چه اتفاقی میفته، به اینجا دقت کنین. اگه فقط از عبارت x = df.x.values استفاده کنیم، یه آرایه ی به ابعاد یک در m مثل [1,2,3] خروجی میده که m به تعداد ایکس های موجود در دیتافریمه. حالا اگه به صورت x = df.x.values.reshape(-1,1) بنویسیمش، میاد هر عنصر رو توی یه سطر جا میده و ابعادش میشه m در یک. ینی مثلا:

و اگه یه T. به آخرش اضافه کنیم این بردار ترانهاده (Transpose) میشه و به فرم [[3] ,[2] ,[1]] درمیاد. حالا ممکنه بپرسین چرا این همه سختی به خرج دادیم و به این فرم درآوردیمش؟! مگه [1,2,3] چش بود؟! (=

نکته توی خط 5 نهفته شده. ینی w (که همون θ₁ عه). به این صورت تعریفش کردیم که اولا رندوم باشه، دوما ابعادش یک در تعداد ستون های x باشه. حالا چرا؟ درواقع میتونستیم همون اول به جای عبارت x.shape[0] بنویسیم یک؛ اما صورتی از مسئله رو درنظر بگیرین که هر ایکسمون به جای یه عدد، چندین عدد باشه. مثلا طبق مثال معروف فروش بستنی بر حسب دما، عوامل دیگه ای هم دخیل باشن، مثل طعم یا رنگ بستنی( اگه ذغالی باشه تو دمای 273- هم میشه خورد??). در اون حالت دیگه یه دونه w برای کل ایکس ها کافی نیست و باید تعدادش به تعداد فیچر (feature) های ایکس باشه.

(مثال پیچیده(!): اگه ایکس 2 در 27 باشه ینی دوتا فیچر و 27 تا نمونه، پس w میشه 1 در 2. ینی دوتا دونه w! دقت کنین که ضرب ماتریس x در بردار w باید قابل انجام باشه!)

در خط بعدی متغیر b (که همون θ₀ عه) رو مساوی صفر گذاشتیم، و یه خط مرموز! دلیل اینکه چرا تمام اعداد ایکس رو تقسیم بر 100 کردم رو بعدا بهتون میگم؛ فعلا بپذیرین (= و درنهایت هم چاپ ابعاد x و y و w و خروجی دادنشون به همراه b.

x, y, w, b = initializer()
plt.scatter(x, y, color='red', s=60)
plt.xkcd()
plt.show()

توی این 4 خط، اول تابع initializer رو صدا زدیم و خروجی هاشو ریختیم توی 4 تا متغیر، و درنهایت نمودار x و y رو به کمک کتابخونه مت پلات لیب چاپ کردیم. (اون خط سوم باعث میشه نمودار کیوت نشون داده بشه! (= ) و خروجی به این صورت درمیاد:

خب. حالا تابع فرض (hypothesis) که همون تابع سیگموید هست رو تعریف میکنیم:

def sigmoid(x):
    return 1/(1 + np.exp(-x))

یادتونه اون بالا مقادیر ایکس رو بر 100 تقسیم کردم؟ جوابش اینجاست! مشکلی که در صورت تقسیم نکردن پیش میاد، اینه که بخش np.exp(-x) فوق العاده کوچیک میشه (ینی مثلا e به توان منفی 100!!!) و این مقدار خیلی کوچیک به عنوان صفر تعبیر میشه و عملا خروجی تابع سیگموید مساوی یک میشه. این یک توی تابع هزینه (که عبارتی مثل log(1 - hθ(x)) داره) باعث میشه به log(0) بخوریم و به سوی بی نهایت و فراتر از آن! پس برای رفع این مشکل سعی میکنیم ایکس ها رو کوچیک کنیم و مثل اینجا تقسیم بر یه عدد کنیم؛ و یا روش دیگه ای به نام فیچر اسکیلینگ (Feature Scaling) که اصولی تره. خب. ادامه میدیم:

def cost(y_h, y):
    m = y_h.shape[0]
    cost = 0
    cost = ((-y * np.log(y_h)) - (1 - y) * np.log(1 - y_h)) / m
    return cost.sum()

اینجا تابع هزینه رگراسیون لجستیک رو به صورت کد درآوردیم. اگه یادتون باشه تابعش این شکلی بود:

اگه هم احیانا فرم شرطی این فرمول رو دیده بودین ایرادی نداره. این همون تابعه؛ فقط اگه y مساوی صفر باشه، عبارت اولی حذف میشه و اگه مساوی یک باشه، عبارت دومی حذف میشه. توی کد بالا هم دوتا متغیر به عنوان ورودی پاس داده شدن که اولی y_h عه و همون حدسیه که سیستم زده (hθ(x)) و دومی y، که جواب اصلی اون نمونه هست. m تعداد نمونه ها رو میگیره و خط بعدی هم مقدار اولیه کاست رو مساوی صفر میذاره. تنها نکته ای که وجود داره اینه که کاست مساوی یه بردار میشه (چون y و y_h هم بردارن!) و ما درنهایت جمع عناصر این بردار رو خروجی میدیم. فقط یه عدد!

حالا باید کدی رو بزنیم که قراره این تابع هزینه رو حداقل کنه. اینجا من کد گرادیان نزولی (Gradient Descent) رو زدم و اگه بخوام کامل توضیحش بدم، بهتره یه مطلب دیگه واسش بنویسم?. میتونین این بخش رو ران کنین و بگذرین ازش!

def GD(w, b, y, iteration, alpha, limit):
    list_iter = []
    list_cost = []
    for i in range(iteration):
        w = w - alpha * np.sum((sigmoid(np.dot(w,x) + b) - y) * x)
        b = b - alpha * np.sum(sigmoid(np.dot(w,x) + b) - y)
        if i % limit == 0:
            list_iter.append(i)
            list_cost.append(cost(sigmoid(w.dot(x) + b), y))
    return w, b, list_cost, list_iter

به طور کلی کاری که این بخش کد انجام میده اینه که مقادیر وزن مثل w و b رو ورودی میگیره و طی یه حلقه، اونا رو یه طوری کم و زیاد میکنه که کاست حداقل بشه. همچنین دوتا لیست خالی هم در ابتداش تعریف کردم که به ازای هرچند بار اجرا شدن حلقه، مقدار تابع کاست و همون تعداد اجرا شدن حلقه رو داخلشون ذخیره کنه که بعدا بتونیم به تصویر بکشیمش. به کمک این کار میتونیم متوجه بشیم که آیا الگوریتم حداقل کردن کاست به درستی کار میکنه یا نه.

حالا از تمام کد های قبلی به این صورت استفاده میکنیم:

print('w and b before Gradient Decsent: ', w, b)
w, b, list_cost, list_iter = GD(w, b, y, 1000, 0.1, 10)
print('w and b after Gradient Decsent: ', w, b)
plt.plot(list_iter,list_cost)
plt.show()

که این خروجی رو میده:

توی این بخش از کد، اولا تابع گرادیان نزولی رو ران کردم و بعدش نمودار اون دوتا لیست کاست ها و شمارش ها رو چاپ کردم.همونطور که مشخصه، کاست بعد از تعدادی شمارش به حداقل رسیده و این ینی الگوریتم به درستی کار میکنه. و از اون مهمتر! ما w و b برای کشیدن اون تابع سیگمویدی که میخواستیم رو داریم! حالا نمودارشو رسم میکنیم:

y_test = sigmoid(w.dot(x) + b)
plt.scatter(x, y_test, s=10, color='b')
plt.scatter(x, y, color='r')
plt.show()

و خروجی رویایی:

کاری که توی چند خط آخر کردم این بود که به ازای نمونه های ایکس، تابع سیگموید رو (به همراه w و b جدید!) فراخوانی کردم و ریختمشون توی y_test. حالا یه بار نمودار x و y_test رو کشیدم و یه بار x و y.

حالا تنها کاری که لازمه اینه که یه شرط بذاریم برای اینکه آیا جواب از 0.5 کوچیکتره یا نه. اگه آره پس خروجی مساوی صفره، وگرنه مساوی یک. در هر صورت، بخش سخت کارو گذروندیم و بهتره مطلب رو از اینی که هست طولانی تر نکنم!?

تبریک میگم! شما اولین کد رگراسیون لجستیکتون رو زدین!?