第2章从零认识线性回归（附高清PDF和PPT）

经过第1章的介绍，我们已经完成了对于Python开发环境的安装与配置。现在，就正式开始介绍第一个算法：线性回归（Linear Regression）。整个线性回归的学习路线如图2-1所示。由于这是介绍的第一个算法，所以我们会介绍很多基本的内容，导致看起来有很多内容。因此，对于整个线性回归的学习将会通过7个小节的内容来进行介绍。但是，值得高兴的是只要完成了前4个步骤的学习就基本上是达到了第一阶段的要求。下面就将正式开始介绍线性回归。

p1-1

图 2-1 线性回归学习路线图

2.1模型的建立与求解

2.1.1理解线性回归模型

通常来说，机器学习中的每一个算法都是为了解决某一类问题而诞生。换句话说，也就是在实际情况中存在一些问题能够通过线性回归来解决，例如对房价的预测。但是有人可能会问，为什么对于房价的预测就应该用线性回归，而不是其它算法呢？其原因就在于常识告诉我们房价都是随着面积的增长而增长，且总体上呈线性增长的趋势。那有没有那种当面积大到一定程度后价格反而降低，因此不符合线性增长的呢？这当然也可能存在，但在实际处理中肯定会优先选择线性回归模型，当效果不佳时我们会再尝试其它算法。因此，当学习过多个算法后，在拿到某个具体的问题时，就需要考虑哪种模型更适合来解决这个问题了。

图 2-2 房价走势图

某市房价的一个走势图如图2-2所示，其中横坐标为面积，纵坐标为价格，且房价整体上呈线性增长的趋势。那假如现在随意告诉你一个房屋的面积，要怎么才能预测（或者叫计算）出其对应的价格呢？

2.1.2建立线性回归模型

一般来说，当我们拿到一个实际问题时，首先会根据问题的背景结合常识选择一个合适的模型。同时，现在常识告诉我们房价的增长更优先符合线性回归这类模型，因此可以考虑建立一个如下的线性回归模型

\hat{y}=h(x)=wx+b\;\;\;\;\;\;(2.1)

$w$ $b$ $w,b$ $x$ ，就能够预测出其对应的房价。

注意：在机器学习中所谓的模型，可以简单理解为一个函数。

2.1.3求解线性回归模型

$h(x)$ 进行求解。在中学时期我们倒是学过如何通过两个坐标点来求解过这两点的直线，可在上述的场景中这种做法显然行不通的（因为所有的点并不在一条直线上），那有没有什么好的解决的办法呢？

$h(x)$ $h(x)$ $h(x)$ $x$ $y$ $h(x)$ 不好入手，那么我们就从“准确”来入手。

可又怎么来定义准确呢？在这里，我们可以通过计算每个样本的真实房价与预测房价之间的均方误差来对“准确”进行刻画

\begin{aligned} & J(w,b)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}} \\ & {{{\hat{y}}}^{(i)}}=h({{x}^{(i)}})=w{{x}^{(i)}}+b \end{aligned} \;\;\;\;\;\;(2.2)

$m$ $x^{(i)}$ $i$ $i$ $y^{(i)}$ $i$ $\hat{y}^{(i)}$ $i$ 个房屋的预测价格。

$J(w,b)$ $\hat{w},\hat{b}$ $J(w,b)$ $J(w,b)$ 就是0了。

$h(x)$ $J(w,b)$ $J(w,b)$ 也有一个专门的术语叫目标函数（Objective Function）或者是代价函数（Cost Function）亦或是损失函数。

至此，我们离第一阶段的学习目标就只差一步之遥，那就是如何通过开源框架来进行建模求解，并进行预测。至于求解过程到底怎样如何进行的，那就是第二阶段的任务了，下面开始完成第一阶段的最后一步。

2.1.4 sklearn简介

如图2-3所示，scikit-learn简称sklearn[1]。它是一个开源的机器学习框架，常用的机器学习算法都可以在里面找到，例如线性回归、逻辑回归、决策树等等。同时，其Python化的设计风格对于Python用户来说也十分友好与易于上手，并且每个算法都给了相应的示例以及API文档。

图 2-3 sklearn介绍图

2.1.5安装sklearn以及其它库

在第1章我们已经介绍了如何配置Python环境，下面就开始在第1章中所新建的Python虚拟环境里安装所需要的包（库）。

1) 打开终端

如果是在Windows环境中，则先点击“开始”，然后找到Anaconda Prompt并打开，最后再激活相应的Python虚拟环境即可。如果是在Linux环境中，则直接打开命令行终端，然后再激活相应的虚拟环境。

2) 安装Python包

为了完成整个线性回归的建模任务以及结果的可视化，需要安装sklearn和matplotlib这两个Python包。如果之前在配置环境的时候，替换了pip源地址为清华镜像，则下载的速度会更快，如图2-4所示。

图 2-4 sklearn安装示意图

由于在安装相应的Python包后PyCharm会建立新的包索引，所以需要等待几分钟，如图2-5所示。

图 2-5 Python包索引图

2.1.6线性回归示例代码

在安装完成相应的Python包后就可以正式的对线性回归进行建模与求解，详细完整代码见Chapter02/01_house_price_train.py文件。

1) 导入包

首先需要将用到的相关Python包进行导入，代码如下：


