跟我一起学scikit-learn16：线性回归算法

线性回归算法是使用线性方程对数据集拟合得算法，是一个非常常见的回归算法。本章首先从最简单的单变量线性回归算法开始介绍，然后介绍了多变量线性回归算法，其中成本函数以及梯度下降算法的推导过程会用到部分线性代数和偏导数；接着重点介绍了梯度下降算法的求解步骤以及性能优化方面的内容；最后通过一个房价预测模型，介绍了线性回归算法性能优化的一些常用步骤和方法。

1.单变量线性回归算法

我们先考虑最简单的单变量线性回归算法，即只有一个输入特征。

1.预测函数

针对数据集x和y，预测函数会根据输入特征x来计算输出值h(x)。其输入和输出的函数关系如下：
$$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$$
这个方程表达的是一条直线。我们的任务是构造一个$h_{\theta }$函数，来映射数据集中的输入特征x和输出值y，使得预测函数$h_{\theta }$计算出来的值与真实值y的整体误差最小。构造$h_{\theta }$函数的关键就是找到合适的${\theta }_0$和${\theta }_1$的值，${\theta }_0$和${\theta }_1$成为模型参数。

假设我们有如下的数据集：

假设模型参数${\theta }_0=1$，${\theta }_1=3$，则预测函数为$h_{\theta }(x)=1+3x$。针对数据集中的第一个样本，输入为1，根据模型函数预测出来的值是4，与输出值y是吻合的。针对第二个样本，输入为2，根据模型函数预测出来的值是7，与实际输出值y相差1。模型的求解过程就是找出一组最合适的模型参数${\theta }_0$和${\theta }_1$，以便能最好地拟合数据集。

怎样来判断最好地拟合了数据集呢？没错，就是使用成本函数（也叫损失函数）。当拟合成本最小时，即找到了最好的拟合参数。

2.成本函数

单变量线性回归算法的成本函数是：
$$J(\theta )=J({\theta }_0,{\theta }_1)=\frac{1}{2m}\sum _{i=1}^m(h(x^{(i)})-y^{(i)}{)}^2$$
其中，$h(x^{(i)})-y^{(i)}$是预测值和真实值之间的误差，故成本就是预测值和真实值之间误差平方的平均值，之所以乘以1/2是为了方便计算。这个函数也称为均方差公式。有了成本函数，就可以精确地测量模型对训练样本拟合的好坏程度。

3.梯度下降算法

有了预测函数，也可以精确地测量预测函数对训练样本的拟合情况。但怎么求解模型参数${\theta }_0$，${\theta }_1$的值呢？这时梯度下降算法就排上了用场。

我们的任务是找到合适的${\theta }_0$，${\theta }_1$，使得成本函数$J({\theta }_0,{\theta }_1)$最小。为了便于理解，我们切换到三维空间来描述这个任务。在一个三维空间里，以${\theta }_0$作为x轴，以${\theta }_1$作为y轴，以成本函数$J({\theta }_0,{\theta }_1)$作为z轴，那么我们的任务就是要找出当z轴上的值最小的时候所对应的x轴上的值和y轴上的值。

梯度下降算法的原理是：先随机选择一组${\theta }_0$，${\theta }_1$，同时选择一个参数$\alpha$作为移动的步长。然后，让x轴上的${\theta }_0$和y轴上的${\theta }_1$分别向特定的方向移动一小步，这个步长的大小就由参数$\alpha$决定。经过多次迭代之后，x轴和y轴上的值决定的点就慢慢靠近z轴上的最小值处，如图所示。

这是个等高线图，就是说在我们描述的三维空间里，你的视角在正上方，看到一圈一圈z轴值相同的点构成的线。在上图中，随机选择的点在$X_0$处，经过多次迭代后，慢慢地靠近圆心处，即z轴上最小值附近。

问题来了，$X_0$（由$[{\theta }_0,{\theta }_1]$描述）怎么知道往哪个方向移动，才能靠近z轴上最小值附近？答案是往成本函数逐渐变小的方向移动。怎么表达成本函数逐渐变小的方向呢？答案是偏导数。

