Logistic Regression Model

โรงงานผลิตไมโครชิพแห่งหนึ่ง เมื่อผลิตไมโครชิพเสร็จแล้ว ก็จะมีกระบวนการตรวจสอบมาตฐานไมโครชิพอยู่ด้วยกัน 6 แบบทดสอบ ก่อนที่จะถูกรับรองมาตรฐานและแพคเข้าบรรจุภัณฑ์

ในแต่ละด่านการทดสอบ ก็จะมีค่าใช้จ่ายในกระบวนการ. สำหรับด่านแรก คิดเป็นเงิน 5% ของค่าใช้จ่ายทั้งหมด ด่านที่สองอีก 5% ดังภาพข้างต้น. ชิพบางอันผ่านการทดสอบทั้งห้าด่านแรก แต่ไปตกด่านที่หก ก็ต้องถูกนำมาแยกชิ้นส่วนเพื่อเข้าสู่จุดเริ่มต้นการผลิตใหม่ จนกว่าจะผ่านการทดสอบทั้งหกด่าน ถึงจะสามารถนำออกสู่ตลาดได้ 

คราวนี้บริษัทต้องการลดต้นทุนการผลิต และต้องการประหยัดเวลาในการทดสอบ โดยอยากทำนายว่าชิพแต่ละอันจะผ่านทั้งหกด่านหรือไม่ ด้วยการพิจารณาข้อมูลที่วัดได้จาการทดสอบที่หนึ่ง กับสอง (เพราะมีค่าใช้จ่ายน้อย และไม่เสียเวลาในการทดสอบมาก) จากข้อมูลสองอย่างนี้ เขาอยากทำนายว่าควรจะส่งเข้าแบบทดสอบต่อไปหรือเปล่า หรือควรจะนำกลับไปผลิตใหม่ทันที และไม่ต้องเสียค่าใช้จ่ายอีก 90% ที่เหลือ.

ตัวอย่างค่าที่วัดได้จากด่านที่หนึ่ง และด่านที่สอง กับผลลัพท์สุดท้ายเมื่อทดสอบไปหมดทั้งหกด่าน

Test 1	Test 2	Final Result
0.051267	0.69956	Accepted
-0.092742	0.68494	Accepted
-0.06394	-0.18494	Accepted
0.18376	0.93348	Rejected
-0.13306	-0.4481	Rejected
-0.10426	0.99196	Rejected

ตัวอย่างข้อมูลที่บริษัทเก็บมาให้ดังนี้ โดยที่ 1 คือผ่าน และ 0 คือไม่ผ่าน

คำถามถัดมาคือ แล้วเราจะเขียนโปรแกรมเพื่อทำนายผลลัพท์จากข้อมูลการทดสอบเพียงสองด่านนี้ได้อย่างไร ?

เริ่มจากวิธีที่ง่ายที่สุดคือการ เอาข้อมูลทั้งหมดที่มีมาพล็อตกราฟ เพื่อหาแพทเทิร์น ดังนี้

	"""Logistic Regression Model"""
	import matplotlib.pyplot as plt
	import numpy as np
	def plot_data(data, window_title):
	"""
	:param data: np array of shape (n,3)
	:param window_title: String
	"""
	fig, ax = plt.subplots()
	fig.canvas.set_window_title(window_title)
	accepted = data[data[:, 2] == 1, :]
	rejected = data[data[:, 2] == 0, :]
	ax.scatter(accepted[:, 0], accepted[:, 1], c='blue', label='Accepted')
	ax.scatter(rejected[:, 0], rejected[:, 1], c='red', label='Rejected')
	ax.set_xlabel("Microchip Test 1")
	ax.set_ylabel("Microchip Test 2")
	ax.legend()
	ax.grid(True)
	plt.show(block=False)
	def run():
	"""
	Entry point
	"""
	data = np.loadtxt("./data/ex2data2.txt", delimiter=",")
	plot_data(data, "Test 1 & 2 data")
	if __name__ == '__main__':
	run()

view raw ThoughtLog_LRM01.py hosted with ❤ by GitHub

ก่อนรันโค้ดนี้ สร้างโฟเดอร์ data ขึ้นมาแล้ววางไฟล์ข้อมูล ex2data2.txt

