系列文章目錄：
如果說感知機是最最最簡單的分類算法，那么線性回歸就是最最最簡單的回歸算法，所以這一篇我們就一起來快活的用兩種姿勢手擼線性回歸吧；
算法介紹
線性回歸通過超平面擬合數據點，經驗誤差一般使用MSE（均平方誤差），優化方法為最小二乘法，算法如下：
- 假設輸入數據為X，輸出為Y，為了簡單起見，這里的數據點為一維數據（更好可視化，處理方式沒區別）；
- MSE公式為：\(\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)；
- 最小二乘法：最小指的是目標是min，二乘指的就是MSE中誤差的二次方，公式為：\(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)；
- 由於目標是查找擬合最好的超平面，因此依然定義變量w和b；
- 對於w和b的求解有兩種方式：
- 列出最小化的公式，利用優化求解器求解：
- 基於已知的X、Y，未知的w、b構建MSE公式；
- 定義最小化MSE的目標函數；
- 利用求解器直接求解上述函數得到新的w和b；
- 對經驗誤差函數求偏導並令其為0推導出w和b的解析解：
- 基於最小化MSE的優化問題可以直接推導出w和b的計算方法；
- 基於推導出的計算方法直接計算求解；
- 列出最小化的公式，利用優化求解器求解：
利用求解器求解
利用求解器求解可以看作就是個列公式的過程，把已知的數據X和Y，未知的變量w和b定義好，構建出MSE的公式，然后丟到求解器直接對w和b求偏導即可，相對來說代碼繁瑣，但是過程更簡單，沒有任何數學推導；
代碼實現
初始化數據集
X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])
定義變量符號
所謂變量指的就是那些需要求解的部分，次數就是超平面的w和b；
w,b = symbols('w b',real=True)
定義經驗誤差函數MSE
RDh = 0
for xi,yi in zip(X,y):
RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)
定義求解函數
此處就是對w和b求偏導；
eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)
求解w和b
ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]
運行結果
完整代碼
from sympy import symbols, diff, solve
import numpy as np
import matplotlib.pyplot as plt
'''
線性回歸擬合wx+b直線；
最小二乘法指的是優化求解過程是通過對經驗誤差（此處是均平方誤差）求偏導並令其為0以解的w和b；
'''
# 數據集 D X為父親身高，Y為兒子身高
X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])
# 構造符號
w,b = symbols('w b',real=True)
# 定義經驗誤差計算公式：(1/N)*sum(yi-(w*xi+b))^2)
RDh = 0
for xi,yi in zip(X,y):
RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)
# 對w和b求偏導：求偏導的結果是得到兩個結果為0的方程式
eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)
# 求解聯立方程組
ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]
print('使得經驗誤差RDh取得最小值的參數為：'+str(ans))
plt.scatter(X,y)
x_range = [min(X)-0.1,max(X)+0.1]
y_range = [w*x_range[0]+b,w*x_range[1]+b]
plt.plot(x_range,y_range)
plt.show()
推導公式求解
與利用優化器求解的區別在於針對\(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)對\(w\)和\(b\)求偏導並令其為0，並推導出w和b的計算公式是自己推導的，還是由優化器完成的，事實上如果自己推導，那么最終代碼實現上會非常簡單（推導過程不會出現在代碼中）；
w和b的求解公式推導
首先，我們的優化目標為：
\[min \frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]
去除公式中無關的常量部分：
\[min \sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]
由於一般w是向量，而b為標量，因此通常會將w和b組成[w b]，x變為[x 1]來統一處理w和b，調整后如下：
\[\sum_{i=1}^{N}(wx_i^T-y_i)^2 \]
上式把平方拆開有：
\[\sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2) \]
對於w（注意此時w為原w和b的組合）求偏導過程如下：
\[\begin{equation*} \begin{split} \frac{\partial \sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2)}{\partial w} &= 2w^Tx_ix_i^T-2x_i^Ty_i \\ &= 0 \\ \end{split} \end{equation*} \]
上式變形后有：
\[2w^Tx_ix_i^T - 2x_i^Ty_i = 0 \\ w^Tx_ix_i^T = x_i^Ty_i \\ w^T = (x_ix_i^T)^{-1}x_i^Ty_i \]
由於此處的w其實是w和b的組合，因此通過這一次推導就得到了w和b兩個求解方法；
代碼實現
構造數據集
X = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])
為X增加元素全為1的一列用於和b的計算
ones = np.ones(X.shape[0]).reshape(-1,1)
X = np.hstack((ones,X))
通過求解公式求解w和b
w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
w,b = w[1:],w[0]
運行結果
完整代碼
完整代碼對於求解部分使用的是偽逆而不是逆，原因在於求解公式中正好構造了偽逆，而偽逆適用性強國求逆，因此使用偽逆代替逆；
import numpy as np
import matplotlib.pyplot as plt
rnd = np.random.RandomState(3) # 為了演示，采用固定的隨機
'''
單變量線性回歸最小二乘法的矩陣實現：矩陣實現的優勢在於numpy本身支持偽逆；
其實就是對於誤差平方和的矩陣形式對於W求導並令其為0，得到w_hat = (X^T*X)^-1*X^T*Y，其中(X^T*X)^-1*X^T稱為偽逆（pseudo inverse，即函數pinv）
因此可以省略中間大量的構造經驗誤差、解偏導方程組等步驟；
'''
class LinearRegression(object):
def __init__(self,X,y):
ones = np.ones(X.shape[0]).reshape(-1,1) # 1用於計算b
self.X = np.hstack((ones,X))
self.y = y
def train(self):
# 注意，雖然一般情況下下面二者是等價的，但是在矩陣無法求逆或某些其他情況下時，二者並不相等
# 相對而言偽逆定義更加寬泛，用處更廣，因此可以的情況下建議使用偽逆
# self.w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
self.w = np.linalg.pinv(self.X) @ self.y
self.w = self.w.reshape(-1)
self.w,self.b = self.w[1:],self.w[0]
return self.w,self.b
def predict(self,x):
return self.w.dot(x)+self.b
def get(self):
return self.X,self.y,self.w,self.b
if __name__ == '__main__':
X0 = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])
model = LinearRegression(X=X0,y=y)
w,b = model.train()
print(f'最小二乘法的矩陣方式結果為：w={w} b={b}')
print(model.predict(np.array([X0[0]])))
plt.scatter(X0,y)
plt.plot([min(X0),max(X0)],[model.predict(min(X0)),model.predict(max(X0))])
plt.show()
最后
對於算法的學習，一定的數學知識是必要的，對於公式的推導可以讓我們對於算法內部運行邏輯有更深的了解；