可以简单地把偏导数理解为斜率。我们要让${\theta }_j$不停地迭代，由当前${\theta }_j$的值，根据$J(\theta )$的偏导数函数，算出$J(\theta )$在${\theta }_j$上的斜率，然后在乘以学习率$\alpha$，就可以让${\theta }_j$往前$J(\theta )$变小的方向迈一小步。

用数学来描述上述过程，梯度下降的公式为：
$${\theta }_j={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )$$
公式中，下标j就是参数的序号，针对单变量线性回归，即0和1。$\alpha$称为学习率，它决定每次要移动的幅度大小，它会乘以成本函数对参数${\theta }_j$的偏导数，以这个结果作为参数移动的幅度。如果幅度太小，就意味着要计算很多次才能到达目的地；如果幅度太大，可能会直接跨过目的地，从而无法收敛。

把成本函数$J(\theta )$的定义代入上面的公式中，不难推导出梯度下降算法公式：
$${\theta }_0={\theta }_0-\frac{\alpha }{m}\sum _{i=1}^m(h(x^{(i)})-y^{(i)})$$
$${\theta }_1={\theta }_1-\frac{\alpha }{m}\sum _{i=1}^m((h(x^{(i)})-y^{(i)})x_i)$$
公式中，$\alpha$是学习率；m是训练样本的个数；$h(x^{(i)})-y^{(i)}$是模型预测值和真实值的误差。需要注意的是，针对${\theta }_0$和${\theta }_1$分别求出了其迭代公式，在${\theta }_1$的迭代公式里，累加器中还需要乘以$x_i$。对公式推导感兴趣的读者，可以参考本章扩展部分的内容。

2.多变量线性回归算法

工程应用中往往不止一个输入特征。熟悉了单变量线性回归算法后，我们来探讨一下多变量线性回归算法。

1.预测函数

上文介绍的线性回归模型里只有一个输入特征，我们推广到更一般的情况，即多个输入特征。此时输出y的值由n个输入特征x1,x2,x3,…,xn决定。那么预测函数模型可以改写如下：
$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n$$
假设x0=1，那么上面的公式可以重写为：
$$h_{\theta }(x)=\sum _{j=0}^n{\theta }_j{x}_j$$
其中，${\theta }_0$，${\theta }_1$，…，${\theta }_n$统称为$\theta$，是预测函数的参数。即一组$\theta$值就决定了一个预测函数，记为$h_{\theta }(x)$，为了简便起见，在不引起误解的情况下可以简写为$h(x)$。理论上，预测函数有无穷多个，我们求解的目标就是找出一个最优的$\theta$值。

思考：当有n个变量$x_1$，$x_2$，…，$x_n$决定y值的时候，训练数据集应该长什么样呢？

（1）向量形式的预测函数
根据向量乘法运算法则，成本函数可重写为：
$$h_{\theta }(x)=\left[\begin{array}{c}{\theta }_0,{\theta }_1,\cdots ,{\theta }_n\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_n\end{array}\right]={\theta }^{T}x$$
此处，依然假设$x_0=1$，$x_0$称为模型偏置（bias）。

为什么要写成向量形式的预测函数呢？一是因为简洁，而是因为在实现算法时，要用到数值计算里的矩阵运算来提高效率，比如Numpy库里的矩阵运算。

（2）向量形式的训练样本
假设，输入特征的个数是n，即$x_1$，$x_2$，…，$x_n$，我们总共有m个训练样本，为了书写方便，假设$x_0=1$。这样训练样本可以写成矩阵的形式，即矩阵里每一行都是一个训练样本，总共有m行，每行有n+1列。

思考：为什么不是n列而是n+1列？答案是：把模型偏置$x_0$也加入了训练样本里。最后把训练样本写成一个矩阵，如下：