จากกราฟ เราบอกได้คร่าวๆว่า ชิพที่มีโอกาสผ่านการทดสอบทั้งหกด่านไปได้ จะมีผลลัพท์ของเทสที่หนึ่งอยู่ระหว่าง -0.5 กับ 0.75 และ ผลลัพท์ของเทสที่สองอยู่ระหว่าค่า -0.5 กับ 0.75 ซึ่งการเขียนโปรแกรมแบบนี้ก็ไม่ได้ยาก ใช้ if - else ไม่กี่ตัวก็สามรถประเมินได้ง่ายๆ แต่การประเมินคร่าวๆนี้ยังมีความผิดพลาดค่อนข้างมาก เพื่อเพ่ิ่มโอกาสผ่านเข้ารอบสำหรับจุดสีน้ำเงิน และลดโอกาสจุดสีแดง เราสามารถวาดเส้นที่ใกล้เคียงได้มากว่าสี่เหลี่ยมเช่นรูปภาพด้านล่าง

คำถามถัดมาคือ เส้นที่วาดลงบนกราฟนั้นเขียนออกมาเป็นสมการได้อย่างไร? ต้องไม่ลืมว่าไม่ใช่สมการ $y=something...x$ แต่เป็น $y=something...x_1 , x_2, ... x_n$(ถ้าใครสังเกตดีๆ ก็จะเห็นเส้นสีดำเขียน 0.000 กำกับไว้ ไว้อ่านโพสนี้จบ ก็จะเข้าใจที่มาของมัน). คราวนี้สมมุติว่าเราได้สมการนั้นมาชื่อว่า $h_{\theta}(x_1,x_2)$ และเราก็ต้องการให้ผลลัพท์อยู่ระหว่าง 0 กับ 1 คือ Accepted หรือ Rejected นั่นเอง

$0 \le h_{\theta}(x_1,x_2) \le 1$

พอป้อนค่าของ Test1 และ test2 เข้าไป เราก็จะได้ผลลัพท์เพื่อที่จะเอาไปทำนายได้เช่น ถ้า $h_{\theta}(x_1,x_2) \le 0.5$ ก็แสดงว่าไม่ควรจะเสียเวลาไปเทสต่อ. คำถามถัดมาอีกคือ $h_{\theta}(x_1,x_2)$ จะต้องมีหน้าตายังไง?

$h_{\theta}(x1,x2)$ สามารถกระจายรูปออกเป็นดังนี้

$h_{\theta}(x1,x2) = \theta_0 + \theta_{1}x_1 + \theta_{2}x_2$

ถ้าให้ $x_1$ แสดงค่าของ Test1 และ $x_2$ แสดงค่าของ Test2 ก็สามารถเขียนออกมาเป็นสมการ

$\begin{aligned}\theta_0 + \theta_{1}x^{(0)}_1 + \theta_{2}x^{(0)}_2 &=1 \\\theta_0 + \theta_{1}x^{(1)}_1 + \theta_{2}x^{(1)}_2 &=1 \\... &=...\\... &=...\\\theta_0 + \theta_{1}x^{(117)}_1 + \theta_{2}x^{(117)}_2 &=0 \end{aligned}$

หรือเขียนในรูปผลคูณของเมทริกซ์

$\begin{bmatrix}1 & x^{0}_1 & x^{0}_2 \\ 1 & x^{1}_1 & x^{2}_2 \\&...&\\1 & x^{117}_1 & x^{117}_2 \\\end{bmatrix}\begin{bmatrix}\theta_0 \\ \theta_1 \\\theta_2 \end{bmatrix}=\begin{bmatrix}y_0 \\ y_1 \\...\\y_{117}\\\end{bmatrix}$

ณ จุดนี้ ปัญหาหลักๆสามอย่างคือ

1. จะหาค่า $\theta_0, \theta_1$ และ $\theta_2$ ได้ยังไง
2. เป็นไปได้ขนาดไหนที่ว่า ลำพัง $\theta$ แค่สามตัวจะช่วยให้วาดกราฟที่เหมาะสมกับรูปแบบของข้อมูลเหล่านั้นได้
3. ถึงแม้จะได้ $\theta$ ออกมาจริงๆ แต่ $y$ อาจมีค่ามากกว่า 1 หรือ น้อยกว่า 0 ก็ได้ ซึ่งไม่ใช่สิ่งที่เราต้องการ

เพื่อจะตอบคำถามทั้งสามข้อได้ เราจะต้องเข้าใจเรื่อง Feature mapping, Function Optimization และ Logistic Function. ถึงจะเอาปัญหาเหล่านั้นอยู่

Logistic Function

Logistic Function ในที่นี้สามรถใช้อีกฟังก์ชั่นหนึ่งที่แทนกันได้คือ Sigmoid Function โดยรวมแล้วก็คืออันเดียวกัน.