3
1
import numpy as np
2
from sklearn.linear_model import LinearRegression
3
import matplotlib.pyplot as plt

2) 制作样本（数据集）

制作用于训练模型的数据集，代码如下：


6
1
def make_data():
2
    np.random.seed(20)
3
    x = np.random.rand(100) * 30 + 50  # 面积
4
    noise = np.random.rand(100) * 50
5
    y = x * 8 - 127 - noise              # 价格
6
    return x,y

上述代码中第2行代码的作用是在真实的房价中加入一定的噪音（误差）。

3) 定义模型并求解

通过LinearRegression模型来对模型的参数进行求解与预测，代码如下：


6
1
def main(x, y):
2
    model = LinearRegression()            # 定义模型
3
    x = np.reshape(x, (-1, 1))
4
    model.fit(np.reshape(x, (-1, 1)), y)# 求解模型
5
    y_pre = model.predict(x) # 预测
6
    print(y_pre)

上述代码中np.reshape(x,(-1,1))表示把x变成[n,1]的形状，至于n到底是多少，将通过np.reshape函数自己推导得出。例如x的shape为[4,5]，如果想把a改成[2,10]形状则可以使用a.reshape([2,10])，或者a.reshape([2,-1])来进行形状的变换。

3) 运行结果

最后，调用定义好的函数运行程序，并输出最后训练得到的参数结果，代码如下：


5
1
if __name__ == '__main__':
2
    x, y = make_data()
3
    main(x, y)      
4
    #参数w=[7.97699647],b=-154.31885006061555
5
    #面积50的房价为： [244.53097351]

$w=7.97,b=-154.32$ $h(x)=7.97x-154.32$ $h(x)$ 来对新的输入进行预测。同时，还能够根据求解后的模型画出对应拟合出的直线，如图2-6所示。

图 2-6 房价预测走势图

到此，便完成了对于线性回归第一阶段的大致学习。

2.1.7小结

在本节内容中，笔者首先通过一个实际的场景介绍了什么是线性回归；接着介绍了如何建立一个简单的线性模型；然后引导大家如何将模型的求解问题，转化为目标函数的最小化过程；最后通过开源框架sklearn搭建了一个简单的线性回归模型并进行了求解。虽然内容不多，但是却体现了整个线性回归算法的核心思想。

2.2 多变量线性回归

2.2.1 理解多变量

在2.1节的内容中，笔者详细的介绍了什么是线性回归以及一个典型的应用场景；同时还介绍了如何通过开源的sklearn来搭建一个简单的线性回归模型，使得对于线性回归的核心思想有了一定的掌握；接下来，我们将继续开始学习线性回归后续的内容。

在这里还是以房价预测为例。尽管影响房价的主要因素是面积，但是其它因素同样也可能影响到房屋的价格。例如房屋到学校的距离、到医院的距离和到大型商场的距离等等（总不能卖你一套深山老林的房子你也要吧）。虽然现实生活中没有这么量化，但是开发商也总是会拿学区房做卖点。所以，这时便有了影响房价的四个因素，而在机器学习中我们将其称之为特征（Feature）。因此，包含有多个特征的线性回归就叫做多变量线性回归（Multiple Linear Regression）。

2.2.2多变量线性回归建模

以波士顿房价数据集为例，其一共包含了13个特征属性。因此可以得到如下线性模型

h(x)={{w}_{1}}{{x}_{1}}+\cdots +{{w}_{13}}{{x}_{13}}+b\;\;\;\;\;\;(2.3)

且同时，其目标函数为

\begin{aligned} & J(W,b)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}} \\ & {{{\hat{y}}}^{(i)}}=h({{x}^{(i)}})={{w}_{1}}x_{1}^{(i)}+\cdots +{{w}_{13}}x_{13}^{(i)}+b \end{aligned}\;\;\;\;\;\;(2.4)

$x^{(i)}_j$ $i$ $j$ $W$ $b$ 为一个标量表示偏置。

$J(W,b)$ 后，便可以求解出模型对应的参数。

2.2.3 多变量回归示例代码

下面依旧以sklearn来进行多变量线性回归模型的建模与求解，完整代码见Chapter02/02_boston_price_train.py文件。

1) 导入数据集

这里直接导入了一个sklearn内置的Boston房价数据集为例进行演示，代码如下：


5
1
def load_data():
2
    data = load_boston()
3
    x = data.data
4
    y = data.target
5
    return x, y

2) 求解与结果

训练模型与输出相应的权重参数和预测值，代码如下：


8
1
def train(x, y):
2
    model = LinearRegression()
3
    model.fit(x,y)
4
    print("权重为：",model.coef_,"偏置为：",model.intercept_)
5
    print("第12个房屋的预测和真实价格：",model.predict(x[12,:].reshape(1,-1)))
6
    
7
#参数w=[7.97699647],b=-154.31885006061555
8
#面积50的房价为： [244.53097351]

根据上述代码便完成了对于多变量线性回归模型的建立与求解，同时也得出了各个特征所对应的权重参数。但由于不易对高维数据进行可视化，所以这里只能从预测的结果来评判模型的好坏，具体的模型评估指标将在第2.4节内容中进行介绍。

2.3 多项式回归

2.3.1 理解多项式

