线性回归

线性回归，就是已知一系列x和y对应的点，通过求出$y=wx+b$（线性，所以是一条直线）去拟合数据点，预测某一个$x_0$对应的$y_0$是多少。

x = np.array([56, 72, 69, 88, 102, 86, 76, 79, 94, 74])
y = np.array([92, 102, 86, 110, 130, 99, 96, 102, 105, 92])

那么，如何求出这条直线？如何判断这条直线对数据的拟合程度好坏？
这里需要引入损失函数。

损失函数

平均绝对误差（MAE）

平均绝对误差（MAE）就是绝对误差的平均值，它的计算公式如下：
$$\text{MAE}(y,\hat{y})=\frac{1}{n}\sum _{i=1}^n|y_i-{\hat{y}}_i|\text{\hspace{1em}(1)}$$
其中，$y_i$ 表示真实值，${\hat{y}}_i$ 表示预测值，$n$ 则表示值的个数。MAE 的值越小，说明模型拥有更好的拟合程度。

def mae_value(y_true, y_pred):
    n = len(y_true)
    mae = sum(np.abs(y_true - y_pred))/n
    return mae

均方误差（MSE）

均方误差（MSE）表示误差的平方的期望值，它的计算公式如下：
$$\text{MSE}(y,\hat{y})=\frac{1}{n}\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2\text{\hspace{1em}(2)}$$
其中，$y_i$ 表示真实值，${\hat{y}}_i$ 表示预测值，$n$ 则表示值的个数。MSE 的值越小，说明预测模型拥有更好的精确度。

def mse_value(y_true, y_pred):
    n = len(y_true)
    mse = sum(np.square(y_true - y_pred))/n
    return mse

在这里，我们已经知道了如何求损失，但是如何才能让损失最小呢？

最小二乘法

最小二乘法是用于求解线性回归拟合参数 $w$ 的一种常用方法。最小二乘法中的「二乘」代表上面的均方误差（MSE），即均方误差最小。

最小二乘法代数推导

均方误差函数为：
$$\begin{gathered} \widehat{y_i}=w{x}_i+b \\ f=\sum _{i=1}^n{(y_i-(w{x}_i+b))}^2\text{\hspace{1em}(3)} \end{gathered}$$
这里要求f的最小值，故求偏导如下：
$$\begin{gathered} \frac{\partial f}{\partial b}=\frac{\partial {\sum }_{i=1}^n(y_i^2-2y_i(w{x}_i+b)+(w{x}_i+b{)}^2)}{\partial b} \\ =\sum _{i=1}^n(-2y_i+2w{x}_i+2b) \\ =-2(\sum _{i=1}^n{y}_i-nb-w\sum _{i=1}^n{x}_i)\text{\hspace{1em}(4a)} \end{gathered}$$
$$\begin{gathered} \frac{\partial f}{\partial w}=\frac{\partial {\sum }_{i=1}^n(y_i^2-2y_i(w{x}_i+b)+(w{x}_i+b{)}^2)}{\partial w} \\ =\sum _{i=1}^n(-2y_i{x}_i+2w{x}_i^2+2x_{i}b) \\ =-2(\sum _{i=1}^n{x}_i{y}_i-b\sum _{i=1}^n{x}_i-w\sum _{i=1}^n{x_i}^2)\text{\hspace{1em}(4b)} \end{gathered}$$
令 $\frac{\partial f}{\partial b}=0$ 以及 $\frac{\partial f}{\partial w}=0$，解得：
$$w=\frac{n{\sum }_{}^{}{x}_i{y}_i-{\sum }_{}^{}{x}_i{\sum }_{}^{}{y}_i}{n{\sum }_{}^{}{x_i}^2-({\sum }_{}^{}{x}_i{)}^2}\text{\hspace{1em}(5b)}$$
$$b=\frac{{\sum }_{}^{}{x_i}^2{\sum }_{}^{}{y}_i-{\sum }_{}^{}{x}_i{\sum }_{}^{}{x}_i{y}_i}{n{\sum }_{}^{}{x_i}^2-({\sum }_{}^{}{x}_i{)}^2}\text{\hspace{1em}(5b)}$$
已经求出了平方损失函数最小时对应的 $w$、 $b$ 参数值，这也就是最佳拟合直线。