จากกราฟ ไม่ว่า $x$จะมากขนาดไหน $y$ ก็ไม่มีทางเกิน 1. เช่นกันกับที่ไม่ว่า $x$ จะติดลบขนาดไหน $y$ ก็ไม่มีทางต่ำกว่า 0. ด้วยลักษณะพิเศษนี้ เราจะเอาคุณสมบัตินี้มาแก้ปัญหาข้อที่ 3 ของเรา ด้วยการพลิกแพลงสมการนิดหน่อย โดยการเอาค่าที่ได้จาก $h_{\theta}(x)$ ไปยัดใส่ในฟังก์ชั่น Sigmoid ดั่งตามนี้

$\begin{aligned}h_{\theta}(x) &=g(x\theta) \\ z &=x\theta \\ g(z) &=\dfrac{1}{1+e^{-z}}\end{aligned}$

แบบนี้แล้ว $h_{\theta}(x)$ ของเราก็จะมีค่าระหว่าง 0 กับ 1 อย่างที่ต้องการ. $h_{\theta}(x)$ บ่งบอกความน่าจะเป็นของผลลัพท์ที่เท่ากับ 1. ยกตัวอย่างเช่น $h_{\theta}(x) = 0.7$ จะหมายถึงโอกาส 70% ที่ผลลัพท์จะเท่ากับ 1 และ 30% เท่ากับ 0.

$h_{\theta}(x) = P(y=1|x;\theta) = 1 - P(y = 0|x;\theta)$

$P(y=1|x;\theta) + P(y = 0|x;\theta)= 1$

และนี่ก็เป็นคำตอบที่ว่าเส้นสีดำที่เขียนกำกับ 0.000 นั้นก็หมายความว่า ถ้า $z = 0 , g(0) = 0.5$ ดังการฟของ Sigmoid, ก็หมายความว่าผลการทำนายคือ Accepted ( >= 0.5). เส้น 0.000 นี้จึงเป็นเส้นที่แบ่งอาณาเขตระหว่าจุดสีน้ำเงินกับสีแดง

Feature mapping

$x_1 , x_2$ ของเราเรียกว่า Feature. Feature Mapping ก็คือการเอาสองค่านั้นไป map แล้วก็ได้ค่าที่ผ่านการ map ไปใช้ต่อ. แล้วทำไมต้องทำ Feature Mapping ด้วย? เหตุผลหลักๆก็คือ ข้อมูลเราอาจมีไม่เพียงพอที่จะนำไปสร้างกราฟที่มีความซับซ้อนเหมือนเส้น 0.000 นั้นได้ หากเราไม่ทำการแมบเลย ข้อมูลที่ใช้ได้ก็จะมีแค่ $x_1 , x_2$ ทำได้อย่างมาก็แค่สร้างสมการประมาณนี้ $y = \theta_0 + \theta_1x_1 + \theta_2x_2$ ก็จะได้แค่สมการเส้นตรงประมาณนี้

แต่ถ้าเราเพิ่ม order มันให้มากขึ้น ก็จะสามารถสร้างกราฟที่มีความซับซ้อนมากขึ้นไปตามลำดับ อย่างเช่น

สิ่งที่เราต้องการถัดมาคือฟังก์ชั่นที่ใช้ในการ map ข้อมูลสองอย่างให้ออกมามี order ที่มีมากขึ้น

$map\_feature(x_1,x_2,6)=\begin{matrix}[x_{1}, & x_{2}, & x^{2}_{1}, & x_{1}x_{2}, & \\ x^{2}_{2}, & x^{3}_{1}, & x^{2}_{1}x_{2}, & x_{1}x^{2}_{2}, & \\ x^{3}_{2}, & x^{4}_{1}, & x^{3}_{1}x_{2}, & x^{2}_{1}x^{2}_{2}, & \\ x_{1}x^{3}_{2}, & x^{4}_{2}, & x^{5}_{1}, & x^{4}_{1}x_{2}, & \\ x^{3}_{1}x^{2}_{2}, & x^{2}_{1}x^{3}_{2}, & x_{1}x^{4}_{2}, & x^{5}_{2}, & \\ x^{6}_{1}, & x^{5}_{1}x_{2}, & x^{4}_{1}x^{2}_{2}, & x^{3}_{1}x^{3}_{2}, & \\ x^{2}_{1}x^{4}_{2}, & x_{1}x^{5}_{2}, & x^{6}_{2}]\end{matrix}$

	def map_feature(X1, X2, degree=6):
	"""
	Map Feature
	:param X1: NP array of shape (n,)
	:param X2: NP array of shape (n,)
	:param degree: integer denote the highest order
	:return:
	"""
	m = X1.shape[0]
	X1 = X1.T
	X2 = X2.T
	out = np.ones((m, 1), dtype=float)
	for i in range(1, degree + 1):
	for j in range(i + 1):
	value = ((X1 ** (i - j)) * (X2 ** j)).reshape(m, 1)
	out = np.append(out, value, axis=1)
	return out
	def run():
	"""
	Entry point
	"""
	data = np.loadtxt("./data/ex2data2.txt", delimiter=",")
	plot_data(data, "Test 1 & 2 data")
	# Number of samples
	m = data.shape[0]
	# Sample length
	n = data.shape[1]
	X = data[:, 0:n - 1].reshape((m, n - 1))
	X = np.insert(X, 0, 1, axis=1)
	X = map_feature(X[:, 1], X[:, 2])
	print("X.shape", X.shape)