假定已知矩形的面积公式，而不知道求解梯形的面积公式，且同时你手上有若干个类似图2-7所示的梯形图。已知梯形的上底和下底，并且上底均等于高。现在让你建立一个模型，当给定一个上述梯形时你能近似的给出其面积。这该如何进行建模呢？

图 2-7 梯形图

2.3.2 多项式回归建模

$h(x)=w_1x_1+w_2x_2+b$ $(x_1)^2$ $x_1x_2$ $(x_2)^2$ $(2.5)$ 所示的模型

h(x)=x_1w_1+x_2w_2+(x_1)^2w_3+x_1x_2w_4+(x_2)^2w_5+b\;\;\;\;\;\;(2.5)

$w_i$ $h(x)$ $x_1x_2,(x_1)^2,(x_2)^2$ 这些项（其次数高于1），因此将其称之为多项式回归（Polynomial Regression）。

但是，只要进行如下替换，便又回到了普通的线性回归：

h(x)={{x}_{1}}{{w}_{1}}+{{x}_{2}}{{w}_{2}}+{{x}_{3}}{{w}_{3}}+{{x}_{4}}{{w}_{4}}+{{x}_{4}}{{w}_{5}}+b\;\;\;\;\;\;(2.6)

$x_3=(x_1)^2,x_4=x_1x_2,x_5=(x_2)^2$ $h(x)$ 的也将是包含5个特征的数据。

2.3.3 多项式回归示例代码

要完成整个多项式回归的建模，首先要对原始的特征进行转换，构造模型所需要的新特征；然后是构造相应的数据集；最后在进行模型的训练与预测。完整代码见Chapter02/03_trapezoid_polynomial_train.py文件。

1) 特征转化

这里首先介绍一下sklearn中的多项式变换模块PolynomialFeatures，它的作用便是对每个特征进行指定的幂次变换，代码如下：


5
1
from sklearn.preprocessing import PolynomialFeatures
2
a = [[3, 4], [2, 3]]  
3
model = PolynomialFeatures(degree=2,include_bias=False)
4
b = model.fit_transform(a)
5
print(b)#输出结果：[[ 3.  4.  9. 12. 16.] [ 2.  3.  4.  6.  9.]]

2) 构造数据集

构造训练模型时所需要用到的数据集，代码如下：


6
1
def make_data():
2
    np.random.seed(10)
3
    x1 = np.random.randint(5, 10, 50).reshape(50, 1)
4
    x2 = np.random.randint(10, 16, 50).reshape(50, 1)
5
    x,y = np.hstack((x1, x2)), 0.5 * (x1 + x2) * x1
6
    return x, y

根据上述代码便得到了一个50行2列的样本数据，其中第一列为上底，第二列为下底。np.hstack的作用是将两个50行1列的向量合并成一个50行2列的矩阵。

3) 建模求解与结果

首先对原始特征进行多项式构造处理，然后进行模型的训练，代码如下：


xxxxxxxxxx
7
1
def train(x, y):
2
    poly = PolynomialFeatures(degree=2, include_bias=False)
3
    x_mul = poly.fit_transform(x)
4
    model = LinearRegression()
5
    model.fit(x_mul, y)
6
    print("权重为：", model.coef_)#[[0.  0.  0.5   0.5   0]]
7
    print("偏置为：", model.intercept_) # [0.]

根据上述代码建模求解后，就能够预测得到上底为 5，下底为 8的梯形面积为32.5（真实面积也是32.2）。并且，根据求解得到的权重和偏置得

\begin{aligned} & h(x)={{x}_{1}}\cdot 0+{{x}_{2}}\cdot 0+{{x}_{3}}\cdot 0.5+{{x}_{4}}\cdot 0.5+{{x}_{5}}\cdot 0+b \\[1ex] & =0.5\cdot {{({{x}_{1}})}^{2}}+0.5\cdot {{x}_{1}}\cdot {{x}_{2}} \\[1ex] & =0.5\cdot {{x}_{1}}({{x}_{1}}+{{x}_{2}}) \end{aligned} \;\;\;\;\;(2.7)

可以发现，此时模型居然已经自己总结（学习）出了梯形的计算公式。

2.3.4 小结

在本节内容中，笔者首先以两个示例来分别介绍了多变量线性回归和多项式回归；然后通过sklearn来对模型进行了求解；最后和读者一起领略到了算法的魅力所在。在2.4节中，我们将开始对模型的评估进行介绍。

2.4 回归模型评估

在2.1到2.3这3节内容中，笔者介绍了如何建模线性回归（包括多变量与多项式回归）、如何通过sklearn搭建模型并求解。但是对于一个求解出来的模型应该怎样来对其进行评估呢？换句话说，这个模型到底怎么样？

$h_1(x)$ $h_2(x)$ ，那么应该选哪一个呢？亦或是在不能可视化的情况下，应该来如何评估模型的好与坏呢？

图2-8 不同模型房价的预测走势图

$n$ $y_i$ $i$ $\hat{y}_i$ $i$ 个样本的预测值。

2.4.1常见回归评估指标

1) 平均绝对误差（MAE）

MAE用来衡量预测值与真实值之间的平均绝对误差，定义如下：

\text{MAE}=\frac{1}{n}\sum\limits_{i=1}^{n}{|}{{y}_{i}}-{{\hat{y}}_{i}}|\;\;\;\;\;\;\;\;\;(2.8)

$\text{MAE}\in [0,+\infty )$ ，越小表示模型越好，实现代码如下所示：