def w_calculator(x, y):
    n = len(x)
    w = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x*x) - sum(x)*sum(x))
    b = (sum(x*x)*sum(y) - sum(x)*sum(x*y))/(n*sum(x*x)-sum(x)*sum(x))
    return w,b
w_calculator(x,y)
# (0.7545842753077117, 41.33509168550616)

最小二乘法矩阵推导

一元线性函数的表达式为 $y(x,w)=w_0+w_{1}x$(原式子里的w设为$w_1$，b设为$w_0$)，表达成矩阵形式为：
$$\left[\begin{array}{c}1,x_1\\ 1,x_2\\ \cdots \\ 1,x_9\\ 1,x_{10}\end{array}\right]\left[\begin{array}{c}{w}_0\\ w_1\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ \cdots \\ y_9\\ y_{10}\end{array}\right]\Rightarrow \left[\begin{array}{c}1,56\\ 1,72\\ \cdots \\ 1,94\\ 1,74\end{array}\right]\left[\begin{array}{c}{w}_0\\ w_1\end{array}\right]=\left[\begin{array}{c}92\\ 102\\ \cdots \\ 105\\ 92\end{array}\right]$$
$$y(x,w)=XW\text{\hspace{1em}(6)}$$
$(6)$ 式中，$W$ 为 $\left[\begin{array}{c}{w}_0\\ w_1\end{array}\right]$，而 $X$ 则是 $[X_1,x]$$X_1=\left[\begin{array}{c}1\\ 1\\ \cdots \\ 1\\ 1,\end{array}\right]$,$x=\left[\begin{array}{c}{x}_1\\ x_2\\ \cdots \\ x_9\\ x_{10}\end{array}\right]$ )矩阵。然后，平方损失函数为：

$$f=\sum _{i=1}^n{(y_i-(w_0+w_1{x}_i))}^2=(y-XW{)}^T(y-XW)\text{\hspace{1em}(7)}$$
$$f=y^{T}y-y^T(XW)-(XW{)}^{T}y+(XW{)}^T(XW)\text{\hspace{1em}(8)}$$
在该公式中 $y$ 与 $XW$ 皆为相同形式的 $(m,1)$ 矩阵，由此两者相乘属于线性关系，所以等价转换如下：
$$\begin{gathered} f=y^{T}y-(XW{)}^{T}y-(XW{)}^{T}y+(XW{)}^T(XW) \\ =y^{T}y-2(XW{)}^{T}y+(XW{)}^T(XW)\text{\hspace{1em}(9)} \end{gathered}$$
$W^T的偏导数是W，Tr(AB)对A或B的偏导数是B或A$
对矩阵求偏导得：
$$\frac{\partial f}{\partial W}=2X^{T}XW-2X^{T}y=0\text{\hspace{1em}(10)}$$

当矩阵 $X^{T}X$ 满秩时， $(X^{T}X{)}^{-1}{X}^{T}X=E$，且 $EW=W$。所以有 $(X^{T}X{)}^{-1}{X}^{T}XW=(X^{T}X{)}^{-1}{X}^{T}y$，并最终得到：
$$W=(X^{T}X{)}^{-1}{X}^{T}y\text{\hspace{1em}(11)}$$

def w_matrix(x, y):
    w = (x.T * x).I * x.T * y
    return w
x= [[1, i] for i in x]
x = np.matrix(x)
y = np.matrix(y)

w_matrix(x, y.reshape(-1, 1))
# matrix([[41.33509169],
#        [ 0.75458428]])

线性回归 Python 实现

以最小二乘法代数方法为例

w,b=w_calculator(x, y)
# 求损失
loss = mae_value(y,wx+b)

在上述例子中，得到的loss为447.69153479025357
接下来，我们尝试将拟合得到的直线绘制到原图中：

x_temp = np.linspace(50, 120, 100)  # 绘制直线生成的临时点

plt.scatter(x, y)
plt.plot(x_temp, x_temp*w + b, 'r')

此时，想要预估x为100对应的y只需要代入公式$y=wx+b$即可得到：116.79351921627732

线性回归 scikit-learn 实现