view raw ThoughtLog_LRM02.py hosted with ❤ by GitHub

พอ map features ออกมาได้แล้ว ทีนี้ $\theta$ ของเราจากเดิมมีแค่ 3 ก็จะกลายไปเป็น ตามจำนวนฟีเจอร์ + 1 คือ

$\theta=\begin{bmatrix}\theta_0 \\ \theta_1 \\\vdots \\ \theta_n\end{bmatrix}$

ณ ตอนนี้ เราได้ข้อมูลทุกอย่างที่จะใช้ในการสร้างโมเดลการทำนายแล้ว นั่นคือ

$sigmoid \begin{pmatrix} \begin{bmatrix}1 & x^{0}_1 & x^{0}_2 & \dots & x^6_2 \\ 1 & x^{1}_1 & x^{2}_2 & \dots & x^6_2 \\&...&\\1 & x^{117}_1 & x^{117}_2 & \dots & x^6_2 \\\end{bmatrix}\begin{bmatrix}\theta_0 \\ \theta_1 \\\vdots \\ \theta_{28}\end{bmatrix}\end{pmatrix}=\begin{bmatrix}y_0 \\ y_1 \\...\\y_{117}\\\end{bmatrix}$

ปรับแต่ง $\theta$ ให้เหมาะสม

ณ ตอนนี้ ถ้าเราบอกว่าเซต $\theta$ ทุกตัวของเราให้เป็น $0$ ทั้งหมด สิ่งที่ได้ก็คือการทำนายว่า ทุกไมโครชิพผ่านการทำนายทั้งหมด ($sigmoid(0) = 0.5$) แต่ในความเป็นจริงไม่ใช่ ผิดไปตั้ง 60 กว่าจุด การทำนายผิดไป 60 กว่าจุดนี้เรียกว่า cost ของการเซต $\theta$ ให้เป็น 0 ทั้งหมด.

งานที่เราต้องทำถัดไปคือ หาว่าจะปรับแต่ง $\theta$ ทั้ง 28 ตัวให้แต่ละตัวมีค่าเท่าไหร่ จึงจะเกิด cost น้อยที่สุด! แต่ก่อนจะไปหาวิธีการปรับแต่ง $\theta$ เราต้องมีวีธีการทำคำนวน cost ของ $\theta$ กันก่อน

Cost Function

ตั้งสมการไว้ก่อน แล้วค่อยๆมาอธิบายกัน

$Cost(h_{\theta}(x),y)=\begin{cases}-log(h_{\theta}(x)) &\text{if}\space y=1 \\ -log(1-h_{\theta}(x)) &\text{if}\space y=0\end{cases}$

เพื่อการเห็นภาพที่ชัดเจนและเข้าใจได้ง่ายถึง ลองเอาสองสมการมาพล๊อตเป็นกราฟก็จะได้

คิดง่ายๆว่า เส้นแรก สีม่วง $-log(h_{\theta}(x)) \space \text{if}\space y=1$ เอาไปใช้คำนวณ cost ของ $h_{\theta}(x)$ เฉพาะในกรณีที่ตัวอย่างนี้มีผลลัพท์จริงๆคือ 1 (ที่ได้ระบุมาแล้วจากข้อมูลเก็บ $y = 1$ ) พูดอีกอย่างหนึ่งคือ ถ้าเราทำนายถูกว่าได้ค่าออกมาเป็น 0.85 ดังนั้นแล้ว cost (หรือความผิดพลาด) ก็จะมีค่าออกมาเป็น $-log(0.85) = 0.1625$ ซึ่งเป็น cost ที่ต่ำมากๆ. ในทางตรงกันข้าม ถ้าเราเกิดทำนายผิด กลายเป็นว่าเราทำนายได้ 0.12 ซะงั้น เราก็จะได้ cost ออกมาเป็น $-log(0.12) = 2.12026$ ซึ่งสูงมาก ถ้าเทียบกับ 0.1625. เช่นเดียวกันกับ สมการที่ 2 คือ $-log(1-h_{\theta}(x)) \space \text{if}\space y=0$ ก็มีหลัการคิดเช่นเดียวกัน

