这篇笔记主要介绍线性回归算法，对在第一篇笔记中介绍过的线性回归算法进行实现。kNN算法主要解决的是分类问题，并且它的结果不具备良好的解释性。线性回归算法主要解决回归问题，它的结果具有良好的可解释性，和kNN算法的介绍过程一样，线性回归算法也蕴含了机器学习中很多重要的思想，并且它是许多强大的非线性模型的基础。
简单线性回归
在第一篇笔记中，我们举过房屋面积大小和价格例子，将其绘制在二维坐标图上，横轴表示房屋面积，纵轴表示房屋价格，这里样本特征数据只有一个，那就是房屋面积，而在kNN算法的分类问题中，二维坐标图上横纵轴表示的都是样本特征数据，这是一个比较明显的区别。如果线性回归问题要在图中表示两种样本特征数据的话就需要三维空间坐标来表示。我们一般将只有一种样本特征数据的线性回归问题称为简单线性回归问题。
回顾线性回归
线性回归其实就是寻找一条直线，最大程度的拟合样本特征和样本输出标记之间的关系。
我们来看上面这张图，大家知道直线的方程是$y=ax+b$，那么点A肯定是在一条直线上，该条直线方程为$y^{(i)}=ax^{(i)}+b$，那么点A的横轴值为$x^{(i)}$，既是A的样本特征值，纵轴值为$y^{(i)}$，既是A的样本输出值。我们假设图中的红线就是拟合直线，方程为$\hat y^{(i)}=ax^{(i)}+b$，也就是将$x^{(i)}$代入这条红线，会得到一个预测的纵轴值$\hat y^{(i)}$。我们希望真值$y^{(i)}$和预测值$\hat y^{(i)}$的差值越小，说明我们的拟合直线拟合的越好。
因为差值有正有负，为了保证都是正数，所以将差值进行平方，之所以不用绝对值，是为了方便求导数。
方程求导的知识可参阅 《机器学习笔记一之机器学习定义、导数、最小二乘》 。
$$ (y^{(i)} - \hat y^{(i)})^2 $$
将所有样本特征都考虑到，既将所有真值和预测值的差值求和：
$$ \sum_{i=1}^m(y^{(i)} - \hat y^{(i)})^2 $$
将$ax^{(i)}+b$代入上面的公式就得到：
$$\sum_{i=1}^m(y^{(i)} - ax^{(i)}-b)^2$$
上面的公式我们称为损失函数（Loss Function），损失函数值越小，我们的拟合直线越好。在Loss函数中，$a$和$b$是变量，所以我们要做的就是找到使Loss函数值最小的$a$和$b$。这个套路是近乎所有参数学习算法常用的套路，既通过分析问题，确定问题的损失函数，通过最优化损失函数获得机器学习的模型。像线性回归、多项式回归、逻辑回归、SVM、神经网络等都是这个套路。
上述的Loss函数是一个典型的最小二乘法的问题，既通过最小化误差的平方和寻找数据的最佳函数匹配。求函数的最小值就会用到导数这个数学工具，具体如何推导上面的Loss函数可以参见第一篇学习笔记，这里不再累赘。最后得出a和b的求解公式为：
$$a=\frac {\overline x \ \overline y-\overline {xy}} {(\overline x)^2-\overline {x^2}}$$
$$b=\overline y-a\overline x$$
实现简单线性回归
我们先在Jupyter Notebook中实现：
import numpy as np |
然后我们使用上面推导出的公式求出$a$和$b$：
a = (np.mean(x)*np.mean(y) - np.mean(x*y))/(np.mean(x)**2 - np.mean(x**2)) |
# 新来一个特征值，利用拟合直线计算输出值 |
封装简单线性回归方法
我们在PyCharm中封装我们自己的简单线性回归方法：
import numpy as np |
然后我们就可以在Jupyter Notebook中使用我们封装的简单线性回归方法：
from myML.SimpleLinearRegression import SimpleLinearRegression |
线性回归算法的评测标准
在讲kNN算法时，我们分类问题的评测标准是基于将样本数据拆分为训练数据和测试数据的前提下的，那么在线性回归算法中也是一样的。我们使用训练数据计算出a和b的值，然后将测试数据代入拟合直线方程算出结果，进行比较。我们通过公式来看一下。
$$ \sum_{i=1}^m(y_{train}^{(i)} - ax_{train}^{(i)} -b)^2 = \sum_{i=1}^m(y_{train}^{(i)} - \hat y_{train}^{(i)})^2 $$
将训练数据代入公式算出a和b，然后代入拟合直线方程算出$\hat y_{test}^{(i)}$：
$$ \hat y_{test}^{(i)} = ax_{test}^{(i)} + b $$
此时我们的衡量标准既为：
$$ \sum_{i=1}^m(y_{test}^{(i)} - \hat y_{test}^{(i)})^2 $$
也就是上面的公式值越小说明我们拟合的越好。
均方误差（MSE）
上面这个公式有一个问题，那就是最终值受$m$的影响，比如某个算法10个样本数据求出的值为80，另一个算法10000个样本数据求出的值为100，也不能表明第一个算法就比第二个算法好，因为样本数据量相差巨大，所以将上面公式改变一下，将值除以$m$：
$$ \frac 1 m \sum_{i=1}^m(y_{test}^{(i)} - \hat y_{test}^{(i)})^2$$
这个衡量标准称为均方误差（MSE, Mean Squared Error）
均方根误差（RMSE）
在对量纲不敏感的情况下，使用均方误差没什么问题，但是在一些对量纲比较敏感的场景下，均方误差就会有问题，因为均方误差的量纲为XX平方，比如房屋面积售价的例子，均方误差的量纲就成了$元^2$，所以就有了均方根误差（RMSE, Root Mean Squared Error）用以统一量纲：
$$ \sqrt {\frac 1 m \sum_{i=1}^m(y_{test}^{(i)} - \hat y_{test}^{(i)})^2}= \sqrt {MSE_{test}} $$
平均绝对误差（MAE）
另外一个能统一量纲的衡量公式为平均绝对误差，既将真实值与预测值的差取绝对值而不是平方：
$$ \frac 1 m \sum _{i=1}^m |y_{test}^{(i)} - \hat y_{test}^{(i)}| $$
实现MSE， RMSE， MAE
实现这三个指标我们使用Scikit Learn中提供的波士顿房价的数据集，我们先使用简单线性回归进行预测：
import numpy as np |
从图中我们可以看到在顶部有一些似乎到最大值的数据，这是因为在真实的数据中有一些类似50万以上这类数据，都会被归为数据集最大值一类，这些数据对我们的预测不但没有帮助，反而会有影响，所以我们将这些数据去掉：
x = x[y < 50] |
下面我们来计算MSE：
y_test_predict = slr.predict(x_test) |
再来计算RMSE：
rmse_test = np.sqrt(mse_test) |
再来看看MAE：
mae_test = np.sum(np.absolute(y_test - y_test_predict)) / len(y_test) |
从结果可以看出RMSE比MAE的值要大，那是因为RMSE是差值先平方然后求和然后再开方，所以如果有某个真实值和预测值之间差距比较大的时候，平方操作就会放大数据的量级。所以一般我们使用RMSE更有实际意义，因为RMSE的值小，说明了最大误差比较小。
封装MSE， RMSE， MAE
我们将这三个衡量指标封装起来，在metrics.py文件中增加三个方法：
import numpy as np |
这样就可以在Jupyter Notebook中方便的使用了：
from myML.metrics import mean_squared_error |
R Squared
在之前的分类问题中，衡量标准的值都是在0/1之间，1表示最好，0表示最差，这个值和量纲无关。但是在线性回归问题中衡量标准的值是带有量纲的，比如某个算法在预测房价的场景中RMSE是5万，但在预测学生分数的场景中RMSE是10分，那么这个算法是在预测房价场景中好呢还是在预测学生分数场景中好呢？这个是无法判断的，这就是RMSE和MAE的局限性。那么为解决这个问题，就出现了一个新的衡量指标R Squared，也就是$R^2$。这个指标也是目前机器学习算法使用比较广泛的一个指标。我们先来看看$R^2$的公式：
$$ R^2 = 1 - \frac {\sum _{i=1}^m(\hat y^{(i)} - y^{(i)})^2} {\sum _{i=1}^m(\bar y - y^{(i)})^2} $$
这个公式的分子其实就是简单线性回归的模型预测产生的错误，既Loss函数。分母中的$\bar y$是均值，其实均值也是一种线性模型，只不过是比较粗糙的线性模型，所以分母是使用均值模型预测产生的错误。那么如果这个比值远远小于1，说明我们的模型的质量远远超出均值模型，那么$R^2$就无限接近于1，说明我们的模型拟合的非常好。如果比值接近1，说明我们的模型和均值模型没差多少，表明我们的模型比较烂，此时$R^2$就会接近0。那么综上，$R^2$的值小于等于1，值越大表示模型越好，当值为负数时表明我们的模型还不如均值模型，也表明了我们分析的数据之间可能根本就没有线性关系。
我们对$R^2$的公式再处理一下，将1后面的分数分子分母各除以$m$：
$$ R^2 = 1 - \frac {\sum _{i=1}^m(\hat y^{(i)} - y^{(i)})^2} {\sum _{i=1}^m(\bar y - y^{(i)})^2} = 1 - \frac {(\sum _{i=1}^m(\hat y^{(i)} - y^{(i)})^2)/m} {(\sum _{i=1}^m(\bar y - y^{(i)})^2)/m} $$
此时，分子其实就是上文中讲到的MSE，而分母就是我们在第二篇笔记中讲过的方差，所以$R^2$又可以写为：
$$ R^2 = 1 - \frac {MSE(\hat y, y)} {Var(y)} $$
实现R Squared
实现$R^2$其实是非常简单的：
r_square = 1 - mean_squared_error(y_test, y_test_predict)/np.var(y_test) |
我们同样将$R^2$的方法封装起来：
import numpy as np |
这样就可以在Jupyter Notebook中方便的使用了：
from myML.metrics import r2_score |
Scikit Learn中的$R^2$用法也是一样的:
from sklearn.metrics import r2_score |
多元线性回归
上面我们讲了简单线性回归，也就是只关注样本数据的一个特征，这一节我们来看看多元线性回归，既关注样本数据的多个特征。
像上图中展示的，简单线性回归我们只关注$x^{(i)}$一个特征，假如$x^{(i)}$是一个特征向量，那此时就变成了关注多个特征的多元线性回归问题，如下图所示：
此时相对与简单线性回归的公式$y=ax+b$，因为$x$成为了特征行向量，所以我们将$a$也看作特征系数列向量，那么公式展开后可以写成这样：
$$ y = a_1x_1+a_2x_2+…+a_nx_n+b $$
我们按惯例，将多元线性回归公式中的特征系数称为$\theta$：
$$ y = \theta_1x_1+\theta_2x_2+…+\theta_nx_n+\theta _0 $$
那么描述每行样本数据中每个特征和其预测值的公式为：
$$ \hat y^{(i)} = \theta _0+\theta_1X_1^{(i)}+\theta_2X_2^{(i)}+…+\theta_nX_n^{(i)} $$
我们对上面的公式再进一步做一下处理，将截距$\theta_ 0$也乘以一个特征$X_0^{(i)}$，不过该特征值恒等于1:
$$ \hat y^{(i)} = \theta _0X_0^{(i)}+\theta_1X_1^{(i)}+\theta_2X_2^{(i)}+…+\theta_nX_n^{(i)} $$
我们将上面公式用矩阵的方式表示为（用到矩阵相乘的知识，1行n列矩阵乘以n行1列矩阵为1行1列矩阵，既一个标量）：
$$ \hat y^{(i)} = \begin{bmatrix} X_0^{(i)}& X_1^{(i)}&X_2^{(i)}& … &X_n^{(i)}\end{bmatrix}\begin{bmatrix}
\theta_ 0 \\
\theta_ 1 \\
\theta_ 2 \\
… \\
\theta_ n \\
\end{bmatrix} $$
我们可以直接将上面公式表示为：
$$ \hat y = X_b \theta $$
我们知道简单线性回归就是使$\sum_{i=1}^m(y^{(i)}-\hat y^{(i)})^2$尽可能小，那多元线性回归也是一样的，只不过多元线性回归中的$\hat y^{(i)}$的公式不一样而已，将上面的公式代入：
$$\sum_{i=1}^m(y-X_b \theta)^2$$
将上面公式里的平方展开：
$$\sum_{i=1}^m(y-X_b \theta)(y-X_b \theta)$$
矩阵的相乘求和可以转换为如下形式（矩阵或向量的相同元素相乘再求和等于该向量或矩阵的转置乘以该向量或矩阵）：
$$(y-X_b \theta)^T(y-X_b \theta)$$
所以多元线性回归就是求使上面公式尽可能小的$\theta$列向量。
上面的公式出现了矩阵乘积转置，看看下面这张手记就能明白如何转换矩阵乘积转置了：
如此，上面的公式可以转换为：
$$ F_{loss}=(y^T-\theta^T X_b^T)(y-X_b \theta) = y^Ty-y^TX_b\theta-\theta^TX_b^Ty+\theta^TX_b^TX_b\theta$$
要求$F_{loss}$函数的最小值既对该函数中的每项$theta$求偏导数：
$$\frac {\partial F_{loss}}{\partial \theta}=\frac {\partial y^Ty}{\partial \theta}-\frac {\partial y^TX_b\theta}{\partial \theta}-\frac {\partial \theta^TX_b^Ty}{\partial \theta}+\frac {\partial \theta^TX_b^TX_b\theta}{\partial \theta}$$
这里用到了矩阵求导的知识，具体可见下面的手记：
上面两张手记简要介绍了矩阵求导的一种方法，根据求导的类型，然后确定是分子布局还是分母布局，最后查阅常用的矩阵求导表确定结果。
- $\frac {\partial y^Ty}{\partial \theta}$：分子是标量，分母是列向量，是标量/向量类型，并且是分母布局，查表满足$\frac{\partial a}{\partial \mathbf{x}}$条件，故可得结果为0。
- $\frac {\partial y^TX_b\theta}{\partial \theta}$：分子是标量，分母是列向量，是标量/向量类型，并且是分母布局，查表满足${\frac {\partial {\mathbf {b}}^{\top }{\mathbf {A}}{\mathbf {x}}}{\partial {\mathbf {x}}}}$条件，故可得结果为$X_b^Ty$。
- $\frac {\partial \theta^TX_b^Ty}{\partial \theta}$：分子是矩阵，分母是列向量，是矩阵/向量类型，并且是分母布局，但是根据矩阵乘积转置的规则，可以将分子转换为$\frac {\partial y^TX_b\theta}{\partial \theta}$，和第二项偏导表达式相同，故结果为$X_b^Ty$。
- $\frac {\partial \theta^TX_b^TX_b\theta}{\partial \theta}$：分子为标量，分布为列向量，是标量/向量类型，并且是分母布局，因为$X_b^TX_b$还是矩阵，所以查表满足${\frac {\partial {\mathbf {x}}^{\top }{\mathbf {A}}{\mathbf {x}}}{\partial {\mathbf {x}}}}$条件，故可得结果为$2X_b^TX_b\theta$。
所以$F_{loss}$函数的最小值为：
$$F(min)_{loss}=\frac {\partial F_{loss}}{\partial \theta}=\frac {\partial y^Ty}{\partial \theta}-\frac {\partial y^TX_b\theta}{\partial \theta}-\frac {\partial \theta^TX_b^Ty}{\partial \theta}+\frac {\partial \theta^TX_b^TX_b\theta}{\partial \theta}=0-2X_b^Ty+2X_b^TX_b\theta$$
便可求得$\theta$的值为：
$$-2X_b^Ty+2X_b^TX_b\theta=0$$
$$\theta=\frac {X_b^Ty}{X_b^TX_b}$$
根据逆矩阵的转换规则可得：
$$\theta=(X_b^TX_b)^{-1}X_b^Ty$$
上面的公式就是多元线性回归的正规方程解 (Normal Equation)。
使用正规方程解实现多元线性回归
在实现之前，我们先明确一下$\theta$，它是一个列向量：
$$\theta=\begin{bmatrix}
\theta_ 0\
\theta_ 1\
\theta_ 2\
…\
\theta_ n\
\end{bmatrix}$$
其中$\theta_0$是多元线性方程的截距（intercept），$\theta_1$到$\theta_n$才是系数（coefficients）。
在PyCharm中，我们创建一个类LinearRegression，然后构造函数如下：
class LinearRegression: |
然后我们来看训练的过程，既fit方法，这里因为是用正规方程解实现的，所以我们将训练方法的名称取为fit_normal()：
def fit_normal(self, X_train, y_train): |
训练的实现过程很简单，首先求出X_b，也就是将X_train矩阵前面填加元素为1的一列，这里我们使用到了Numpy的np.hstack方法，也就是在水平方向组合矩阵，另外使用了np.ones这个快捷创建元素为1的矩阵的方法。然后实现了上文中推导出的正规方程解，其中np.linalg.inv是对矩阵求逆的方法。
然后来看预测的过程：
# 给定待预测数据集X_predict，返回表示X_predict的结果向量 |
首先有两个断言增加健壮性，然后同样是计算出X_b，最后乘以训练出的$\theta$就得到了预测的结果。
最后同样增加一个评分的方法，使用我们之前实现的$R^2$方法来计算评分：
# 根据测试数据集X_test和y_test确定当前模型的准确度 |
在Jupyter Notebook中使用LinearRegression
我们使用Scikit Learn中提供的波士顿房价数据来验证我们实现的基于正规方程解的多元线性回归。
import numpy as np |
我们再来看看Scikit Learn中提供的线性回归的使用方式：
from sklearn.linear_model import LinearRegression |
我们可以看到使用Scikit Learn中的线性回归和我们实现的线性回归结果是一样的，但这里要注意的是Scikit Learn中的线性回归实现方式并不是正规方程解，只是因为样本数据量比较少，所以得到了相同的结果。下一篇笔记会介绍另外一种实现多元线性回归的方式。
线性回归的可解释性
我们将全量的波士顿房价进行预测，求出系数来看一下：
from sklearn.linear_model import LinearRegression |
我们看到系数有正有负，表示了正相关和负相关，既正系数越大的特征对结果的影响越大，负系数越小对结果影响越大，我们来看看影响波士顿房价的特征都是什么，先将系数按索引从小到大排序：
np.argsort(lin_reg.coef_) |
通过boston.DESCR我们可以知道正系数最大的特征是RM，既房屋面积，第二大的正系数特征是CHAS，既房屋临Charles河的距离，离河越近房屋价格越高，从这两项都可以看出这是合理的情况。再来看看最小的系数特征NOX，它的系数是负的，这个特征是房屋周围的一氧化碳浓度，浓度约小房价越高，可见这都是合理的。综上，就解释了线性回归的可解释性，它可以让我们更好的采集收集样本数据，提高数据质量，提升预测准确性。
总结
这篇笔记介绍了线性回归和多元线性回归的概念、公式推导、评价标准以及实现。是解决一些预测类问题的基本知识。下篇笔记将学习机器学习算法中重要的一个算法梯度下降。
申明：本文为慕课网liuyubobobo老师《Python3入门机器学习 经典算法与应用》课程的学习笔记，未经允许不得转载。