$$X=\left[\begin{array}{ccccc}{x}_0^{(1)}& x_1^{(1)}& x_2^{(1)}& \cdots & x_n^{(1)}\\ x_0^{(2)}& x_1^{(2)}& x_2^{(2)}& \cdots & x_n^{(2)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_0^{(m)}& x_1^{(m)}& x_2^{(m)}& \cdots & x_n^{(m)}\end{array}\right],\theta =\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ ⋮\\ {\theta }_n\end{array}\right]$$

理解训练样本矩阵的关键在于理解这些上标和下标的含义。其中，带括号的上标表示样本序号，从1到m；下标表示特征序号，从0到n，其中$x_0$为常数1。比如，$x_j^{(i)}$表示第i个训练样本的第j个特征的值。而$x^{(i)}$只有上标，则表示第i个训练样本所构成的列向量。

熟悉矩阵乘法的话不难得出结论，如果要一次性计算出所有训练样本的预测值$h_{\theta }(X)$，可以使用下面的矩阵运算公式：
$$h_{\theta }(X)=X\theta$$
从这个公式也可以看到矩阵形式表达的优势。实际上，在scikit-learn里，训练样本就是用这种方式表达的，即使用$m\times n$维的矩阵来表达训练样本，可以回顾一下scikit-learn里模型的fit()函数的参数。

2.成本函数

多变量线性回归算法的成本函数：
$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h(x^{(i)})-y^{(i)}{)}^2$$
其中，模型参数$\theta$为n+1维的向量，$h(x^{(i)})-y^{(i)}$是预测值和实际值的差。这个形式和单变量线性回归算法的类似。

成本函数有其对应的矩阵形式：
$$J(\theta )=\frac{1}{2m}(X\theta -\vec{y}{)}^T(X\theta -\vec{y})$$
其中，X为$m\times (n+1)$维的训练样本矩阵；上标T表示转置矩阵；$\vec{y}$表示由所有的训练样本的输出$y^{(i)}$构成的向量。这个公式的优势是：没有累加器，不需要循环，直接使用矩阵运算，就可以一次性计算出对特定的参数$\theta$下模型的拟合成本。

思考：矩阵运算真的不需要循环吗？
这里所说的不需要循环，是指不需要在算法实现层使用循环，但在数值运算库如Numpy里，实现的矩阵运算还是要用到循环。虽然都是循环，但是有差别，一是在数值运算库里实现的循环效率更高，而是矩阵运算的循环可以使用分布式来实现。一个大矩阵运算可以拆成多个子矩阵运算，然后在不同的计算机上执行运算，最终再把运算结果汇合起来。这种分布式计算对大型矩阵运算来说是一种必要的手段。

3.梯度下降算法

根据单变量线性回归算法的介绍，梯度下降的公式为：
$${\theta }_j={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )$$
公式中，下标j是参数的序号，其值从0到n；$\alpha$为学习率。把成本函数代入上式，利用偏导数计算法则，不难推导出梯度下降算法的参数迭代公式：
$${\theta }_j={\theta }_j-\frac{\alpha }{m}\sum _{i=1}^m((h(x^{(i)})-y^{(i)})x_j^{(i)})$$

我们可以对比一下单变量线性回归函数的参数迭代公式。实际上和多变量线性回归函数的参数迭代公式是一模一样的。惟一的区别就是因为$x_0$为常数1，在单变量线性回归算法的参数迭代公式中省去了。

这个公式怎么样用编程语言来实现呢？在编写机器学习算法的时候，一般步骤如下。

（1）确定学习率：$\alpha$太大可能会使成本函数无法收敛，太小则计算太多，机器学习算法效率就比较低。

（2）确定参数起始点：比如让所有的参数都以1作为起始点，即${\theta }_0=1$，${\theta }_1=1$，…，${\theta }_n=1$。这样就得到了我们的预测函数$h_{\theta }(x)={\sum }_{i=1}^m{x}^{(i)}$。根据预测值和成本函数，就可以算出在参数起始位置的成本。需要注意的是，参数起始点可以根据实际情况灵活选择，以便让机器学习算法的性能更高，比如选择比较靠近极点的位置。