ตอนนี้เราได้ Cost Function เป็นที่เรียบร้อยแล้ว ต่อไปเราก็จะมาหาวิธีการลดค่าให้เหลือน้อยที่สุด โดยการลองเปลี่ยนค่า $\theta$ ไปเรื่อยๆ จนกว่าจะได้ Cost Function ที่น้อยที่สุด

เนื่องจากการแยก Cost Function ออกเป็น 2 สมการ ทำให้ยุ่งยากต่อการคำนวณ เราจึงรวมสองสมการนั้นเข้าไว้ด้วยกัน ดังนี้

$Cost(h_{\theta}(x),y)=-ylog(h_{\theta}(x)) -(1-y)log(1-h_{\theta}(x))$

ทำไมจึงเลือก Cost Function ให้เป็นแบบนั้น มีตัวเลือกอื่นนอกจากนี้อีกไหม?

เรื่องนี้เริ่มจะออกนอกสโคปของโพสนี้แล้ว แต่สั้นๆคือได้มาจาก Principle of Maximum likelihood estimation. เป็นศาสตร์ทางวิชาสถิติที่เอาไว้ใช้หาตัวแปรของข้อมูล

กลับมาสู่เรื่องสำคัญอีกครั้งคือ วิธีการลด Cost Function ให้เหลือน้อยที่สุด. ในที่นี้เราจะเอาแคลคูลัสเข้ามาช่วย โดยการดิฟสมการ Cost Function เพื่อเอาผลลัพท์มาปรับค่า $\theta$ นั่นเอง

Gradient Desent

จากสมการ Cost Function

$Cost(h_{\theta}(x),y)=-ylog(h_{\theta}(x)) -(1-y)log(1-h_{\theta}(x))$

$x,y$ จะอยู่ในรูปของ row หนึ่งๆ ($[\theta,x_1,x_2,y]$) เราจะเขียนแยกออกมาเพื่อให้เข้าใจได้ง่ายขึ้น และที่เราทำเป็นเมทริกซ์ก็เพื่อให้ง่ายและรวดเร็วในการให้คอมพิวเตอร์คำนวน ถ้าจะคำนวนหาจำนวน cost ทั้งหมดทุกตัวอย่างเราก็จะได้

$J(\theta)=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}y^{(i)}log(h_{\theta}(x^{(i)})) + (1-y)log(1-h_{\theta}(x^{(i)})) \Big]$

ถ้า $\theta$ ครั้งที่หนึ่งของเราได้ Cost สมมุติว่าอยู่ทีประมาณ 198 ดังนั้น เราต้องการหา $\theta$ ครั้งที่สอง เพื่อจะได้ Cost ที่ต่ำกว่า 198 วิธีหนึ่งที่จะการันตีว่าการหา $\theta$ ครั้งถัดไป ได้ Cost น้อยกว่าของเดิมแน่นอน นั่นก็คือ

$\theta_{j} := \theta_{j} - \alpha\dfrac{ \partial}{\partial \theta_{j}}J(\theta_{j})$

และแน่นอน เราจะไม่ปล่อยให้ $\dfrac{ \partial}{\partial \theta_{j}}J(\theta_{j})$ ลอยนวล ว่าแล้วจะรอช้าอยู่ใย ม่ะจัดเลย

$\begin{aligned}J(\theta) &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[y^{(i)}log(h_{\theta}(x^{(i)})) + (1-y)log(1-h_{\theta}(x^{(i)}))] \Big] \\given \space z &=x\theta\\given \space A &=ylog(\dfrac{1}{1+e^{-z}}) \\given \space B &=(1-y)( log( 1 - \dfrac{1}{1 + e^{-z}}) )\\J(\theta) &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ A + B] \Big]\\ A &=y\bigg[ log(1) - log(1+e^{-z}) \bigg] \\\space &=y\bigg[ 0 - log(1+e^{-z}) \bigg]\\\space &=-ylog(1+e^{-z}) \\B &=(1-y)( log( \dfrac{1 + e^{-z}-1}{1 + e^{-z}}) )\\\space &=(1-y)( log( \dfrac{e^{-z}}{1 + e^{-z}}) ) \\\space &=(1-y)( log( e^{-z}) - log(1 + e^{-z}) ) \\\space &=(1-y)( -z - log(1 + e^{-z}) ) \\\space &=-z - log(1 + e^{-z}) + yz + ylog(1 + e^{-z}) \\A + B &=-ylog(1+e^{-z}) -z - log(1 + e^{-z}) + yz + ylog(1 + e^{-z}) \\\space &=- \cancel{ylog(1+e^{-z})}-z - log(1 + e^{-z}) + yz + \cancel{ylog(1 + e^{-z})}\\\space &=-z - log(1 + e^{-z}) + yz \\\space &=yz - (z + log(1 + e^{-z})) \\\space &=yz - (log(e^{z}) + log(1 + e^{-z})) \\\space &=yz - log(e^{z}(1 + e^{-z})) \\\space &=yz - log(e^{z}+ 1) \\\space &=yx\theta - log(e^{x\theta}+ 1) \\J(\theta) &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ yx\theta - log(e^{x\theta}+ 1)] \Big]\\\dfrac{\partial J(\theta)}{\partial \theta}&=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}\dfrac{\partial [ yx\theta - log(e^{x\theta}+ 1)]}{\partial \theta}\Big]\\\space &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ yx - \dfrac{e^{x\theta}x}{e^{x\theta}+ 1}] \Big]\\\space &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ yx - \dfrac{x}{1 + e^{-x\theta}}] \Big]\\\space &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ yx - h_{\theta}(x)x] \Big]\\\space &=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}[ y - h_{\theta}(x)]x \Big]\\\therefore \space \dfrac{\partial J(\theta)}{\partial \theta_{j}}&=\dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}( h_{\theta}(x^{(i)}) - y^{(i)})x^{(i)}_{j}\Big]\\\end{aligned}$

