Patchouli的機器學習系列教程三：線性回歸——『器』

02-02

我們這次要實踐的內容就是上節課中提供的馬拉松數據：Patchouli的機器學習系列教程三：線性回歸——『道』

import pylab as pltnimport numpy as npnimport pandas as pdnfrom pandas import DataFramenb=pd.read_csv(marathon.csv)nx=DataFrame(b[year])ny=DataFrame(b[hour])nplt.plot(x, y, rx)n

導入並展示數據。

1.無腦坐標下降法

如往常一樣，設計目標函數，隨便定一組參數，不斷優化至目標函數最小的點。不過與電影推薦不同的是我們不再按照梯度一點點下降，而是直接令目標函數的偏導數等於零，求出參數的最大似然值，這種方法被稱為坐標下降法，也就是在『道』篇中的 $c^* ,m^*,(sigma^2)^*$ 。

其中 $c^* = frac{sum_{i=1}^n (y_i - m^* x_i)}{n}$

我們隨機設定一個差不多的 $m$ ，然後就可以算出 $c^*$ 了

m = -0.4nc = (y - m*x).mean()nprint(c)n

再根據算出的 $c^*$ ，用同樣的方式計算出 $m^*$ 。

$m^* = frac{sum_{i=1}^n (y_i - c)x_i}{sum_{i=1}^n x_i^2}$

並直接擬合出回歸直線：

m = ((y - c)*x).sum()/(x**2).sum()nprint(m)nx_test = np.linspace(1890, 2020, 130)[:, None]nf_test = m*x_test + cnplt.plot(x_test, f_test, b-)nplt.plot(x, y, ro)n

emmmmmmmmmmmm

我們發現擬合的效果並不是很好，畢竟是隨機寫的參數，擬合成這樣也是沒辦法。

我們增加一些擬合次數：

for i in np.arange(10):n m = ((y - c)*x).sum()/(x*x).sum()n c = (y-m*x).sum()/y.shape[0]nf_test = m*x_test + cnplt.plot(x_test, f_test, b-)nplt.plot(x, y, rx)n

結果好像還是不太好，目標函數的值還是很大，那麼我們就以目標函數的變化率為基礎，直到目標函數找到最小值為止（畢竟不可能真的等於零）：

f,ax=plt.subplots(figsize=(8,8))nax.plot(x,y,r*,markersize=8)nm=-0.4nc=(y-m*x).mean()nE1=10000 #好大好大的數nE2=((y-m*x-c)**2).sum()nwhile ((E1-E2)/E1)>0.0001: n E1=E2n for i in np.arange(100):n m = ((y-c)*x).sum()/(x*x).sum()n c = (y-m*x).sum()/y.shape[0]n E2=((y-m*x-c)**2).sum()n print("after iteration",(((y-m*x-c)**2).sum()))nx_test = np.linspace(1890, 2020, 131)[:,None]nf_test = m*x_test + cnax.plot(x_test,f_test,b-,markersize=8)n

在經過大約26000次迭代後，我們終於擬合出了一條符合規律的直線，以上就是無腦線性回歸的基本知識，接下來我們繼續改進我們的演算法。

2.特解法

我之前提到過我們可以把線性方程 $y=mx+c$ 寫成向量的內積形式：

令 $mathbf{x} = begin{bmatrix} 1 x end{bmatrix}， mathbf{w} = begin{bmatrix} c m end{bmatrix}$

則 $mathbf{x}cdotmathbf{w} = 1 times c + x times m = mx + c$

那麼我們便可以設計出 $X$ 這個包含所有隨機變數的矩陣：

$mathbf{X} = begin{bmatrix} mathbf{x}_1^top mathbf{x}_2^top vdots mathbf{x}_n^top end{bmatrix} = begin{bmatrix} 1 & x_1 1 & x_2 vdots & vdots 1 & x_n end{bmatrix}$

w = np.zeros(shape=(2,1))nw[0] = cnw[1] = mnX = np.hstack((np.ones_like(x), x))n

則： $E(m,c)= E(mathbf{w}) = sum_{i=1}^n (y_i - f(mathbf{x}_i; mathbf{w}))^2$

其中 $f(mathbf{x}_i; mathbf{w}) = mathbf{x}_i^top mathbf{w}$ ；那麼 $mathbf{f}(mathbf{x}; mathbf{w}) = begin{bmatrix} f(mathbf{x}_1; mathbf{w}) f(mathbf{x}_2; mathbf{w}) vdots f(mathbf{x}_n; mathbf{w}) end{bmatrix}$

$y_i$ 也可以被表示成 $mathbf{y} = begin{bmatrix} y_1 y_2 vdots y_n end{bmatrix}$