（3）计算参数的下一组值：根据梯度下降参数迭代公式，分别同时计算出新的${\theta }_j$的值。然后用新的$\theta$值得到新的预测函数$h_{\theta }(x)$。再根据新的预测函数，代入成本函数就可以算出新的成本。

（4）确定成本函数是否收敛：拿新的成本和旧的成本进行比较，看成本是不是变得越来越小。如果两次成本之间的差异小于误差范围，即说明已经非常靠近最小成本了，就可以近似地认为我们找到了最小成本。如果两次成本之间的差异在误差范围之外，重复步骤（3）继续计算下一组参数$\theta$，直到找到最优解。

3.模型优化

线性回归模型常用的优化方法，包括增加多项式特征以及数据归一化处理等。

1.多项式与线性回归

当线性回归模型太简单导致欠拟合时，我们可以增加特征多项式来让线性回归模型更好地拟合数据。比如有两个特征$x_1$和$x_2$，可以增加两个特征的乘积$x_1\times x_2$作为新特征$x_3$。同理，我们也可以增加$x_1^2$和$x_2^2$分别作为新特征$x_4$和$x_5$。

在scikit-learn里，线性回归是由类sklearn.learn_model.LinearRegression实现的，多项式由类sklearn.preprocessing.PolynomialFeatures实现。那么要怎样添加多项式特征呢？我们需要用一个管道把两个类串起来，即用sklearn.pipeline.Pipeline把这两个模型串起来。

比如下面的函数就可以创建一个多项式拟合：

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
def polynomial_model(degree=1):
    polynomial_features = PolynomialFeatures(degree=degree,include_bias=False)
    linear_regression = LinearRegression(normalize=True)
    # 这是一个流水线，先增加多项式阶数，然后再用线性回归算法来拟合数据
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("linear_regression", linear_regression)])
    return pipeline

一个Pipeline可以包含多个处理节点，在scikit-learn里，除了最后一个节点外，其他的节点都必须实现fit()方法和transform()方法，最后一个节点只需要实现fit()方法即可。当训练样本数据送进Pipeline里进行处理时，它会逐个调用节点的fit()方法和transform()方法，最后调用最后一个节点的fit()方法来拟合数据。管道的示意图如下所示：

2.数据归一化

当线性回归模型有多个输入特征时，特别是使用多项式添加特征时，需要对数据进行归一化处理。比如，特征$x_1$的范围在[1,4]之间，特征$x_2$的范文在[1,2000]之间，这种情况下，可以让$x_1$除以4来作为新特征$x_1$，同时让$x_2$除以2000来作为新特征$x_2$，该过程称为特征缩放（feature scaling）。可以使用特征缩放来对训练样本进行归一化处理，处理后的特征范围在[0,1]之间。

为什么要进行数据归一化处理？以及归一化处理有哪些注意事项？
归一化处理的目的是让算法收敛更快，提升模型拟合过程中的计算效率。进行归一化处理后，当有个新的样本需要计算预测值时，也需要先进行归一化处理，再通过模型来计算预测值，计算出来的预测值要再乘以归一化处理的系数，这样得到的数据才是真正的预测数据。

在scikit-learn里，使用LinearRegression进行线性回归时，可以指定normalize=True来对数据进行归一化处理。具体可以查阅scikit-learn文档。

4.示例：使用线性回归算法拟合正弦函数

本节用线性回归算法来模拟正弦函数。
首先生成200个在区间$[2\pi ,2\pi ]$内的正弦函数上的点，并给这些点加上一些随机的噪声。

import numpy as np
n_dots = 200
X = np.linspace(-2 * np.pi, 2 * np.pi, n_dots)
Y = np.sin(X) + 0.2 * np.random.rand(n_dots) - 0.1
X = X.reshape(-1, 1)
Y = Y.reshape(-1, 1);