ตามคอนเซปแล้ว cost จะต้องเป็นค่าตัวเลขหนึ่งๆ เพื่อการนำไปเทียวว่า cost ปัจจุบัน กับ cost อันที่ผ่านมา อันไหนมากกน้อยกว่ากัน. ในทางกลับกัน Gradient descent จะต้องมีจำนวนตัวแปรเท่ากับ $\theta$ เพื่อที่ว่าเราจะต้องเอาไปเพิ่มลดค่าของ $\theta$ แต่ละตัวเพื่อจะได้นำไปหา Cost ใหม่ในครั้งต่อไป. โค๊ดข้างล่างนี้คือการหาค่า Cost และ Gradient

	def sigmoid(z):
	return np.divide(1, np.add(1, np.exp(-z)))
	def hypothesis(theta, X):
	return sigmoid(X.dot(theta)).reshape((X.shape[0], 1))
	def cost_function(theta, X, y):
	"""
	Compute the cost of a particular choice of theta.
	"""
	training_example = X.shape[0]
	h_value = hypothesis(theta, X)
	cost = np.sum(-y * np.log(h_value) - (1 - y) * np.log(1 - h_value))
	return np.divide(cost, training_example)
	def gradient(theta, X, y):
	"""
	Calculate a derivative value of a particular set of theta
	"""
	m = X.shape[0]
	cost = np.sum((hypothesis(theta, X) - y) * X, axis=0)
	return np.divide(cost, m).reshape(X.shape[1], 1)
	def run():
	"""
	Entry point
	"""
	data = np.loadtxt("./data/ex2data2.txt", delimiter=",")
	plot_data(data, "Test 1 & 2 data")
	# Number of samples
	m = data.shape[0]
	# Sample length
	n = data.shape[1]
	X = data[:, 0:n - 1].reshape((m, n - 1))
	X = np.insert(X, 0, 1, axis=1)
	X = map_feature(X[:, 1], X[:, 2])
	y = data[:, n - 1].reshape((m, 1))
	print("X.shape", X.shape)
	nx = X.shape[1]
	initial_theta = np.zeros((nx, 1))
	cost = cost_function(initial_theta, X, y)
	grad = gradient(initial_theta, X, y)
	print("cost: ", cost)
	print("grad: ", grad)

view raw ThoughtLog_LRM03.py hosted with ❤ by GitHub

ณ จุดนี้ เราไม่เพียงได้ข้อมูลที่พร้อมแล้ว (Feature mapping) เรายังได้สมการที่พร้อมในการคำนวนอีก ที่เหลือก็แค่ใส่ลูปให้มัน หาCost-> หาGradient -> อับเดท$\theta$ วนไปซักพัก สุดท้ายเราก็จะได้ Cost ที่น้อยทีสุด และได้ $\theta$ ออกมาชุดหนึ่งเืพือนำไปใช้ทำนาย

โดยปกติแล้ว เราจะเขียนลูปอับเดท $\theta$ ของเราเอง ซึ่งเราต้องกำหนด จำนวนรอบในการ iterate และ learning rate ($\alpha$ ) ว่าพอได้ Gradient มาแล้ว เราจะทำการขยับเข้าใกล้จุด convex ในอัตราเร็วเท่าไหร่

$\theta_{j} := \theta_{j} - \alpha\dfrac{ \partial}{\partial \theta_{j}}J(\theta_{j})$