xxxxxxxxxx
2
1
def MAE(y, y_pre):
2
    return np.mean(np.abs(y - y_pre))

2) 均方误差（MSE）

MSE用来衡量预测值与真实值之间的误差平方，定义如下：

\text{MSE}=\frac{1}{n}\sum\limits_{i=1}^{n}{{{({{y}_{i}}-{{{\hat{y}}}_{i}})}^{2}}}\;\;\;\;\;\;\;\;(2.9)

$\text{MSE}\in [0,+\infty )$ ，越小表示模型越好，实现代码如下所示：


xxxxxxxxxx
2
1
def MSE(y, y_pre):
2
    return np.mean((y - y_pre) ** 2)

3) 均方根误差（RMSE）

RMSE是在MSE的基础之上开根号而来，其定义如下：

\text{RMSE}=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}{{{({{y}_{i}}-{{{\hat{y}}}_{i}})}^{2}}}}\;\;\;\;\;\;\;\;(2.10)

$\text{RMSE}\in [0,+\infty )$ ，越小表示模型越好，实现代码如下所示：


xxxxxxxxxx
2
1
def RMSE(y, y_pre): 
2
    return np.sqrt(MSE(y, y_pre))

4) 平均绝对百分比误差（MAPE）

MAPE和MAE类似，只是在MAE的基础上做了标准化处理，其定义如下：

\text{MAPE}=\frac{100%}{n}\sum\limits_{i=1}^{n}{\left| \frac{{{y}_{i}}-{{{\hat{y}}}_{i}}}{{{y}_{i}}} \right|}\;\;\;\;\;\;\;\;(2.11)

$\text{MAPE}\in [0,+\infty )$ ，越小表示模型越好，实现代码如下所示：


xxxxxxxxxx
2
1
def MAPE(y, y_pre):
2
    return np.mean(np.abs((y - y_pre) / y))

$R^2$ 评价指标

$R^2$ 是线性回归模型中sklearn默认采用的评价指标，其定义如下：

{{R}^{2}}=1-\frac{\sum\limits_{i=1}^{n}{{{({{y}_{i}}-{{{\hat{y}}}_{i}})}^{2}}}}{\sum\limits_{i=1}^{n}{{{({{y}_{i}}-\bar{y})}^{2}}}}\;\;\;\;\;\;\;\;(2.12)

${{R}^{2}}\in (-\infty ,1]$ $\overline{y}$ 表示真实值的平均值，实现代码如下所示：


xxxxxxxxxx
4
1
def R2(y, y_pre):
2
    u = np.sum((y - y_pre) ** 2)
3
    v = np.sum((y - np.mean(y_pre)) ** 2)
4
    return 1 - (u / v)

2.4.2 回归指标示例代码

有了这些评估指标后，在对模型训练时就可以选择其中的一些指标对模型的精度进行评估。这里以前面波士顿房价的预测结果为例进行示例，完整代码见Chapter02/04_metrics_boston_price.py文件。


xxxxxxxxxx
4
1
def train(x, y):
2
    model = LinearRegression()
3
    model.fit(x, y)
4
    y_pre = model.predict(x)

2.4.3小结

在本节中，笔者详细的介绍了如何评价一个回归模型的好坏，以及一些常用的评估指标和实现方式。最后，笔者还通过波士顿房价预测示例来展示了评价指标的用法。到此，对于线性回归模型在整个阶段一部分的内容就介绍完了。在2.5节的内容中，笔者将继续介绍如何通过梯度下降算法来求解目标函数，以及这个目标函数的由来。

2.5梯度下降

$J(W,b)$ ，但并没有给出严格的数学定义。同时，在求解的过程中也是直接通过开源框架所实现的，也不知道其内部的真正原理。因此，在这一节内容中我们将会仔细地学习目标函数的求解过程以及最小二乘法。

$A$ $B$ 处呢?

图2-9 跳跃方向图

$Y$ $\bf{V1}$ $X$ $\it{\bf{V2}}$ $\bf{l}$ $\bf{l}$ 所在的方向往前跳第一步，然后接着再选类似的方向往前跳第二步直到谷底。可为什么都会这样选呢？答：这还用问一看就知，不信你自己看。

2.5.1方向导数与梯度

$f(x)$ $x_0$ $f(x)$ $x=x_0$ $|f^{\prime}(x_0)|$ $f(x)$ $\Delta x$ $\Delta y$ $z=f(x,y)$ $z$ $A$ $z$ 在该点的方向导数

\frac{\partial f}{\partial \mathbf{l}}=\{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\}\cdot \{cos\alpha ,cos\beta \}=|gradf|\cdot |\mathbf{l}|\cdot \cos \theta\;\;\;\;\;\;\;\;(2.13)

$\bf{l}$ $\alpha,\beta$ $\bf{l}$ $x$ $y$ $\theta$ $\bf{l}$ 的夹角。

$\theta$ $0$ 。由此可知，只有当某点处方向导数的方向与梯度的方向一致时，方向导数在该点才会取得最大的变化率。

$z=x^2+y^2+5$ $A$ $(-3,3,23)$ $\partial z/\partial x=2x,\partial z/\partial y=2y$ $A$ $(-6,6)$ $A$ $(\sqrt{2}/2,-\sqrt{2}/2)$ $\Delta z$ 。