其中，reshape()函数的作用是把Numpy的数组转换成符合scikit-learn输入格式的数组（这里是把一个n维向量转换成一个n*1维的矩阵），否则scikit-learn会报错。

接着我们使用PolynomialFeatures和Pipeline创建一个多项式拟合模型：

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
def polynomial_model(degree=1):
    polynomial_features = PolynomialFeatures(degree=degree,include_bias=False)
    linear_regression = LinearRegression(normalize=True)
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("linear_regression", linear_regression)])
    return pipeline

分别用2/3/5/10阶多项式来拟合数据集：

from sklearn.metrics import mean_squared_error
degrees = [2, 3, 5, 10]
results = []
for d in degrees:
    model = polynomial_model(degree=d)
    model.fit(X, Y)
    train_score = model.score(X, Y)
    mse = mean_squared_error(Y, model.predict(X))
    results.append({"model": model, "degree": d, "score": train_score, "mse": mse})
for r in results:
    print("degree: {}; train score: {}; mean squared error: {}"
          .format(r["degree"], r["score"], r["mse"]))

算出每个模型拟合的评分，此外，使用mean_squared_error算出均方根误差，即实际的点和模型预点之间的距离，均方根误差越小说明模型拟合效果越好——上述代码的输出结果为：

degree: 2; train score: 0.14988484858514306; mean squared error: 0.4346164668605285
degree: 3; train score: 0.27617240286045885; mean squared error: 0.37005268328809543
degree: 5; train score: 0.89946485219823; mean squared error: 0.05139801432804186
degree: 10; train score: 0.9941202722467253; mean squared error: 0.0030059768938090893

从输出结果可以看出，多项式阶数越高，拟合评分越高，均方根误差越小，拟合效果越好。

最后我们把不同模型的拟合效果在二维坐标上画出来，可以清楚地看到不同阶数的多项式的拟合效果：

import matplotlib.pyplot as plt
from matplotlib.figure import SubplotParams
plt.figure(figsize=(12,6),dpi=200,subplotpars=SubplotParams(hspace=0.3))
for i,r in enumerate(results):
    fig = plt.subplot(2,2,i+1)
    plt.xlim(-8,8)
    plt.title("LinearRegression degree={}".format(r["degree"]))
    plt.scatter(X,Y,s=5,c='b',alpha=0.5)
    plt.plot(X,r["model"].predict(X),'r-')
plt.show()

我们使用SubplotParams调整了子图的竖直间距，并且使用subplot()函数把4个模型的拟合情况都画在同一个图形上。上述代码的输出结果如下图所示：

思考：在[-2π，2π]区间内，10阶多项式对数据拟合得非常好，我们可以试着画出这10阶模型在[-20,20]的区域内的曲线，观察一下该模型的曲线和正弦函数的差异。代码如下：

plt.figure(figsize=(12,6),dpi=200)
X = np.linspace(-20,20,2000).reshape(-1, 1)
Y = np.sin(X).reshape(-1, 1)
model_10 = results[3]["model"]
plt.xlim(-20,20)
plt.ylim(-2,2)
plt.plot(X,Y,'b-')
plt.plot(X,model_10.predict(X),'r-')
dot1 = [-2*np.pi,0]
dot2 = [2*np.pi,0]
plt.scatter(dot1[0],dot1[1],s=50,c='r')
plt.scatter(dot2[0],dot2[1],s=50,c='r')
plt.show()

从图中可以看出，10阶多项式模型只有在区间[-2π,2π]之间对正弦曲线拟合较好，在此区间以外，两者相差甚远。此案例告诉我们，每个模型都有自己的适用范围，在满足适用范围的基本前提下，要尽可能寻找拟合程度最高的模型来使用。

5.示例：预测房价

本节使用scikit-learn自带的波士顿房价数据来训练模型，然后用模型来预测房价。

1.输入特征

房价和哪些因素有关？很多人可能对这个问题特别敏感，随时可以列出很多，如房子面子、房子地理位置、周边教育资源、周边商业资源、房子朝向、年限、小区情况等。在scikit-learn的波士顿房价数据集里，它总共收集了13个特征，具体如下：