แต่ในที่นี้ เราจะมาลองใช้ Optimizaiton library ของ Scipy เพื่อมาทำงานนี้ให้เรา ซึ่งก็ควรจะได้ผลลัพท์ไม่ต่างจากที่เราเขียนเอง เพียงแต่ไม่ต้องระบุ learning rate และจำนวน iteration. ทั้งหมดนี้ทำได้โดยเพิ่มโค๊ดข้างล่างต่อจาก method run

	def cost_function_decorate(theta):
	return cost_function(theta, X, y)
	def gradient_decorate(theta):
	return gradient(theta, X, y).flatten()
	theta = fmin_ncg(cost_function_decorate, initial_theta, gradient_decorate)
	print("theta", theta)

view raw ThoughtLog_LRM04.py hosted with ❤ by GitHub

เหตุผลที่ต้องสร้าง cost_function_decorate กับ gradient_decorate ก็เพื่อให้สอดคล้องกับ signature ที่ fmin_ncg ต้องการ. เมื่อเราได้ $\theta$ มาแล้ว ก็สามารถเอาไปทำนายข้อมูลใหม่ๆได้ เช่นต่อไปนี้ ถ้าเรารู้ค่าที่ได้จากการทดสอบด่านที่หนึ่งและสอง ก็เอาสองค่านั้นไปเข้า map_feature ของเรา จากนั้นก็เอาไปใส่ $h_{\theta}(x)$ พร้อมกับชุดของ$\theta$ ที่ได้จากการเทรนด์ (hypothesis(theta, X)) ก็จะได้ค่าออกมา ถ้าผลลัพท์มากกว่า 0.5 ก็แน่นอนว่าชิพอันนี้ได้ผ่านเข้ารอบต่อไป

โดยหลักๆแล้วคอนเซปพื้นฐานก็มีเพียงแค่นี้ โค๊ดทั้งหมดตามนี้

	"""Logistic Regression Model"""
	import matplotlib.pyplot as plt
	import numpy as np
	from scipy.optimize import fmin_ncg
	def plot_data(data, window_title):
	"""
	:param data: np array of shape (n,3)
	:param window_title: String
	"""
	fig, ax = plt.subplots()
	fig.canvas.set_window_title(window_title)
	accepted = data[data[:, 2] == 1, :]
	rejected = data[data[:, 2] == 0, :]
	ax.scatter(accepted[:, 0], accepted[:, 1], c='blue', label='Accepted')
	ax.scatter(rejected[:, 0], rejected[:, 1], c='red', label='Rejected')
	ax.set_xlabel("Microchip Test 1")
	ax.set_ylabel("Microchip Test 2")
	ax.legend()
	ax.grid(True)
	plt.show(block=False)
	def map_feature(X1, X2, degree=6):
	"""
	Map Feature
	:param X1: NP array of shape (n,)
	:param X2: NP array of shape (n,)
	:param degree: integer denote the highest order
	:return:
	"""
	m = X1.shape[0]
	X1 = X1.T
	X2 = X2.T
	out = np.ones((m, 1), dtype=float)
	for i in range(1, degree + 1):
	for j in range(i + 1):
	value = ((X1 ** (i - j)) * (X2 ** j)).reshape(m, 1)
	out = np.append(out, value, axis=1)
	return out
	def sigmoid(z):
	return np.divide(1, np.add(1, np.exp(-z)))
	def hypothesis(theta, X):
	return sigmoid(X.dot(theta)).reshape((X.shape[0], 1))
	def cost_function(theta, X, y):
	"""
	Compute the cost of a particular choice of theta.
	"""
	training_example = X.shape[0]
	h_value = hypothesis(theta, X)
	cost = np.sum(-y * np.log(h_value) - (1 - y) * np.log(1 - h_value))
	return np.divide(cost, training_example)
	def gradient(theta, X, y):
	"""
	Calculate a derivative value of a particular set of theta
	"""
	m = X.shape[0]
	cost = np.sum((hypothesis(theta, X) - y) * X, axis=0)
	return np.divide(cost, m).reshape(X.shape[1], 1)
	def plot_contour(theta, X, y):
	"""Plot Boundary"""
	fig, ax = plt.subplots()
	accepted = X[y[:, 0] == 1, 1:3]
	rejected = X[y[:, 0] == 0, 1:3]
	ax.scatter(accepted[:, 0], accepted[:, 1], c='blue', label='Accepted')
	ax.scatter(rejected[:, 0], rejected[:, 1], c='red', label='Rejected')
	ax.set_xlabel("Microchip Test 1")
	ax.set_ylabel("Microchip Test 2")
	ax.legend()
	ax.grid(True)
	mc1 = X[:, 1]
	mc2 = X[:, 2]
	margin = 0.2
	min_mc1 = np.min(mc1) - margin
	min_mc2 = np.min(mc2) - margin
	max_mc1 = np.max(mc1) + margin
	max_mc2 = np.max(mc2) + margin
	delta = 0.025
	mc1_range = np.arange(min_mc1, max_mc1, delta)
	mc2_range = np.arange(min_mc2, max_mc2, delta)
	mc1_range, mc2_range = np.meshgrid(mc1_range, mc2_range)
	Z = np.zeros(mc1_range.shape)
	for i in range(mc1_range.shape[0]):
	Z[:, i] = map_feature(mc1_range[:, i], mc2_range[:, i]).dot(theta)
	CS = ax.contour(mc1_range, mc2_range, Z, colors='k')
	ax.clabel(CS, inline=1, fontsize=10)
	plt.show(block=False)
	def run():
	"""
	Entry point
	"""
	data = np.loadtxt("./data/ex2data2.txt", delimiter=",")
	plot_data(data, "Test 1 & 2 data")
	# Number of samples
	m = data.shape[0]
	# Sample length
	n = data.shape[1]
	X = data[:, 0:n - 1].reshape((m, n - 1))
	X = np.insert(X, 0, 1, axis=1)
	X = map_feature(X[:, 1], X[:, 2])
	y = data[:, n - 1].reshape((m, 1))
	print("X.shape", X.shape)
	nx = X.shape[1]
	initial_theta = np.zeros((nx, 1))
	cost = cost_function(initial_theta, X, y)
	grad = gradient(initial_theta, X, y)
	print("cost: ", cost)
	print("grad: ", grad)
	def cost_function_decorate(theta):
	return cost_function(theta, X, y)
	def gradient_decorate(theta):
	return gradient(theta, X, y).flatten()
	trained_theta = fmin_ncg(cost_function_decorate, initial_theta, gradient_decorate)
	print("theta", trained_theta)
	plot_contour(trained_theta, X, y)
	if __name__ == '__main__':
	run()