如图2-10所示，要想每次都以最快的速度下降，则每次都必须向梯度的反方向向前跳。

图2-10 负梯度方向图

2.5.2梯度下降算法

介绍这么多总算是把梯度这回事儿给说清楚了，那么如何用具体的数学表达式进行描述呢？总不能一个劲儿的喊它“跳”对吧。为了方便后面的表述以及将读者带入到一个真实求解的过程中，这里先将图2-9中的字母替换成模型中的参数表述。

$J(w_1,w_2)=w_1^2+w_2^2+2w_2+5$ $b$ $w_1,w_2$ $A$ 为初始权重值。下面就一步步的通过梯度下降法来进行求解。

$A=(w_1,w_2)=(-2,3)$ $J(-2,3)=24$ $A$ $-grad\ J=-(2{{w}_{1}},2{{w}_{2}}+2)=(1,-2)$ 。

图2-11 梯度下降图

$\bf{OQ}$ $\bf{AP}$ $\alpha$ $A$ $A$ $\bf{AP}$ $A$ $P$ $\bf{AP}$ $A$ $P$ 点的位置。

图2-12 梯度计算图

$\bf{AP},\bf{OA},\bf{OP}$ 三者的关系为

\mathbf{OP}=\mathbf{OA}-\mathbf{PA} \;\;\;\;\;\;\; \;\;\;\;(2.14)

$(2.14)$ 改写成

\mathbf{OP}=\mathbf{OA}-\alpha \cdot grad\ J\;\;\;\;\;\;\;(2.15)

$\bf{OP},\bf{OA}$ $w_1,w_2$ 更新后与跟新前的值，所以便可以得出梯度下降的更新公式

W=W-\alpha \cdot \frac{\partial J}{\partial W}\;\;\;\;\;\;\;\;(2.16)

$W=(w_1,w_2)$ $\partial J/\partial W$ $\alpha$ $A$ 在第一次跳跃后的着落点为

\begin{aligned} & {{w}_{1}}={{w}_{1}}-0.1\times 2\times {{w}_{1}}=-2-0.1\times 2\times (-2)=-1.6 \\ & {{w}_{2}}={{w}_{2}}-0.1\times (2\times {{w}_{2}}+2)=3-0.1\times (2\times 3+2)=2.2 \\ \end{aligned} \;\;\;\;\;\;\;\;(2.17)

$(-2,3)$ $(-1.6,2.2)$ $J(w_1,w_2)$ $A$ $P$ 之后，又可以再次利用梯度下降算法进行跳跃，直到跳到谷底（或附近），如图2-13所示。

图2-13 梯度下降示意图

最后，根据上述原理，还可以通过实际的代码将整个过程展示出来，完整代码见Chapter02/05_gradient_descent_visualization.py文件。


x
1
def cost_function(w1, w2):
2
    J = w1 ** 2 + w2 ** 2 + 2 * w2 + 5
3
    return J
4
def compute_gradient(w1, w2):
5
    return [2 * w1, 2 * w2 + 2]
6
7
def gradient_descent():
8
    w1, w2 = -2, 3
9
    jump_points = [[w1, w2]]
10
    costs = [cost_function(w1, w2)]
11
    step = 0.1
12
    print("P:({},{})".format(w1, w2), end=' ')
13
    for i in range(20):
14
        gradients = compute_gradient(w1, w2)
15
        w1 = w1 - step * gradients[0]
16
        w2 = w2 - step * gradients[1]
17
        jump_points.append([w1, w2])
18
        costs.append(cost_function(w1, w2))
19
        print("P{}:({},{})".format(i + 1, round(w1, 3), round(w2, 3)), end=' ')
20
    return jump_points, costs
21
22
if __name__ == '__main__':
23
    jump_points, costs = gradient_descent()
24
    plot_surface_and_jump_points(jump_points, costs)
25
    
26
# 部分结果：P:(-2,3) P1:(-1.6,2.2) P2:(-1.28,1.56).....P20:(-0.023,-0.954)

$(-0.023,-0.954)$ 附近，如图2-14所示。此致，就介绍完了如何通过编码来实现梯度下降算法的求解过程，等后续完成线性回归模型的推导后，再来自己编码完成线性回归模型的参数求解过程。

图2-14 梯度下降可视化图

2.5.3 小结

在本节中，笔者通过一个跳跃的例子详细给大家介绍了什么是梯度，以及为什么要沿着梯度的反方向进行跳跃；然后通过图示导出了梯度下降的更新公式

W=W-\alpha \cdot \frac{\partial J}{\partial W}\;\;\;\;\;\;\;\;(2.18)

$\alpha$ 也是一个十分重要的参数，这将在第3章中进行介绍。

2.6正态分布

2.6.1一个问题的出现

$\sum{{{(\hat{y}-y)}^{2}}}$ $\hat{y}$ $y$ $\theta$ ${{x}_{1}},{{x}_{2}},...,{{x}_{n}}$ $n$ $e_{i}={x}_{i}-\theta$ ，按照勒让德提出的方法，累计误差为

E(\theta )=\sum\limits_{i=1}^{n}{e_{i}^{2}}=\sum\limits_{i=1}^{n}{{{({{x}_{i}}-\theta )}^{2}}}\;\;\;\;\;\;\;\;(2.19)