• CRIM：城镇人均犯罪率。
• ZN：城镇超过25000平方英尺的住宅区域的占地比例。
• INDUS：城镇非零售用地占地比例。
• CHAS：是否靠近河边，1为靠近，0为远离。
• NOX：一氧化氮浓度
• RM：每套房产的平均房间个数。
• AGE：在1940年之前就盖好，且业主自住的房子的比例。
• DIS：与波士顿市中心的距离。
• RAD：周边高速公路的便利性指数。
• TAX：每10000美元的财产税率。
• PTRATIO：小学老师的比例。
• B：城镇黑人的比例。
• LSTAT：地位较低的人口比例。
  从这些指标里可以看到中美指标的一些差异。当然，这个数据是在1993年之前收集的，可能和现在会有差异。不要小看了这些指标，实际上一个模型的好坏和输入特征的选择关系密切。大家可以思考一下，如果要在中国预测房价，你会收集哪些特征数据？这些特征数据的可获得性如何？收集成本多高？

我们先导入数据：

from sklearn.datasets import load_boston
boston = load_boston()
X = boston.data
y = boston.target
X.shape

输出如下：

(506, 13)

表明这个数据集有506个样本，每个样本有13个特征。整个训练样本放在一个506*13的矩阵里。可以通过X[0]来查看一个样本数据：

X[0]

输出如下：

array([6.320e-03, 1.800e+01, 2.310e+00, 0.000e+00, 5.380e-01, 6.575e+00,
       6.520e+01, 4.090e+00, 1.000e+00, 2.960e+02, 1.530e+01, 3.969e+02,
       4.980e+00])

还可以通过boston.features_names来查看这些特征的标签：

boston.feature_names

输出如下：

array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')

我们可以把特征和数值对应起来，观察一下数据。

2.模型训练

在scikit-learn里，LinearRegression类实现了线性回归算法。在对模型进行训练之前，我们需要先把数据集分成两份，以便评估算法的准确性。

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.2,random_state=3)

由于数据量比较小，我们只选了20%的样本来作为测试数据集。接着，训练模型并测试模型的准确性评分：

import time
from sklearn.linear_model import LinearRegression
model = LinearRegression()
start = time.process_time()
model.fit(X_train,y_train)
train_score = model.score(X_train,y_train)
cv_score = model.score(X_test,y_test)
print("elaspe:{0:.6f};train_score:{1:0.6f};cv_score:{2:.6f}"
      .format(time.process_time()-start,train_score,cv_score))

我们顺便统计了模型的训练时间，除此之外，统计模型对训练样本的准确性得分（即对训练样本拟合的好坏程度）train_score，还测试了模型对测试样本的得分sv_score。运行结果如下：

elaspe:0.000000;train_score:0.723941;cv_score:0.794958

从得分情况来看，模型的拟合效果一般，还有没有办法来优化模型的拟合效果呢？

3.模型优化

首先观察一下数据，特征数据的范围相差比较大，最小的在${10}^{-3}$级别，而最大的在${10}^2$级别，看来我们需要先把数据进行归一化处理。归一化处理最简单的方式是，创建线性回归模型时增加normalize=True参数：

model = LinearRegression(normalize=True)

当然，数据归一化处理只会加快算法收敛速度，优化算法训练的效率，无法提升算法的准确性。

怎么样优化模型的准确性呢？我们回到训练分数上来，可以观察到模型针对训练样本的评分比较低（train_score:0.723941），即模型对训练样本的拟合成本比较高，这是一个典型的欠拟合现象。回忆我们之前介绍的优化欠拟合模型的方法，一是挖掘更多的输入特征，而是增加多项式特征。在我们这个例子里，通过使用低成本的方案——即增加多项式特征来看能否优化模型的性能。增加多项式特征，其实就是增加模型的复杂度。

