使用numpy跟sympy實作Linear regression « Terrence的宅宅幻想
這一陣子在上機器學習的課，對於線性回歸（linear regression）的演算法一直不是很理解
這幾天嘗試用numpy去實作才總算真的能夠稍微理解，今天這裡就用簡單線性回歸（Simple linear regression）為例來做個筆記
線性回歸(linear regression)
個人的理解就是想要找出一群資料(Data Set)的線性關係(回歸方程式)，以對未來未知的新進資料進行推估或預測．用在機器學習領域則可以做簡單的二元分類．做法是把資料丟到回歸方程式看出來的結果是否大於某個門檻值(threshold)
用簡單線性回歸作說明，假設有一組(x, y)資料代表散落在平面上的一堆點
拉回正題，簡單線性回歸因為只有一個feature（或稱自變數），所以其回歸方程式可以再二維平面上簡單表示成
這回歸方程式可以表示x與y之間的線性關係的一條線，可以用來推估新的x進來的時候，他的y可能是多少
現在問題來了，這條線(方程式)怎麼來的
網路上還蠻容易找到已經套好的“公式解”直接求a跟b的值，但是我想更深入探討求a跟b過程的精神
參考上圖每筆資料的應變數y跟用回歸方程式y=ax+b求出來的y之間都存在著一個距離一般稱之為殘差（residual）
對於我們來說希望這條線到所有點的距離總和也就是殘差的和要最小，這樣這條線就最能代表這群資料的線性關係
對於這種求差距最小的方法常用的是最小平方法（least squares method）
最小平方法（least squares method）
資料(x, y)跟回歸方程式求出來的點(x, y^=ax+b)的差距(y - y^)可以用以下公式表示
而因為是“差”所以有可能有正有負，因此取平方保證差距是個正值以確保之後加總
在此之前我一直有個疑問，為什麼不用絕對值就好呢？
查了一些資料跟複習了大一微積分之後才多少有點領略
使用平方（squares）取正值有兩個原因
1. 放大差距 2. 可以利用微積分的性質替我們求最佳解
針對第一點，統計學上在某些情況下差距只取絕對值的時候會難以看出資料間的離散程度，但是透過平方除了可以取正值之外最重要的一個因素就是“放大差異
而最後我們可以得到所有點的殘差和(sum residual)的公式如下
對於回歸方程式的系數a,b我們的最終目標是讓上面這條公式算出來的差距總和“最小”，這樣就確保y=ax+b這條線到所有點的距離總和最小
而要怎麼讓上面這個公式有最小值呢？以前的時候我百思不得其解，回頭翻大一微積分的時候終於比較能理解做法跟使用微積分的原因，答案是用微積分的“梯度（Gradient）”
在談梯度之前先聊聊為什麼用微積分，在做線性回歸(linear regression)的時候由於求的解是一條直線，如果單純用絕對值去算殘差和的話他的表示式仍是直線
直線沒辦法知道極限值(Limit)在哪裡，要求最小值(也可以說是極限值)比較難以有直觀的解．但方程式經過平方（squares）後就不一樣了．它會變成一條曲線．
有曲線之後就會有山谷出現，有山谷自然就會有最低點或最高點(極限值)，透過微分可以簡單替我們找出一條曲線什麼時候會有最低點
舉個簡單例子，參考以下公式及其代表曲線
上面這條曲線的最低點（極限值），透過微分方程式＝0的時候可以找到極限值這個性質去找最低點
微分的解是x=1的時候我們可以找到該曲線的極限值
拉回正題，回頭看我們的殘差和公式
上面這條公式要如何求最佳解呢，套用梯度的定義，我們分別針對a跟b作微分可以得到梯度向量
我們將得兩組微分方程式可以對a,b解聯立方程式求得a跟b的最佳解，這就是最小平方法了
舉個實際的例子，假設有三個點(1, 2), (2, 4), (3, 3)，求代表這三個點的直線回歸方程式
上面三筆資料帶入殘差和的公式可以得到
展開後可以得到
現在的目標就是要求上面這條公式帶入什麼樣的a跟b會有最小值，現在分別對a跟b微分可以到兩條方程式
對a微分
對b微分
只要將上面兩式子分別等於0就可以用聯立方程式求出a,b
最後得到a=0.5跟b=2這條線就是我們想找的最佳解
使用numpy跟sympy實作
用numpy可以方便我們作數值運算，用sympy則可以替我們處理方程式的計算
首先產生測試資料，用y=x這條方程式為基準去隨機產生一些點，上下震盪幅度約為2
import numpy as np import matplotlib.pyplot as plt #Target function line X = np.linspace(0, 10) f = lambda x: x #y=x F = np.vectorize(f) Y = F(X) #random data by F(X) + random residual(upper bound=2) num = 15 #number of data random_sign = np.vectorize(lambda x: x if np.random.sample() > 0.5 else -x) data_X = np.linspace(1, 9, num) data_Y = random_sign(np.random.sample(num) * 2) + F(data_X) plt.plot(X, Y, 'b') # render blue line plt.plot(data_X, data_Y, 'ro') # render data point plt.show()
上面的結果大概會長的像下面一樣
簡單解釋一下linspace這個函式會在一個區間產生連續的一群資料，除了可以當產生資料的工具之外
用在畫圖上也可以用來畫方程式線，vectorize是把一個把python函式包裝成一個可以對numpy的資料進行運算的新function
>>> import numpy as np
>>> X = np.linspace(0, 10)
>>> X
array([ 0. , 0.20408163, 0.40816327, 0.6122449 ,
0.81632653, 1.02040816, 1.2244898 , 1.42857143,
1.63265306, 1.83673469, 2.04081633, 2.24489796,
2.44897959, 2.65306122, 2.85714286, 3.06122449,
3.26530612, 3.46938776, 3.67346939, 3.87755102,
4.08163265, 4.28571429, 4.48979592, 4.69387755,
4.89795918, 5.10204082, 5.30612245, 5.51020408,
5.71428571, 5.91836735, 6.12244898, 6.32653061,
6.53061224, 6.73469388, 6.93877551, 7.14285714,
7.34693878, 7.55102041, 7.75510204, 7.95918367,
8.16326531, 8.36734694, 8.57142857, 8.7755102 ,
8.97959184, 9.18367347, 9.3877551 , 9.59183673,
9.79591837, 10. ])
>>> f = lambda x: x + 1 #把X所有的值+1當作Y的值
>>> F = np.vectorize(f)
>>> Y = F(X)
>>> Y
array([ 1. , 1.20408163, 1.40816327, 1.6122449 ,
1.81632653, 2.02040816, 2.2244898 , 2.42857143,
2.63265306, 2.83673469, 3.04081633, 3.24489796,
3.44897959, 3.65306122, 3.85714286, 4.06122449,
4.26530612, 4.46938776, 4.67346939, 4.87755102,
5.08163265, 5.28571429, 5.48979592, 5.69387755,
5.89795918, 6.10204082, 6.30612245, 6.51020408,
6.71428571, 6.91836735, 7.12244898, 7.32653061,
7.53061224, 7.73469388, 7.93877551, 8.14285714,
8.34693878, 8.55102041, 8.75510204, 8.95918367,
9.16326531, 9.36734694, 9.57142857, 9.7755102 ,
9.97959184, 10.18367347, 10.3877551 , 10.59183673,
10.79591837, 11. ])
現在開始進行linear regression的實作、如下
#using sympy from sympy import * def linear_regression(X, Y): a, b = symbols('a b') residual = 0 #residualSum 求殘差 for i in range(num): residual += (Y[i] - (a * X[i] + b)) ** 2 #展開觀看方程式內容 print expand(residual) #對a微分 f1 = diff(residual, a) #對b微分 f2 = diff(residual, b) print f1 print f2 #求聯立方程式的解 res = solve([f1, f2], [a, b]) return res[a], res[b] a, b = linear_regression(data_X, data_Y)
用sympy的symbols可以替我們產生可以當作方程式未知數的“符號”並且可以進行數學幾何運算
之後再用solve就可以很簡單的算出方程式=0的解
之後再用求到的a,b畫出回歸方程式的線
LR_X = X h = lambda x: a*x + b H = np.vectorize(h) LR_Y = H(LR_X) plt.plot(X, Y, 'b') # render blue line plt.plot(LR_X, LR_Y, 'g') # render green line plt.plot(data_X, data_Y, 'ro') plt.show()
其結果會長得像如下
藍色的線是我們用來產生資料的原始方程式，綠色的線則是經由回歸算出來的回歸方程式
由於資料點不多，回歸方程式跟產生資料的方程式差的有點多
這就牽扯到機器學習跟統計的問題了，資料要夠多才能誤差比較小預測比較精準
下面一個資料比較多的結果，可以看出來回歸方程式已經非常接近原始產生資料的方程式了
這次只用簡單線性回歸作介紹，但上面這個方法可以在套用到多個feature，不管有多少個自變數方法都一樣
1. 算殘差和 2. 對每個自變數作微分 3. 解聯立方程式
希望下次有機會再用多變數進行分享
A = np.vstack([data_X, np.ones(len(data_X))]).T a, b = np.linalg.lstsq(A, data_Y)[0]