$E(\theta)$ $\bar{x}=1/n\sum x_{i}$ 。由于算术平均是一个历经考验的方法，而以上的推理说明从另一个角度也说明了最小二乘法的优良性。这使得当时的人们对于最小二乘法有了更强的信心。从这里可以看出，这种做法的逻辑是，首先认为算术平均这种做法好但不知道为什么，然后有人提出了一种衡量误差的方法最小二乘，接着对误差最小化求解后发现其解正是算术平均，所以肯定了最小二乘的有用性。但事实上就是既没有说清楚算术平均为什么好，反而用算术平均的结果来肯定了最小二乘的作用。

$f(x)$ $x=0$ $f(x)$ $|x|$ 增大而减小，如图2-15所示。于是许多天文学家和数学家开始了寻找误差分布曲线的尝试，但最终都没能给出有用的结果。

图2-15 理想状态下误差分布图

2.6.2 正态分布

1801年1月，天文学家朱塞普·皮亚齐发现了一颗从未见过的光度为8等的星在移动，这颗现在被称作谷神星（Ceres）的小行星在夜空中出现 6 个星期，扫过八度角后就在太阳的光芒下没了踪影无法观测。由于留下的观测数据有限难以计算出它的轨迹，所以天文学家们也因此无法确定这颗新星是彗星还是行星。不过这个问题很快成了学术界关注的焦点。高斯当时已经是很有名望的年轻数学家，这个问题引起了他的兴趣。高斯以其卓越的数学才能创立了一种崭新的行星轨道的计算方法，一个小时之内就计算出了谷神星的轨道，并预言了他在夜空中出现的时间和位置。1801年12月31日夜，德国天文爱好者奥伯斯，在高斯预言的时间里，用望远镜对准了这片天空，果然不出所料，谷神星出现了。

高斯为此名声大震，但是高斯当时拒绝透露计算轨道的方法，原因可能是高斯认为自己的方法的理论基础还不够成熟。直到1809年高斯系统地完善了相关的数学理论后，才将他的方法公布于众，而其中使用的数据分析方法，就是以正态误差分布为基础的最小二乘法。那高斯是如何推导出误差分布为正态分布的呢？

$\theta$ ${{x}_{1}},{{x}_{2}},...,{{x}_{n}}$ $n$ $e_i=x_i - \theta$ $e_i$ $f(x)$ $n$ 个误差同时出现的概率记为

L(\theta )=L(\theta ;{{x}_{1}},{{x}_{2}},...,{{x}_{n}})=f({{e}_{1}})f({{e}_{2}})\cdots f({{e}_{n}})\;\;\;\;\;\;(2.20)

$L(\theta)$ $\hat{\theta}$ $\theta$ $(2.21)$ $\hat{\theta}$ 值。

L(\hat{\theta })=\arg \underset{\theta }{\mathop{\max }}\,L(\theta )\;\;\;\;\;\;\;\;(2.21)

$L(\theta)$ $\hat{\theta}$ $\theta$ 的极大似然估计。在这里高斯首次给出了极大似然的思想，这个思想后来被统计学家费希尔系统的发展成为参数估计中的极大似然估计理论。同时，所谓最大似然估计是指在已知样本结果的情况下，推断出最有可能使得该结果出现的参数的过程。也就是说最大似然估计一个过程，它用来估计出某个模型的参数，而这些参数能使得已知样本的结果最可能发生。

$f(x)$ $f(x)$ ，并证明所有的概率密度函数中，唯一满足这个性质的就是

f(x)=\frac{1}{\sqrt{2\pi }\sigma }{{e}^{-\frac{{{x}^{2}}}{2{{\sigma }^{2}}}}}\;\;\;\;\;\;\;\;(2.22)

$\sigma>0$ 为常数，而这也就是正态分布。

$e_i$ ，由式可知其对应是似然估计为

L(\theta )=\prod\limits_{i=1}^{n}{f}({{e}_{i}})=\frac{1}{{{(\sqrt{2\pi }\sigma )}^{n}}}\exp \left\{ -\frac{1}{2{{\sigma }^{2}}}\sum\limits_{i=1}^{n}{e_{i}^{2}} \right\}\;\;\;\;\;\;\;\;(2.23)

$L(\theta)$ $\sum\nolimits_{i=1}^{n}{e_{i}^{2}}$ 取值最小，显然这正好就是最小二乘法的要求。可以看出，高斯这种做法的初始动机仍旧是以算术平均作为一种“公理”，然后以此为基础做出假设找到一种符合人们常识的误差密度函数，即正态分布。最后再通过最大似然估计来印证了最小二乘法。

2.7 目标函数推导

经过前面几节内容的介绍，我们知道了什么是线性回归、怎么转换求解问题、如何通过sklearn进行建模求解以及梯度下降法的原理与推导。同时，在2.6节中还通过一个故事来交代了最小二乘法的来历，以及误差服从高斯分布的事实。下面，我们就来完成线性回归中的最后两个任务，线性回归的推导以及Python代码的实现。

2.7.1 目标函数

$(x^{(i)},y^{(i)})$ ${{\hat{y}}^{(i)}}={{W}^{T}}{{x}^{(i)}}+b$ ，则有

{{y}^{(i)}}={{\hat{y}}^{(i)}}+{{e}^{(i)}}\;\;\;\;\;\;\;\;(2.24)

