Numpy跟Linear regression 多維度篇 « Terrence的宅宅幻想
這篇是接續使用numpy跟sympy實作Linear regression所做的筆記，在上篇中以simple linear regression為中心做了一些學習．這次想針對多維度的線性回歸演算法做個學習筆記，當中主要參考臺大林軒田教授：機器學習的基石課程內容．
簡單回顧一下，假設存在一筆資料(x, y)，其回歸方程式可用 來表示，推估資料跟實際資料的殘差平方和可以用
來表示，對於尋找一個回歸方程式的解我們希望他的殘差平方和要最小．
在簡單線性回歸的時候，因為只有一個自變數，因此回歸方程式可以簡單表示成y=ax+b，但是在多維的時候會有超過一個以上的feature，在這邊要稍微改變一下表達方式，假設今天有n個feature對於每個變數x，給予其相對應的權重w，新的回歸方程式會長的如下
就是推估值變成w1*x1 + w2*x2 + ... + wn*xn 這樣的形式，多個變數要怎麼套用使用numpy跟sympy實作Linear regression的方法求解呢，上篇提到的
1. 算殘差和 2. 對每個自變數作微分 3. 解聯立方程式
國中數學時告訴我們一個變數需要一個方程式求解，兩個變數需要兩個方程式求解，n個變數就需要n個方程式求解
也就是n個變數我們最少就需要n筆資料才能用聯立方程式找出解，當資料少於n的時候怎麼辦呢，這時候就要用到線性代數的性質了
假設有m筆資料，每筆資料feature長度是n，他的可以表達成如下的聯立方程式
其中b是一個常數，在做線性回歸演算法的時候這個常數很重要，他相當於簡單線性回歸時y=ax+b的常數角色，之前幾次實驗都忘了把常數加進去造成結果有相當的落差．
最後根據線性代數的推導，我們可以把input資料全部整理成一個大大矩陣X，把feature跟輸出結果y整理成一個長長的向量W跟Y
然後我們的殘差和經過各種收納推導可以表示成如下
這邊簡化了一些內容，實際的推導過程請參考機器學習的基石課程內容，當初被推導過程卡了很久，但是總結來說當前階段只需要先理解結論的用法也就夠了，上面的數學公式經過一些推演之後可以得到殘差的最小平方和的解可以用如下表示
在計算多維的線性回歸的時候不再需要辛苦的列舉聯立方程式，只要
1. 將輸入資料(自變數)做成矩陣Ｘ(要記得補上常數項) 2. 將輸出y做成向量Y 3. 帶入公式即可很快的求解
不過上面公式有一個弔詭的地方就是需要“反矩陣”運算，初級線性代數告訴我們在某些情況下反矩陣才會存在
也就是前面提到的資料量m大於feature數量n的時候，當反矩陣不存在該公式則無解
為了應付這種問題我們需要使用廣義反矩陣(pseudo inverse)
http://mathworld.wolfram.com/Moore-PenroseMatrixInverse.html
X的廣義反矩陣符號可以用加號表示
公式修改成如下:
回到Numpy實作linear regression的問題，在numpy除了使用numpy.linalg.lstsq之外可以很簡單的用上面的公式自己實作多維線性回歸，這邊用兩個feature也就是三維空間做個簡單的例子
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #準備資料，以x,y為自變數，z為應變數 x = [45.9, 41.3, 10.8, 48.9, 32.8, 19.6, 2.1, 2.6, 5.8, 24, 35.1, 7.6, 32.9, 39.6, 20.5, 23.9, 27.7, 5.1, 15.9, 16.9, 12.6, 3.5, 29.3, 16.7, 27.1, 16, 28.3, 17.4, 1.5, 20, 1.4, 4.1, 43.8, 49.4, 26.7, 37.7, 22.3, 33.4, 27.7, 8.4, 25.7, 22.5, 9.9, 41.5, 15.8, 11.7] y = [69.3, 58.5, 58.4, 75, 23.5, 11.6, 1, 21.2, 24.2, 4, 65.9, 7.2, 46, 55.8, 18.3, 19.1, 53.4, 23.5, 49.6, 26.2, 18.3, 19.5, 12.6, 22.9, 22.9, 40.8, 43.2, 38.6, 30, 0.3, 7.4, 8.5, 5, 45.7, 35.1, 32, 31.6, 38.7, 1.8, 26.4, 43.3, 31.5, 35.7, 18.5, 49.9, 36.8] z = [9.3, 18.5, 12.9, 7.2, 11.8, 13.2, 4.8, 10.6, 8.6, 17.4, 9.2, 9.7, 19, 24.4, 11.3, 14.6, 18, 12.5, 5.6, 15.5, 9.7, 12, 15, 15.9, 18.9, 10.5, 21.4, 11.9, 9.6, 17.4, 9.5, 12.8, 25.4, 14.7, 10.1, 21.5, 16.6, 17.1, 20.7, 12.9, 8.5, 14.9, 10.6, 23.2, 14.8, 9.7] fig = plt.figure() ax = fig.add_subplot(111, projection='3d') #將輸入資料(自變數)做成矩陣Ｘ X = np.mat([x, y]).T #補上常數 X = np.insert(X, 0, np.ones(46), 1) #將輸出結果做成向量Y Y = np.mat([z]).T #np.linalg.lstsq(X,Y) #用於比對結果np.linalg.lsts #套用公式：廣義反矩陣 * 輸出向量Y，求得解 W = X.I * Y #製作程回歸方程式 y^= b + w1x1 + w2x2 LR = lambda x, y: W[0, 0] + W[1, 0] * x + W[2, 0] * y #畫圖輸出結果 x2 = y2 = np.arange(0, 60.0) X, Y = np.meshgrid(x2, y2) zs = np.array([LR(x2,y2) for x2,y2 in zip(np.ravel(X), np.ravel(Y))]) Z = zs.reshape(X.shape) ax.plot_surface(X, Y, Z) ax.scatter(x, y, z) ax.set_xlabel('X Label') ax.set_ylabel('Y Label') ax.set_zlabel('Z Label') plt.show()
跑出來的結果如下