我们编写创建多项式模型的函数：

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
def polynomial_model(degress=1):
    polynomial_features = PolynomialFeatures(degree=degree,include_bias=False)
    linear_regression = LinearRegression(normalize=True)
    pipeline = Pipeline([("polynomial_features",polynomial_features),
                         ("linear_regression",linear_regression)])
    return pipeline

接着，我们使用二阶多项式来拟合数据：

model = polynomial_model(degree=2)
start = time.process_time()
model.fit(X_train,y_train)
train_score = model.score(X_train,y_train)
cv_score = model.score(X_test,y_test)
print("elaspe:{0:.6f};train_score:{1:0.6f};cv_score:{2:.6f}"
      .format(time.process_time()-start,train_score,cv_score))

输出结果是：

elaspe:0.437500;train_score:0.930547;cv_score:0.860465

训练样本分数和测试分数都提高了，看来模型确实得到了优化。我们可以把多项式改为3阶看一下效果：

elaspe:0.343750;train_score:1.000000;cv_score:-105.483692

改为3阶多项式后，针对训练样本的分数达到了1，而针对测试样本的分数确实负数，说明这个模型过拟合了。

思考：我们总共有13个输入特征，从一阶多项式变为二阶多项式，输入特征个数增加了多少个？
参考：二阶多项式共有：13个单一的特征，$C_{13}^2=78$个两两配对的特征，13个各自平方的特征，共计104个特征。比一阶多项式的13个特征增加了91个特征。（如果有其他答案，欢迎在评论区留言）

4.学习曲线

更好的方法是画出学习曲线，这样对模型的状态以及优化的方向就一目了然。

import matplotlib.pyplot as plt
from common.utils import plot_learning_curve
from sklearn.model_selection import ShuffleSplit
cv = ShuffleSplit(n_splits=10,test_size=0.2,random_state=0)
plt.figure(figsize=(18,4),dpi=200)
title = 'Learning Curves (degree={0})'
degrees = [1,2,3]
start = time.process_time()
for i in range(len(degrees)):
    plt.subplot(1,3,i+1)
    plot_learning_curve(plt,polynomial_model(degrees[i]),title.format(degrees[i]),
                        X,y,ylim=(0.01,1.01),cv=cv)
    print('elaspe:{0:.6f}'.format(time.process_time()-start))

其中，common.utils包里的plot_learning_curve()函数是对sklearn.model_selection.learning_curve()函数的封装，代码如下：

from sklearn.model_selection import learning_curve
import numpy as np

def plot_learning_curve(plt, estimator, title, X, y, ylim=None, cv=None,
                        n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):
    """
    Generate a simple plot of the test and training learning curve.

    Parameters
    ----------
    estimator : object type that implements the "fit" and "predict" methods
        An object of that type which is cloned for each validation.

    title : string
        Title for the chart.

    X : array-like, shape (n_samples, n_features)
        Training vector, where n_samples is the number of samples and
        n_features is the number of features.

    y : array-like, shape (n_samples) or (n_samples, n_features), optional
        Target relative to X for classification or regression;
        None for unsupervised learning.

    ylim : tuple, shape (ymin, ymax), optional
        Defines minimum and maximum yvalues plotted.

    cv : int, cross-validation generator or an iterable, optional
        Determines the cross-validation splitting strategy.
        Possible inputs for cv are:
          - None, to use the default 3-fold cross-validation,
          - integer, to specify the number of folds.
          - An object to be used as a cross-validation generator.
          - An iterable yielding train/test splits.

        For integer/None inputs, if ``y`` is binary or multiclass,
        :class:`StratifiedKFold` used. If the estimator is not a classifier
        or if ``y`` is neither binary nor multiclass, :class:`KFold` is used.

        Refer :ref:`User Guide <cross_validation>` for the various
        cross-validators that can be used here.

    n_jobs : integer, optional
        Number of jobs to run in parallel (default 1).
    """
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    train_sizes, train_scores, test_scores = learning_curve(
        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()

    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, 'o--', color="r",
             label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
             label="Cross-validation score")

    plt.legend(loc="best")
    return plt