$e^{(i)}$ $i$ $W$ $x^{(i)}$ $b$ $e^{(i)}$ $0$ 的高斯分布[3]，于是有

f({{e}^{(i)}})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left( -\frac{{{({{e}^{(i)}}-0)}^{2}}}{2{{\sigma }^{2}}} \right)\;\;\;\;\;\;\;\;(2.25)

$(2.24)$ $(2.25)$ 有

f({{e}^{(i)}})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left( -\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right)\;\;\;\;\;\;\;(2.26)

$(2.26)$ $y^{(i)}$ $y^{(i)}$ $\hat{y}^{(i)}$ $W,b,{{x}^{(i)}}$ ${{y}^{(i)}}$ $x^{(i)},W,b$ 下的条件分布），于是有

p({{y}^{(i)}}|{{x}^{(i)}};W,b)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left( -\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right)\;\;\;\;\;\;\;\;(2.27)

$y^{(i)}$ $\hat{y}^{(i)}$ $\sigma^2$ $\hat{y}^{(i)}$ $W,b$ $W,b$ 能够使得已知的真实值最容易发生呢？此时就要用到极大似然估计来进行参数估计

L(\theta )=\prod\limits_{i=1}^{m}{p}({{y}^{(i)}}|{{x}^{(i)}};W,b)=\prod\limits_{i=1}^{m}{\frac{1}{\sqrt{2\pi }\sigma }}\exp \left( -\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right)\;\;\;\;\;\;\;\;(2.28)

$(2.28)$ 的两边同时取自然对数

\begin{aligned} & \log L(W,b)=\log \left\{ \prod\limits_{i=1}^{m}{\frac{1}{\sqrt{2\pi }\sigma }}\exp \left( -\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right) \right\} \\ & =\sum\limits_{i=1}^{m}{\log }\left\{ \frac{1}{\sqrt{2\pi }\sigma }\exp \left( -\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right) \right\} \\ & =\sum\limits_{i=1}^{m}{\left\{ \log \frac{1}{\sqrt{2\pi }\sigma }-\frac{{{({{y}^{(i)}}-{{{\hat{y}}}^{(i)}})}^{2}}}{2{{\sigma }^{2}}} \right\}} \\ & =m\cdot \log \frac{1}{\sqrt{2\pi }\sigma }-\frac{1}{{{\sigma }^{2}}}\frac{1}{2}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-{{{\hat{y}}}^{(i)}} \right)}^{2}}} \end{aligned}\;\;\;\;\;\;\;\;(2.29)

$\max L(W,b)$ $\max\log L(W,b)$ ，所以

\begin{aligned} & \max \log L(W,b)\Leftrightarrow \min \frac{1}{{{\sigma }^{2}}}\frac{1}{2}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-{{{\hat{y}}}^{(i)}} \right)}^{2}}} \\ & \Leftrightarrow \min \frac{1}{2}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-{{{\hat{y}}}^{(i)}} \right)}^{2}}} \end{aligned}\;\;\;\;\;\;\;\;(2.30)

于是得目标函数

J(W,b)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-{{{\hat{y}}}^{(i)}} \right)}^{2}}}\;\;\;\;\;\;\;(2.31)

2.7.2 求解梯度

$y^{(i)}$ $i$ $\hat{y}^{(i)}$ $i$ $W$ $W_j$ $X$ $m\times n$ $m$ $n$ $x^{(i)}$ $i$ $x^{(i)}_j$ $j$ 维特征。

J(W,b)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-{{{\hat{y}}}^{(i)}} \right)}^{2}}}=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)}^{2}}}\;\;\;\;\;\;\;\;(2.32)

${{W}_{j}}$ 的梯度求解过程为

\begin{aligned} & \frac{\partial J}{\partial {{W}_{j}}}=\frac{\partial }{\partial {{W}_{j}}}\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-({{W}_{1}}x_{1}^{(i)}+{{W}_{2}}x_{2}^{(i)}\cdots {{W}_{n}}x_{n}^{(i)}+b) \right)}^{2}}} \\ & =\frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}-({{W}_{1}}x_{1}^{(i)}+{{W}_{2}}x_{2}^{(i)}\cdots {{W}_{n}}x_{n}^{(i)}+b) \right)}\cdot (-x_{j}^{(i)}) \\ & =\frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)}\cdot (-x_{j}^{(i)}) \end{aligned}\;\;\;\;\;\;\;\;(2.33)

$b$ 的梯度求解过程为

\begin{aligned} & \frac{\partial J}{\partial b}=\frac{\partial }{\partial b}\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)}^{2}}} \\ & =-\frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)} \end{aligned} \;\;\;\;\;\;\;\;(2.34)

此时便得到了目标函数关于参数的梯度计算公式

\begin{aligned} & J(W,b)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)}^{2}}} \\ & \frac{\partial J}{\partial {{W}_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)}\cdot (-x_{j}^{(i)}) \\ & \frac{\partial J}{\partial b}=-\frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}-({{W}^{T}}{{x}^{(i)}}+b) \right)} \end{aligned} \;\;\;\;\;\;\;\;(2.35)

在推导得到每个参数的梯度更新公式后，就能够根据梯度下降算法来迭代的对参数值进行更新，直到目标函数收敛。到此，对于整个线性回归部分阶段二的内容就介绍完了，下面继续来看最后一个阶段的内容。