view raw ThoughtLog_LRM05_final.py hosted with ❤ by GitHub

พอรันก็จะได้แบบนี้

ถ้าคิดว่ากำลังจะได้เห็นบทลงท้าย ก็ต้องขออภัยมา ณ ที่นี้ด้วย เพราะเรื่องยังไม่จบเพียงเท่านี้ ใครอ่านมาถึงตรงนี้ก็จะรู้ว่ายิ่งเขียนยิ่งรั่วนะครับฮ่าๆๆ หากสังเกตดีๆ กราฟที่เราได้มานี้มัน Overfitting. มากเกินไป คือมันอาจทำนายผลลัพท์จากข้อมูลที่เราใช้เทรนได้แม่นยำก็จริง แต่มันอาจจะแย่เอามากๆกับข้อมูลที่มันไม่เคยเห็น ดังภาพสมมุติด้านล่าง

สีสเปร์ยที่พ่นไปแดงๆ อาจจะเป็นข้อมูลอีกหลายพันชิ้นที่เรายังไม่เคยเห็นก็ได้ แต่เผอิญโมเดลของเราไป fit กับ training data มากเกินไป มันก็ไม่ดี คำถามถัดมาคือ เราจะทำการแก้ปัญหานี้ได้ยังไง

Regularized Linear Regression

คือการเพ่ิ่ม Regularization เข้ากับ Cost function เพื่อไปคานอำนาจกับ $\theta$ ไม่ให้ข้อมูลที่เอามาใช้เทรนด์ไปมีผลกระทบต่อ $\theta$ มากเกินไป เลยส่งผลให้ $\theta$ ค่อนข้าง fit กับข้อมูลในแบบภาพรวมมากกว่าการพยายามเจาะจงเก็บทุกๆเม็ด

$J(\theta)=- \dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}y^{(i)}log(h_{\theta}(x^{(i)})) + (1-y)log(1-h_{\theta}(x^{(i)})) \Big] + \dfrac{\lambda}{2m}\sum^n_{j=1}\theta^2_j$

ซึ่งก็ทำให้ Gradient decent กลายมาเป็น

$\dfrac{\partial J(\theta)}{\partial \theta_{j}}=\Big(\dfrac{1}{m}\Big[\displaystyle\sum^{m}_{i=1}( h_{\theta}(x^{(i)}) - y^{(i)})x^{(i)}_{j}\Big]\Big) + \dfrac{\lambda}{m}\theta_j \space \text{for j}\ge 1$

ส่วนคอนเซปเกือบทั้งหมดก็ยังเหมือนเดิม โค้ดชุดนี้ไม่ค่อยสะอาดเท่าไหร่นะครับ ลองโฟกัสที่ค่าของ ramda และสังเกตผลลัพท์ของโมเดลที่ได้ Regularization Code

Reference: https://www.coursera.org/learn/machine-learning