输出的学习曲线如下图所示：

从学习曲线中可以看出，一阶多项式欠拟合，因为针对训练样本的分数比较低；而三阶多项式过拟合，因为针对训练样本的分数达到1，却看不到交叉验证数据集的分数。针对二阶多项式拟合的情况，虽然比一阶多项式的效果好，但从图中可以明显地看出来，针对训练数据集的分数和针对交叉验证数据集的分数之间的间隔比较大，这说明训练样本数量不够，我们应该去采集更多的数据，以提高模型的准确性。

6.拓展阅读

本节内容涉及到较多的数学知识，特别是矩阵和偏导数运算法则。如果阅读起来有困难，可以先跳过。如果有一定数学基础，这些知识对理解算法的实现细节及算法的效率有较大的帮助。

1.梯度下降迭代公式推导

关于梯度下降算法迭代公式的推导过程，可以参考博客：http://blog.kamidox.com/gradient-descent.html，或者直接搜索“线性回归算法kamidox.com”。博客里详细介绍了公式推导过程中用到的偏导数运算法则。

2.随机梯度下降算法

本章介绍的梯度下降算法迭代公式称为批量梯度下降算法（Batch Gradient Descent，简称BGD），用它对参数进行一次迭代运算，需要遍历所有的训练数据集。当训练数据集比较大时，其算法的效率会比较低。考虑另外一个算法：
$${\theta }_j={\theta }_j-\alpha ((h(x^{(i)})-y^{(i)})x_j^{(i)})$$
这个算法的关键点是把累加器去掉，不去遍历所有的数据集，而是改成每次随机地从训练数据集中取一个数据进行参数迭代计算，这就是随机梯度下降算法（Stochastic Gradient Descent，简称SGD）。随机梯度下降算法可以大大提高模型训练的效率。

思考：为什么随机取一个样本进行参数迭代是可行的？
从数学上证明批量梯度下降算法和随机梯度下降算法的等价性涉及到复杂的数学知识。这里有个直观的解释可以帮助理解两者的等价性。回到成本函数的定义：
$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h(x^{(i)})-y^{(i)}{)}^2$$
我们说过，这里累加后除以2是为了计算方便，那么除以m是什么意思呢？答案是平均值，即所有训练数据集上的点到预测函数的距离的平均值。再回到随机选取训练数据集里的一个数据这个做法来看，如果计算次数足够多，并且是真正随机，那么随机选取出来的这组数据从概率的角度来看，和平均值是相当的。打个比方，储钱罐里有1角的硬币10枚，5角的硬币2枚，1元的硬币1枚，总计3元、13枚硬币。随机从里面取1000次，把每次取出来的硬币币值记录下来，然后将硬币放回储钱罐里。这样最后去算这1000次取出来的钱的平均值（1000次取出来的币值总和除以1000）和储钱罐里每枚硬币的平均值（3/13元）应该是近似相等的。

3.标准方程

梯度下降算法通过不断地迭代，从而不停地逼近成本函数的最小值来求解模型的参数。另外一个方法是直接计算成本函数的微分，令微分算子为0，求解这个方程，即可得到线性回归的解。
线性回归算法的成本函数：
$$J(\theta )=\frac{1}{2m}\sum _{i=0}^n(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$$
成本函数的“斜率”为0的点，即为模型参数的解。令$\frac{\partial }{\partial \theta }J(\theta )=0$，求解这个方程最终可以得到模型参数：
$$\theta =(X^{T}X{)}^{-1}{X}^{T}y$$
方程求解过程可参阅https://en.wikipedia.org/wiki/Linear_least_squares_（mathematics）#Derivation_of_the_normal_equations。
这就是我们的标准方程。它通过矩阵运算，直接从训练样本里求出参数θ的值。其中X为训练样本的矩阵形式，它是m×n的矩阵，y是训练样本的结果数据，它是个m维列向量。