2.7.3 矢量化计算

$(2.35)$ 进行矢量化，代码如下：


xxxxxxxxxx
4
1
y_hat=np.matmul(X, W) + b
2
J(W,b)= 0.5 * (1 / m) * np.sum((y – y_hat) ** 2)
3
grad_w=-1/m*np.matmul(X.T,(y-y_hat))
4
grad_b=-1/m*np.sum(y-y_hat)

$(2.35)$ ${{W}_{j}}$ $W$ $X$ $W$ $X$ 的转置放在前面。

2.7.4 从零实现线性回归

下面依旧以波士顿房价预测模型为例，然后通过手动编写代码进行模型的建模与求解，完整代码见Chapter02/06_boston_price_train.py文件。

1) 定义预测函数

$y={{w}_{1}}{{x}_{1}}+{{w}_{2}}{{x}_{2}}\cdots {{w}_{n}}{{x}_{n}}+b$ 。为了同时对输入的所有样本进行计算，一般以两个矩阵相乘的方式进行，代码如下：


xxxxxxxxxx
2
1
def prediction(X, W, bias):#预测
2
    return np.matmul(X, W) + bias  # [m,n] @ [n,1] = [m,1]

2) 定义目标函数

为了方便在训练结束后（或者训练时）观察目标函数的收敛情况，所以通常会计算每一次参数更新后的损失值，代码如下：


xxxxxxxxxx
4
1
def cost_function(X, y, W, bias):  # 代价函数
2
    m, n = X.shape
3
    y_hat = prediction(X, W, bias)
4
    return 0.5 * (1 / m) * np.sum((y - y_hat) ** 2)

3) 定义梯度下降

在这里，需要完成模型训练时的核心部分，也就是整个梯度下降的过程，代码如下：


xxxxxxxxxx
8
1
def gradient_descent(X, y, W, bias, alpha):
2
    m, n = X.shape
3
    y_hat = prediction(X, W, bias)
4
    grad_w = -(1 / m) * np.matmul(X.T, (y - y_hat))
5
    grad_b = -(1 / m) * np.sum(y - y_hat) 
6
    W = W - alpha * grad_w  
7
    bias = bias - alpha * grad_b
8
    return W, bias

可以看到，对于梯度的求解就是式矢量化后的形式，然后运用梯度下降算法对参数更新即可。

4) 训练模型

在定义完训练过程中需要的相关函数后，接着就可以初始化相关权重和变量通过梯度下降算法来进行参数的更新，代码如下：


xxxxxxxxxx
8
1
def train(X, y, ite=200):
2
    m, n = X.shape  # 506,13
3
    W, b, alpha, costs = np.random.randn(n, 1),0.1,0.2
4
    for i in range(ite):
5
        costs.append(cost_function(X, y, W, b))
6
        W, b = gradient_descent(X, y, W, b, alpha)
7
    y_pre = prediction(X, W, b)
8
    return costs

最后，可以通过如下所示的过程来完成整个模型的训练，代码如下：


xxxxxxxxxx
7
1
if __name__ == '__main__':
2
    x, y = load_data()
3
    train_by_sklearn(x, y)
4
    costs = train(x, y)
5
    #输出结果：
6
    #MSE:  21.894831181729206
7
    #MSE:  21.901070160392404

同时，在这里也加入了通过sklearn训练后的模型的MSE评估值。可以看出，我们自己实现的模型与sklearn中的线性回归模型在MSE上几乎没有任何差别。大约在25次迭代后目标函数就开始进入收敛状态，如图2-16所示。

图2-16 目标函数收敛图

公众号后台回复“跟我一起机器学习”即可获得高清PDF与教学PPT！

2.7.5 小结

通过本节内容的学习，对于线性回归的主要内容也算到此结束了。对于其它细枝末节的地方（例如学习率的选择、特征标准化等），在后续模型改善部分再进行介绍。

总结一下，如图2-17所示，在本章中笔者首先介绍了线性回归模型第一阶段的主要内容，至此你可以大致知道拿到一个问题如何通过开源的sklearn库进行线性回归建模。虽然阶段一的内容并没有从数学的角度来完成对于线性回归的论证，但是其包含了学习一个算法并快速入门的路线。然后笔者接着介绍了梯度下降算法的介绍、误差的正态分布特性和目标函数的推导，完成了线性回归第二阶段的学习。这个阶段一般都是算法从数学层面上的理论依据，尤其是统计机器学习更是依赖于这个过程，但是这部分通常来说都有一定的难度，例如SVM的求解过程。最后，笔者通过本节的内容完成了线性回归模型源码的实现。

图2-17 学习层次图

本次内容就到此结束，感谢您的阅读！如果你觉得上述内容对你有所帮助，欢迎分享至一位你的朋友！若有任何疑问与建议，请添加笔者微信'nulls8'或加群进行交流。青山不改，绿水长流，我们月来客栈见！

引用

[1] Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.

[2] 靳志辉; 数学文化, 4 (2013), pp. 36-47.

[3] Andrew Ng, Machine Learning, Stanford University, CS229, Spring 2019.

[4] Machine Learning – MT 2016 3. Maximum Likelihood

[5] 示例代码 https://github.com/moon-hotel/MachineLearningWithMe/

第2章 从零认识线性回归（附高清PDF和PPT）

于2021年10月9日2021年10月9日由 空字符 发布