因而 $原式=(mathbf{y}-mathbf{f})^2$

由於有如下矩陣乘法定理：

$a = sum_{i=1}^k b_i^2 = mathbf{b}^top mathbf{b}$

因此我們的誤差函數最終可以寫成： $E(mathbf{w}) = (mathbf{y} - mathbf{f})^top (mathbf{y} - mathbf{f})$

其中 $mathbf{f} = mathbf{X}mathbf{w}$

f = np.dot(X, w)nresid = (y-f)nE = np.dot(resid.T, resid)n

目標函數由此得到了計算，我們終於迎來了最後的問題，也就是『道』篇最後一個公式

$mathbf{X}^top mathbf{X} mathbf{w} = mathbf{X}^top mathbf{y}$

我們可以使用numpy自帶的solve解方程函數自動解方程：

w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))n

這時候我們再把變數二維集合 $mathbf w$ 帶入原方程中：

m = w[1]; c=w[0]nf_test = m*x_test + cnplt.plot(x_test, f_test, b-)nplt.plot(x, y, rx)n

得到了和上面無腦回歸法一樣的結果。可以看到，我們包括『道』篇那麼多數學推理，拋去圖形展示的代碼，真正的代碼一共連十行都沒有。這就是我曾經說的，人工智慧程序員並不需要寫多少代碼，但是一行千金，在傳統碼農變成勞動密集型行業之後，人工智慧作為知識密集型產業依然保留著上古時期程序員們的風采。

我們本文用到的線性回歸方法同樣可以應用回昨天的電影數據，觀察一個電影不同的屬性對其評分的影響，比如年份，死人數，電影長度等等：

movies=DataFrame(pd.read_csv(movies.csv))nX = movies[[Year, Body_Count, Length_Minutes]]nX[Eins] = 1 ny = movies[[IMDB_Rating]]nw = pd.DataFrame(data=np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y)),index = X.columns,columns=[regression_coefficient]) n

其中w是這樣的：

可以發現，IBDM的評分保持著很大的穩定性，整體的波動並不大。

其中與年份成負相關，越晚出現的電影評分越低；

與死亡人數成負相關，死亡人數越多，評分越低，但是影響程度不大；

與電影長度成正比，電影越長，評分越高。

這份數據由於樣本太少，參考意義其實並不大，假如有能拿到大量數據的朋友不妨以此法檢驗一下（順便給我來一份），本教程僅負責教授方法。

學有餘力的同學可以繼續向下看：

在計算特徵向量時用到的方程是這個：

$mathbf{X}^top mathbf{X} mathbf{w} = mathbf{X}^top mathbf{y}$

這個方程簡明扼要，但是在程序運行環節還是有些太冗餘了，畢竟要把： $mathbf{X}^top mathbf{X}$ 算出來是要把 $mathbf X$ 的所有元素平方，再除到等式那邊去，不僅浪費算力，還會降低數據的精確度。精通數學的人會很快想到一種替代的解決辦法：QR因式分解。

QR是一種常用的矩陣分解手段，它可以將一個矩陣拆解成一個正交矩陣 $Q$ 和一個上三角矩陣 $R$ ，QR分解也因此得名。正交矩陣的特性 $Q^top Q=I$ 使得我們可以省略掉好多討厭的計算環節，況且這個也不消你算，numpy提供給你既有的函數了，我們根據QR分解，可以做出如下推導：

$begin{align*} mathbf{X}^top mathbf{X} boldsymbol{mathbf w} & = mathbf{X}^top mathbf{y} (mathbf{Q}mathbf{R})^top (mathbf{Q}mathbf{R}) boldsymbol{mathbf w} & = (mathbf{Q}mathbf{R})^top mathbf{y} mathbf{R}^top (mathbf{Q}^top mathbf{Q}) mathbf{R} boldsymbol{mathbf w} & = mathbf{R}^top mathbf{Q}^top mathbf{y} mathbf{R}^top mathbf{R} boldsymbol{mathbf w} & = mathbf{R}^top mathbf{Q}^top mathbf{y} mathbf{R} boldsymbol{mathbf w} & = mathbf{Q}^top mathbf{y} end{align*}$

根據最後的結果： $mathbf{R} boldsymbol{mathbf w} = mathbf{Q}^top mathbf{y}$ 再用解方程函數處理，就可以得到特徵向量 $mathbf w$ 了

import scipy as spnQ, R = np.linalg.qr(X)nw = sp.linalg.solve_triangular(R, np.dot(Q.T, y))nw = pd.DataFrame(w, index=X.columns)n

得到的是一樣的結果。