scikit-learn 把线性回归的过程整合到了LinearRegression() 类里，只需要填入需要的参数即可。

sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)
- fit_intercept: 默认为 True，计算截距项。
- normalize: 默认为 False，不针对数据进行标准化处理。
- copy_X: 默认为 True，即使用数据的副本进行操作，防止影响原数据。
- n_jobs: 计算时的作业数量。默认为 1，若为 -1 则使用全部 CPU 参与运算。

from sklearn.linear_model import LinearRegression

# 定义线性回归模型
model = LinearRegression()
model.fit(x.reshape(len(x), 1), y)  # 训练, reshape 操作把数据处理成 fit 能接受的形状

# 得到模型拟合参数
model.intercept_, model.coef_

得到的结果(41.33509168550615, array([0.75458428]))和上述python得到结果一致。

# 想要预测x=100对应的y
model.predict([[100]])

梯度下降法

为了求解$w$的极小值还可以引入一种叫「梯度下降」的求解方法。梯度下降法是一种十分常用且经典的最优化算法，通过这种方法我们就能快速找到函数的最小值。

梯度下降法的原理

什么是「梯度」？梯度是一个向量，它表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。对于一元函数而言，梯度就是指在某一点的导数。而对于多元函数而言，梯度就是指在某一点的偏导数组成的向量。

函数在沿梯度方向变化最快，所以「梯度下降法」的核心就是，我们沿着梯度下降方向去寻找损失函数的极小值（梯度的反方向）。过程如下图所示。

这里的损失函数依然是均方误差MSE：
$$J=\frac{f}{n}=\frac{1}{n}\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2\text{\hspace{1em}(1)}$$
$$J=\frac{1}{n}\sum _{i=1}^n(y_i-(w{x}_i+b){)}^2\text{\hspace{1em}(2)}$$
求解参数和截距项对应的梯度：
$$\frac{\partial J}{\partial w}=-\frac{2}{n}\sum _{i=1}^n{x}_i(y_i-(w{x}_i+b))\text{\hspace{1em}(3a)}$$

$$\frac{\partial J}{\partial b}=-\frac{2}{n}\sum _{i=1}^n(y_i-(w{x}_i+b))\text{\hspace{1em}(3b)}$$
当我们得到梯度的方向，然后乘以一个常数$\alpha$ ，就可以得到每次梯度下降的步长(上图箭头的长度)。最后，通过多次迭代，找到梯度变化很小的点，也就对应着损失函数的极小值了。其中，常数$\alpha$往往也被称之为学习率 Learning Rate。

在下文中用lr(Learning Rate)代表常数$\alpha$。

每次迭代，用w和b的初始值去减掉梯度*lr
$$w=w-\frac{\partial J}{\partial w}\ast lr\text{\hspace{1em}(4a)}$$
$$b=b-\frac{\partial J}{\partial b}\ast lr\text{\hspace{1em}(4b)}$$

梯度下降法求解线性回归

这里的x、y为数据点。

import numpy as np
import pandas as pd
def mse_value(y_true, y_pred):
    n = len(y_true)
    mse = sum(np.square(y_true - y_pred))/n
    return mse

def gradient_w(X,y,z):
    # 梯度计算
    gradient = 2*np.dot(X.T, (z- y)) / y.shape[0]
    return gradient

def gradient_b(X,y,z):
    # 梯度计算
    one = np.ones(y.shape[0])
    gradient = 2*np.dot(one, (z- y)) / y.shape[0]
    return gradient

def gradient_descent():

    w = 0  # 初始参数为 0
    b = 0  # 初始参数为 0
    
    lr = 0.000000001 # 设置合理学习率
    num_iter = 300 # 设置合理迭代次数
    l_list = []  # 保存损失函数值
    for i in range(num_iter):  # 梯度下降迭代
        z=x*w+b
        w=w-gradient_w(x,y,z)*lr
        b=b-gradient_b(x,y,z)*lr
        l = mse_value(y,z)  # 计算损失函数值
        l_list.append(l)
    return w, b
#w, b,l_y=gradient_descent()

绘制损失函数如下，发现曲线在迭代次数为300时趋于平缓，得合理迭代次数为300