机器学习理论及案例分析(part2)--回归

本节全部代码，已上传至https://github.com/GoatBishop/IntroductionToDataScienceC7

博客阅读索引：博客阅读及知识获取指南

下一节：聚类

---

文章目录

---

回归

回归是对一个或多个自变量和因变量之间的关系进行建模，求解的一种统计方法。早在大学的统计学课本上，我们就学习过回归，回归模型看似简单，但确是最重要的数学模型之一，很多模型都是在它的基础上建立的，任何一个复杂模型，其内部可能隐藏着许许多多回归模型。

下面，我们通过一个一元线性回归的例子引出线性回归模型：

举个例子

我们从R自带数据集中获取到了1920年汽车速度与刹车距的数据，现在，我们想研究速度与刹车距离之间到底有什么样的关系。

#section/example01.py
#读取数据
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

data_cars = pd.read_csv("../data/cars.csv", 
                        usecols = ["speed", "dist"])

表1展示了车速-刹车距离的部分数据，其中speed为汽车速度，$dist$为刹车距离：

表1 车速-刹车距离的部分数据

	speed	dist
1	4	2
2	4	10
3	7	4
4	7	22
5	8	16
6	9	10
7	10	18
8	10	26
9	10	34
10	11	17

首先，我们想对这个数据集有一个直观的认识，所以我们先对其进行可视化，即绘制车速与刹车距离散点图，如图1所示。

#section/example01.py
#绘制散点图
speed = data_cars["speed"]
dist = data_cars["dist"]

plt.scatter(speed, dist)
plt.title("Scatter plot of vehicle speed and braking distance")
plt.show()

图1 车速与刹车距离散点图

观察上图，我们可以看到，车速和刹车距离之间貌似呈现出一种线性相关关系，数据点在一条看不见的直线上下随机散落，因而，我们可以考虑用下面的公式来描述他们之间的关系：

$y = \alpha + \beta x \tag{7.1.1}$
其中，$y$为刹车距离，$x$为车速，$\alpha$和$\beta$为回归系数。

我们建立模型的目的，就是为了进行预测，为了使预测更为精准，我们就需要使预测结果与真实结果之间的距离越小，如果我们使用二范数的平方和的$1/2$来度量这里的距离，那么用数学语言描述，就是使如下公式达到最小：
$L(\alpha, \beta) = \frac{1}{2} \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 = \frac{1}{2} \sum_{i=1}^n ||y_i - (\alpha + \beta x_i)||_2^2 \tag{7.1.2}$

其中，n为样本数，$\hat{y}_i$为第$i$样品刹车距离的预测值，$y_i$为第$i$样品刹车距离的真实值，$x_i$为第$i$样品的车速。

从图像直观描述，就是使寻找一条直线，使得所有样本点到直线的距离之和最短，如图2所示：

图2 样本点到直线最短距离合示范

在机器学习中，我们称(7.2)式为损失函数，我们可以通过直接法和迭代法两种方式对该损失函数进行优化，进而得到使损失函数最小的回归系数。下面，我们分别用这两种方法进行求解。

• 直接法

直接法，就是直接给出优化问题的最优解，并不是所有的优化问题都可以用直接法得到最优解，如果要使用直接法，损失函数需要满足两个条件：

1. 损失函数为凸函数；
2. 损失函数为解析解，即通过严格的公式所求得的解。

凸函数定义(补)：

设函数$f(x)$为凸函数，当且仅当对定义域中任意两点$x_1$,$x_2$和任意实数$\lambda \in [0, 1]$，总有：
$f(\lambda x_1 + (1-\lambda) x_2) \le \lambda f(x_1) + (1-\lambda) f(x_2)$
通过图像直观来看，凸函数上任意两点连接而成的虚线，永远在凸函数曲线的上方，如图3所示：

图3 凸函数曲线图

由上述定义，我们可以了解到，如果一个函数是凸函数那么我们可以求解出它的全局最优解，但如果一个函数是非凸的，那么我们可能求出的是它的局部最优解，如图4所示：

图4 非凸函数曲线图

由上图可知，如果函数是非凸，那么最终我们得到的可能是局部最优点$x_a$，而不是全局最优点$x_k$.

很明显，我们的损失函数$L(\alpha, \beta)$满足以上两个条件，现在，我们直接对各个回归系数求偏导，并令其等于0：
$\frac{\partial{L}}{\partial{\alpha}} = -\sum_{i=1}^n (y_i - (\alpha + \beta x_i)) = 0 \tag{7.1.3}$

$\frac{\partial{L}}{\partial{\beta}} = -\sum_{i=1}^n x_i (y_i - (\alpha + \beta x_i)) = 0 \tag{7.1.4}$

对式(7.1.3)和式(7.1.4)加以推导，我们可以得到最优解：
$\beta = \frac{n \sum_{i=1}^n x_i y_i - \sum_{i=1}^n x_i \sum_{i=1}^n y_i}{n\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2} \\ \alpha = \overline{y} - \beta \overline{x} \tag{7.1.5}$
现在，我们利用python对其进行求解。

#section/example01.py

import sympy
#设定回归系数
alpha, beta = sympy.symbols("alpha beta")

#设定损失函数
L = 0.5*np.sum((dist - beta*speed - alpha)**2)

#求偏导
print(sympy.diff(L, alpha))
#50.0*alpha + 770.0*beta - 2149.0
print(sympy.diff(L, beta))
#770.0*alpha + 13228.0*beta - 38482.0
f1 = sympy.diff(L, alpha)
f2 = sympy.diff(L, beta)

#求解线性方程组
outcome = sympy.solve([f1, f2], [alpha, beta])
print(outcome)
#{alpha: -17.5790948905109, beta: 3.93240875912409}

我们求出$\alpha$和$\beta$的最优解分别为-17.579和3.932。

我们将直线$y = 3.932 x + -17.579$，绘制在车速与刹车距离散点图中，如图5所示。

#section/example01.py

#绘图
alpha_num = outcome[alpha]
beta_num = outcome[beta]
#得到预测值
dist_pre = beta_num*speed + alpha_num
plt.scatter(speed, dist)
plt.plot(speed, dist_pre, c = "r")
plt.title("Fitting results")
plt.show()

图5 拟合图

从图像上来看，我们较好的拟合的原始数据。

• 迭代法

我们知道，并不是所有的损失函数都可以使用直接法对损失函数进行优化的，很多情况中，我们的损失函数为非凸函数，或者方程过于复杂难以求出其解析解，因此，在有些实际问题中，我们可以选择使用迭代法求解。迭代法是一种不断用变量的旧值递推新值的过程，即迭代的用旧值修正对最优解的估计。这里，我们使用迭代法中的小批量梯度下降法解决问题。

对于一般问题，小批量梯度下降法的过程如下：

确定要求解的模型参数,模型参数可以1个，也可以是多个；
定义损失函数L，确定迭代递推关系，用于更新模型参数；
确定模型参数初始值,学习率,迭代次数,以及批量大小。

count = 1
while (count <= 迭代次数):
    count += 1

对于我们的问题，小批量梯度下降法的过程如下：

(1)确定要求解的模型参数为$\alpha$，$\beta$；

(2)定义小批量梯度下降法的损失函数：
$L(\alpha, \beta) = \frac{1}{2m} \sum_{i=1}^m ||y_i - (\alpha + \beta x_i)||_2^2 \tag{7.1.6}$
其中,$m$为批量大小，即每次迭代从原始数据集中取多少数据代入损失函数的计算，一般$m$取2的幂次，这里我们定义批量大小为16；

(3)求解梯度，并定义递推关系：
$\alpha_{t+1} = \alpha_t - \gamma_1 \nabla_{\alpha} L(\alpha, \beta) = \alpha_t -(-\frac{1}{m} \sum_{i=1}^m (y_i - (\alpha + \beta x_i))) \tag{7.1.7}$

$\beta_{t+1} = \beta_t - \gamma_2 \nabla_{\beta} L(\alpha, \beta) = \beta_t -(-\frac{1}{m} \sum_{i=1}^m x_i (y_i - (\alpha + \beta x_i))) \tag{7.1.8}$

在[0, 10)中按照均匀分布随机产生alpha和beta的初始值，定义学习率$\gamma_1$和$\gamma_2$均为0.02，定义迭代次数为20000次；

(4)迭代，迭代完成输出最后的模型参数。

我们用python求解过程如下：

#迭代法
import random

#定义递推关系,更新迭代变量
def update_var(old_alpha, old_beta, y, x, learning_rate):
    len_x = len(x)
    alpha_delta = np.sum(-(y - old_beta*x - old_alpha))/len_x
    beta_delta = np.sum(-x*(y - old_beta*x - old_alpha))/len_x
    new_alpha = old_alpha - learning_rate*alpha_delta
    new_beta = old_beta - learning_rate*beta_delta
    
    return (new_alpha, new_beta)

#迭代
def iterative_func(y, x, start_alpha, start_beta,
                   learning_rate, iterative_num, 
                   sample_num):
    alpha_list = []
    beta_list = []
    alpha = start_alpha
    beta = start_beta
    num_list = list(range(1, len(y)))
    for i in range(iterative_num):
        alpha_list.append(alpha)
        beta_list.append(beta)
        random.shuffle(num_list)
        
        index = num_list[:sample_num]
        alpha, beta = update_var(alpha, beta, 
                                 y[index], x[index], learning_rate)
    return (alpha_list, beta_list)

#在[0, 10)之间按照均匀分布随机产生alpha和beta的初始值
start_alpha = np.random.random()*10
start_beta = np.random.random()*10

#设置学习率为0.002,迭代次数为2000次,批量大小为16
learning_rate = 0.002
iterative_num = 20000
sample_num = 16

alpha_list, beta_list = iterative_func(dist, speed, start_alpha, start_beta,
                                     learning_rate, iterative_num,
                                     sample_num)

print("alpha: {}, beta:{}".format(alpha_list[-1], beta_list[-1]))
#alpha: -17.833867412345022, beta:3.9820296889683062

可以看到，我们迭代20000次后，$\alpha=-17.834$，$\beta = 3.982$，我们将迭代结果输出，并将$\alpha$和$\beta$的迭代过程绘制出来，如图6所示。

import csv
parameter_data = zip(alpha_list, beta_list)
with open("./outcome/gradient_descent_parameter.csv", 'w', newline = '') as f:
    csv_writer = csv.writer(f)
    csv_writer.writerow(["alpha","beta"])
    csv_writer.writerows(parameter_data)
    
plt.subplot(121)
plt.plot(alpha_list)
plt.title("alpha change process")
plt.subplot(122)
plt.plot(beta_list)
plt.title("beta change process")
plt.show()

图6 迭代图

我们看到$\alpha$在迭代20000次后参数逐渐趋于平稳，可以预想到如果再迭代10000次，那么图像末端会近乎于水平延伸；$\beta$在15000次后，貌似在一条直线上下波动，且无法稳定，这可能是由于我们针对$\beta$的学习率设置过大，导致模型参数估计值在最优值附近来回震荡，此时，我们可以设置其学习率随着迭代次数增加逐渐递减，来缓解这种现象。

尽管如此，我们通过迭代法得到的模型参数估计值，依然是不错的，其值与直接法计算的模型参数估计值非常接近。

我们得到模型参数(回归系数)的估计值后，需要对模型进行评价，在统计学上，我们通常用判定系数$R^2$说明回归直线对数据拟合的好坏，其数学表达为：
$R^2 = 1- \frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{\sum_{i=1}^n (y_i - \overline{y}_i)^2} \tag{7.1.9}$
判定系数$R^2$测度了回归直线对观测数据的拟合程度，$R^2$越接近1，则说明各个数据点与回归直线越接近，回归模型的预测值越接近真实值；$R^2$越小，则说明各个数据点离回归直线越远，回归模型的预测值越偏离真实值，模型拟合效果越差。

现在，我们对用直接法求得的回归系数估计值进行评价，即计算其$R^2$.

dist_pre = beta_num*speed + alpha_num
dist_mean = np.mean(dist)
R_2 = np.sum((dist_pre - dist_mean)**2)/np.sum((dist - dist_mean)**2)
print(R_2)
#0.651079380758251

由结果可知，判定系数$R^2=0.651$，看起来模型拟合的还算不错。

现在，我们可以利用我们的模型对结果进行预测了，我们依然采用直接法求得的回归系数估计值，对新的样本进行预测。假设新记录的一批车辆的车速分别为10, 20, 30，利用python，我们可以得到刹车距离的预测值为21.745，61.069，100.393。

new_speed = pd.Series([10, 20, 30])
new_dist_pre = beta_num*new_speed + alpha_num
print(new_dist_pre)
#0    21.7449927007299
#1    61.0690802919708
#2    100.393167883212
#dtype: object

通过这个例子，我们已经大致了解了一元线性回归模型的建立，参数估计，评价以及预测，现在，我们对这些步骤及概念进行拓展，并展开一般性的讨论。

回归的基本思想

多元回归基本原理

在上面的例子中，我们仅有一个自变量，但是在现实生活中，我们要分析的不只是一个，而是多个自变量与因变量之间的关系，对于这种预测任务，我们可以采用多元线性回归，假设我们有k个自变量，我们用下面的公式来描述自变量与因变量之间的关系：
$y=\beta_0 + \beta_1 x_1 + \cdot \cdot \cdot + \beta_k x_k \tag{7.1.10}$

如果我们有一组样本量为n的样本$\{ (x_{i1}, x_{i2}, \cdots , x_{ik}, y_i):i = 1,2, \cdots,n \}$，我们$y_i$与$( x_{i1}, x_{i2}, \cdots , x_{ik})$的对应关系写成矩阵形式，即：
$Y= X \beta \tag{7.1.11}$
其中：
$Y=\left( \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{matrix} \right), \beta = \left( \begin{matrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{matrix} \right) \\ X = \left( \begin{matrix} 1 & x_{11} & x_{12} & ... & x_{1k} \\ 1 & x_{21} & x_{22} & ... & x_{2k} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & ... & x_{nk} \end{matrix} \right)$

我们的目的是让机器从数据中学习出一个最优模型，从而进行预测，但是机器不会自动的学习，它需要被注入灵魂，需要在人的指导下去工作。所以，我们需要告诉机器如何学习，即定义一系列"策略"让机器围绕着这种"策略"不断地学习，更新自己。比如，在多元线性回归的问题中，我们的"策略"如下：

• 定义损失函数，描述模型好坏
• 利用算法最小化损失函数

损失函数

与一元线性回归类似，在多元线性回归中，我们将损失函数设定为二范数的平方和，即：
$L(\beta_0, \beta_1, \cdots, \beta_k) = \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 \\ = \sum_{i=1}^n ||y_i - (\beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik})||_2^2 \tag{7.1.12}$
由公式(7.1.12)可知，只要损失函数$L(\beta_0, \beta_1, \cdots, \beta_k)$达到最小，那么预测值与真实值之间的距离就是最小的，此时，模型对训练集的拟合效果是最好的。

自然而然，我们会寻求一种算法，使得损失函数达到最小，即：
$min_{\beta} \; L(\beta_0, \beta_1, \cdots, \beta_k) = min_{\beta} \; \sum_{i=1}^n ||y_i - (\beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik})||_2^2$

最小二乘法

对于回归问题，我们可以直接利用最小二乘法得到最优解，对于式(7.1.12)我们也可以写成矩阵形式：
$L(\beta) = (Y - X \beta)^T (Y - X \beta) \tag{7.1.13}$
对$\beta$求偏导，即：
$\frac{\partial L(\beta)}{\partial \beta} = -2X^T(Y-X\beta) \tag{7.1.14}$
当$X^TX$为满秩矩阵时，令式(7.14)其为0，并求解得$\beta$为：
$\beta = (X^T X)^{-1} X^TY \tag{7.1.15}$

需要注意的是，有时自变量之间会存在不完全多重共现性(即自变量之间存在高度相关，但是不存在完全共线性)，从而使$(X^TX)^{-1}$方差过大；有时$X^TX$不是满秩矩阵(此时无法用式(7.1.15)计算$\beta$)，导致$\beta$有多个解。此时我们可以使用梯度下降法，或者添加正则项的方式解决这些问题。

梯度下降法

对于上述找到一组$\beta$使得损失函数$L(\beta_0, \beta_1, \cdots, \beta_k)$最小化的问题，也叫优化问题，在优化问题求解时，我们经常会使用泰勒展开来近似代替目标函数，梯度下降法就是利用一阶泰勒展开，使$\beta$一步步逼近最优，从而最小化损失函数的方法。

在梯度下降法中，我们重新设定损失函数$L(\beta)$为：
$L(\beta_0, \beta_1, \cdots, \beta_k) = \frac{1}{2n} \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 \\ = \frac{1}{2n} \sum_{i=1}^n ||y_i - (\beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik})||_2^2 \tag{7.1.16}$

损失函数$L(\beta)$在$\beta$的一阶泰勒展开式为：
$L(\beta + \Delta \beta) = L(\beta) + (\nabla L(\beta))^T\Delta \beta + o(\Delta \beta) \tag{7.1.17}$
对式(7.1.17)进行变换：
$L(\beta + \Delta \beta) - L(\beta) = (\nabla L(\beta))^T\Delta \beta + o(\Delta \beta) \tag{7.1.18}$
若有：
$(\nabla L(\beta))^T\Delta \beta < 0 \tag{7.1.19}$

则有：
$L(\beta + \Delta \beta) - L(\beta) < 0 \tag{7.1.20}$
损失函数递减，式(7.1.19)和式(7.1.20)说明，当增量$\Delta \beta$选择适当时，损失函数就可以递减。可证明，当$\Delta \beta = -\nabla L(\beta)$时，损失函数下降最快，即增量$\Delta \beta$选择与梯度相反方向时，损失函数下降最快。$\nabla L(\beta)$的定义为：
$\nabla L(\beta) = \left( \frac{\partial{L}}{\partial{\beta_0}},\frac{\partial{L}}{\partial{\beta_1}},\cdots,\frac{\partial{L}}{\partial{\beta_k}} \right) \tag{7.1.21}$
可设置一个接近0的正数$\gamma$，控制机器学习的速度，我们也称其为学习率，则增量为：
$\Delta \beta = - \gamma \nabla L(\beta)$
则：
$L(\beta + \Delta \beta) - L(\beta) = (\nabla L(\beta))^T\Delta \beta = - \gamma (\nabla L(\beta))^T \nabla L(\beta) \le 0$
由随机初始值$\beta_{(0)}$开始，$\beta$的迭代公式如下：
$\beta_{(t+1)} = \beta_{(t)} - \gamma \nabla L(\beta_{(t)}) \tag{7.1.22}$
我们可以设置终止条件为梯度$\nabla L(\beta_{(t+1)})$小于某个接近0的阈值，当迭代停止时，若损失函数为凸函数，则认为此时取到全局最优解，若为非凸函数，则不一定能取到全局最优解，可能最终得到的是局部最优解。因为我们优化的是关于多元线性回归模型的损失函数，所以函数是凸的，因此，我们可以不必考虑局部最优和全局最优不一致的问题。

算法实现具体步骤:

输入数据，将数据划分为训练集和测试集
定义损失函数式(7.12)
产生参数初始的beta0, beta1,...,betak
设定学习率
通过式(7.13)计算各个参数的梯度,设置阈值
while (梯度 > 阈值):
    在训练集上选择样品训练模型，利用迭代公式(7.14)更新参数

得到满足要求的模型，并将测试集数据代入。

在梯度下降法中，由于每次迭代带入损失函数中样品个数的不同，我们又可将其分为批量梯度下降法(BGD)，小批量梯度下降法(MBGD)和随机梯度下降法(SGD)。

具体来说，批量梯度下降法在每一步迭代中，会将训练集中所有样品带入损失函数，这样每一步的迭代都会消耗大量的时间，算法的求解过程会非常漫长。为了解决这个问题，有些人就考虑，我们是否可以在每次迭代中，不使用全部训练集数据，而每次只选取一部分样品，甚至只选取一个样品来计算损失函数。使用一部分样品的方法就是我们在本节开头的例子中使用的小批量梯度下降法，而只选取一个样品的方法，就是随机梯度下降法，在sklearn中我们可以调用linear_model.SGDRegressor类去使用它。

辛烷值预测案例

某石化企业的催化裂化汽油精制脱硫装置运行4年，积累了大量历史数据。从催化裂化汽油精制装置中，我们得到了325个数据样本，因变量为产品中的辛烷值(RON)，自变量有310个，它们分别为生产前的原材料含量和各种操作变量。这些数据含有少量缺失值，我们已经利用缺失值填补技术，将其填补好了，在这里就不加赘述。

在特征量非常大，甚至快要逼近样本量的情况下，我们建立基于最小二乘法的多元线性回归模型明显是有问题的，但是为了说明某些问题，我们就先这样建模。

首先，导入sklrean中的model_selection.train_test_split划分训练集和测试集，并利用sklearn的linear_model.LinearRegression类实现基于最小二乘法的多元回归，最后利用score方法得到模型得分，score方法计算了$R^2$，即判定系数，其最好的可能得分是1.0，但也可以因为模型过差得到负数。

#section1/example02.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
#读取数据
data_ron = pd.read_csv("../data/data_ron.csv")
print(data_ron.shape)
#(325, 311)

#划分训练集和测试集
y = data_ron["RON"]
X = data_ron.iloc[:, 1: (len(data_ron) + 1)]

X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size = 0.2, random_state = 10)

#构建模型
model = LinearRegression()
model.fit(X_train, y_train)
#模型得分
print(round(model.score(X_train, y_train), 4))
#1.0

可以看到，我们在训练集上的模型得分居然是1.0，这是不是说明我们的模型构建的足够好呢？不急，我们再看看模型在测试集上的得分。

print(round(model.score(X_test, y_test), 4))
#-68.3224

由结果可知，模型在测试集上的得分为负数，已经足够差了。

这种在训练集上表现好，但在测试集上表现差的现象，我们称为过拟合，与其相对的就是在训练集上表现差，但在测试集上表现好，这种我们称之为欠拟合，下面，我们将对其进行详细讲解。

过拟合与欠拟合

泛化

我们在上面的案例中，通过训练数据构建了多元回归模型，并对模型没有见过的测试集进行预测(在计算得分时python内部进行了预测，score方法的具体计算过程，大家可以查看sklearn官方文档，或者源码)。如果我们的回归模型能够对测试集进行准确预测，那么我们就说它具有较好的泛化能力。

过拟合与欠拟合

一般来说，我们的模型在训练集上表现会优于测试集，但是如果我们过度追求在训练集上的准确性，而将模型构建的过于复杂，那么就会出现过拟合的情况。这时，我们的模型会过分关注于训练集的细节，从而在训练集上得到非常优秀的预测结果，但是这种模型无法泛化到新的数据上，一旦将新的样本带入模型，那么模型的表现会非常的差。

在上面的案例中，由于我们使用过多的因变量，而导致训练集得分很高，测试集得分很低，这就是出现了过拟合的情况。可以说，当回归模型的特征数量足够大，甚至大于样本量时，它可以对任意训练集中的因变量进行完美的预测。

遇到过拟合时，我们可以通过增大样本量的方法，减缓过拟合的状况，当我们的样本量足够大时，数据就会存在更多的变化，在不发生过拟合的情况下，我们可以构造出比数据缺乏时更复杂的模型，从而做出更精准的预测。同时，我们也可以采取降低模型复杂度的方式来避免过拟合的发生。

当我们的模型过于简单时，可能会出现欠拟合的情况，这说明了我们的训练不够充分，模型没有抓住训练集中数据的信息。如果在上面的案例中，我们的模型出现了欠拟合，那么在训练集上，迫性的得分应该较低，而测试集上的模型得分应该较高。

遇到欠拟合时，我们应该适量增加模型的复杂度，来避免欠拟合的状况。

在日常生活中，我们应该注意权衡模型复杂度，找到一个最佳平衡点，即当复杂度小于平衡点时，可能会出现欠拟合，当复杂度大于平衡点时，可能会出现过拟合，如图7所示：

图7 复杂度-误差图

简化回归模型

我们对辛烷值预测案例中的模型进行简化，通过减少自变量的方式降低模型的复杂度。通过特征工程技术，我们将因变量筛选到29个，降维的具体过程不进行赘述。

现在，我们对降维后的数据建立多元回归模型，具体步骤还是和之前一样，只是数据的由$325 \times 349$的矩阵，变成了$325 \times 30$的矩阵，在这里我们只贴模型得分部分的代码和结果。

#section1/example02_2.py

print(round(model.score(X_train, y_train), 4))
#0.9619
print(round(model.score(X_test, y_test), 4))
#0.9588

可以看到，我们的模型在训练集和测试集上表现都很不错，过拟合的情况消失了。

模型的评价

我们在上面的例子中已经介绍了一种评价模型的指标，即判定系数$R^2$，现在，我们再介绍几种模型评价的指标。

MSE和RMSE

在回归分析中，最常用的评价模型的指标就是均方差$MSE$以及均方根误差$RMSE$，由名字我们就可以了解到$RMSE$就是$MSE$开根号后的结果，它们的定义如下：
$MSE = \frac{1}{n} \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 \tag{7.1.23}$

$RMSE = \sqrt{MSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n}(y_i - \hat{y}_i)^2} \tag{7.1.24}$

MAE和MAPE

平均绝对误差$MAE$是将预测误差取绝对值后的平均误差，平均百分比误差$MAPE$则消除了因变量单位的影响，反映了误差大小的相对值，它们的定义如下：
$MAE = \frac{\sum_{i=1}^{n}|y_i - \hat{y}_i|}{n} \tag{7.1.25}$

$MAPE = \frac{\sum_{i=1}^{n} \left( |\frac{y_i - \hat{y}_i}{y_i}| \times 100 \right) }{n} \tag{7.1.26}$

备注：MAPE的公式在贾俊平统计学中的绝对值不包含分母，这里的MAE在贾俊平统计学中是MAD。

偏差-方差权衡(仔细研究一下)

在回归分析中，我们假定误差项$\epsilon_i=y_i - \hat{y}_i$满足以下条件：(原文是残差项之间相互独立)

• 均值为0，即$E(\epsilon_i) = 0$
• 同方差，即$Var(\epsilon_i) = \sigma^2$
• 任意两个误差项之间不相关，即$Cov(\epsilon_i, \epsilon_j)=0, i \not= j$

那么模型的均方误差$MSE$的期望可以分解为：
$E(MSE)^ = (irreducible \; error \; variance) + (model \; bias)^2 + (model \; variance) \tag{7.1.27}$
由式(7.1.27)可以看出，测试均方误差的期望被分解成了3个部分，不可约误差的方差，模型偏差的平方，模型方差。让我们来一项一项的解读这3个术语。

我们构造回归模型，得到预测值$\hat{y}$，$\hat{y}$的精确程度可取决于两个量，一个是可约误差，一个是不可约误差，即：

模型误差 = 可约误差 + 不可约误差

当我们选取的模型$\hat{f}$不是最优模型$f$的最优估计时，这种因为选取问题而产生的一部分误差，就是可约误差，如果提高$\hat{f}$的精度，选取更适合的模型，这种误差就能被降低，因此这种误差是可约的。但是，当我们已经得到最优模型$f$的最优估计时，预测依然存在误差，那么这部分误差就是不可约的，无论我们的模型多么精准，我们都不能减少不可约误差。

假定$\hat{f}$和$X$是固定的，易证明：
$E(Y - \hat{Y})^2 = E[f(X) + \epsilon - \hat{f}(X)]^2 = [f(X) - \hat{f}(X)]^2 + Var(\epsilon) \tag{7.1.28}$
其中，$Y$为真实值，$X$为一组因变量，$\hat{Y}$为预测值，$\hat{f}$为对最优模型的估计模型，$f$为最优模型，$\epsilon$为误差项.

式(7.1.28)中$[f(X) - \hat{f}(X)]^2$为可约误差的平方，$Var(\epsilon)$为不可约误差的方差。

由上面的信息可知，式(7.1.27)也可以写为：
$E(MSE) = Var(\epsilon)+ (bias(\hat{f}(X)))^2 + Var(\hat{f}(X)) \tag{7.1.29}$

第二项模型偏差的平方反应了当前选取的模型的函数形式，能够在多大程度上逼近真实数据。如图8所示，y与x之间实际存在着非线性的关系，但是我却用线性回归模型去拟合它，必然会出现较大的偏差：

图8 线性模型拟合非线性关系数据图

此时，我们无论用多少训练数据来训练该线性回归模型，都会存在较大偏差：

图9 高偏差低方差图

可以看到，无论训练数据怎样改变，我们的线性回归模型依然有着较高的偏差，同时，从图像上来看模型的位置基本没有发生变化，也就是说我们的模型基本没有发生改变，即模型具有低方差性。

第三项模型的方差反映了用不同训练集估计最佳模型$f$时，估计模型的改变程度。如果用不同训练集训练数据得到的估计模型$\hat{f}$变化不大，则可以说模型方差较小，反之，如果训练集仅有微小变化，就能很大程度的改变估计模型$\hat{f}$，则可以说模型有较大方差。如图10所示，虽然我们的模型具有低偏差，但是模型随着训练集的改变会发生较明显的改变，说明其有高方差：

图10 低偏差高方差图

理想情况下，我们希望得到一个低方差且低偏差的模型，但是这是很难实现的，如果模型具有低方差，那么它很可能因为比较简单而具有高偏差，一个高偏差的模型可能出现欠拟合的状况。如果模型具有低偏差，那么他可能因为比较复杂而具有高方差，一个高方差的模型则可能会有过拟合的状况。所以，在建模的时候，我们需要控制模型整体的偏差和方差，选取一个合适的模型，避免出现欠拟合和过拟合。

评价回归模型

现在，利用刚才学习的指标，我们对案例中建立的回归模型进行评价。我们使用sklearn中的metrics计算MSE，RMSE，MAE，自编函数计算MAPE。再利用mlxtend中的bias_variance_decomp方法计算模型的方差和偏差。

#MSE/RMSE/MAE/MAPE
from sklearn import metrics
def mape(y_true, y_pre):
    n = len(y_true)
    mape = (sum(np.abs((y_true - y_pre)/y_true))/n)*100
    return mape

y_hat = model.predict(X_test)
MSE = metrics.mean_squared_error(y_test, y_hat)
RMSE = metrics.mean_squared_error(y_test, y_hat)**0.5
MAE = metrics.mean_absolute_error(y_test, y_hat)
MAPE = mape(y_test, y_hat)


print("MSE:{:.4f}, RMSE:{:.4f}, MAE:{:.4f}, MAPE:{:.4f}".format(MSE, RMSE, MAE, MAPE))
#MSE:0.0370, RMSE:0.1924, MAE:0.1511, MAPE:0.1710

#偏差和方差
from mlxtend.evaluate import bias_variance_decomp
#不进行转换会报错
X_train = np.array(X_train)
y_train = np.array(y_train)
X_test = np.array(X_test)
y_test = np.array(y_test)

mse, bias, var = bias_variance_decomp(model, X_train, y_train, X_test, y_test,loss='mse', 
                                     num_rounds=30, random_seed=1)

print("mse:{:.4f}, bias:{:.4f}, var:{:.4f}".format(mse, bias, var))
#mse:0.3947, bias:0.3443, var:0.0504

由结果可知，测试集中MSE为0.0370, RMSE为0.1924, MAE为0.1511, MAPE为0.1710，均方误差的期望为0.3947，模型的偏差为0.3443，模型的方差为0.0504。

bias_variance_decomp.py源码如下(可删除)：

import numpy as np

def _draw_bootstrap_sample(rng, X, y):
    sample_indices = np.arange(X.shape[0])
    bootstrap_indices = rng.choice(sample_indices,
                                   size=sample_indices.shape[0],
                                   replace=True)
    return X[bootstrap_indices], y[bootstrap_indices]


def bias_variance_decomp(estimator, X_train, y_train, X_test, y_test,
                         loss='0-1_loss', num_rounds=200, random_seed=None):

    supported = ['0-1_loss', 'mse']
    if loss not in supported:
        raise NotImplementedError('loss must be one of the following: %s' %
                                  supported)

    rng = np.random.RandomState(random_seed)

    all_pred = np.zeros((num_rounds, y_test.shape[0]), dtype=np.int)

    for i in range(num_rounds):
        X_boot, y_boot = _draw_bootstrap_sample(rng, X_train, y_train)
        pred = estimator.fit(X_boot, y_boot).predict(X_test)
        all_pred[i] = pred

    if loss == '0-1_loss':
        main_predictions = np.apply_along_axis(lambda x:
                                               np.argmax(np.bincount(x)),
                                               axis=0,
                                               arr=all_pred)

        avg_expected_loss = np.apply_along_axis(lambda x:
                                                (x != y_test).mean(),
                                                axis=1,
                                                arr=all_pred).mean()

        avg_bias = np.sum(main_predictions != y_test) / y_test.size

        var = np.zeros(pred.shape)

        for pred in all_pred:
            var += (pred != main_predictions).astype(np.int)
        var /= num_rounds

        avg_var = var.sum()/y_test.shape[0]

    else:
        avg_expected_loss = np.apply_along_axis(
            lambda x:
            ((x - y_test)**2).mean(),
            axis=1,
            arr=all_pred).mean()

        main_predictions = np.mean(all_pred, axis=0)

        avg_bias = np.sum((main_predictions - y_test)**2) / y_test.size
        avg_var = np.sum((main_predictions - all_pred)**2) / all_pred.size

    return avg_expected_loss, avg_bias,

正则项

在建模时，我们总是希望自己的模型能够尽量在训练集上取得较高的精度，又希望模型有好的泛化能力，这时，我们可以通过给损失函数添加正则项实现这一目标。正则项的使用可以减缓过拟合的状况，也可以帮助我们选择模型的特征。

岭回归

岭回归通过在基本线性回归损失函数的基础上，添加$L_2$正则化项，来避免过拟合的发生，岭回归损失函数为：
$L_{L_2} = L(\beta) + R_{L_2} = \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 + \lambda \sum_{j=1}^{k} ||\beta_j||_2^2 \tag{7.1.30}$
其中$\lambda$为正则化项前的调节因子，$k$为回归系数的个数。

岭回归中，由于$L_2$正则项的特性，模型会保留所有的系数，但是回归系数会被限制在一定范围内，从而减小特征对结果的影响，即以牺牲偏差来换取更小方差。

当我们选取的调节因子$\lambda$很大时，回归系数会接近于0，正则项过分约束模型。当选取的调节因子$\lambda$很小时，回归系数基本不变，正则项对模型基本没有产生约束。

Lasso回归

Lasso回归在基本线性回归损失函数的基础上添加$L_1$正则化项，由于$L_1$正则项的特性

，Lasso可以使某些特征的系数为0，从而起到特征选择的作用，Lasso回归的损失函数为：
$L_{L_1} = \frac{1}{2n} L(\beta) + R_{L_1} = \frac{1}{2n} \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 + \lambda \sum_{j=1}^{k} ||\beta_j||_1 \tag{7.1.31}$
需要注意的是，由于Lasso回归的损失函数中含有绝对值函数，所以我们不能对其进行二次求导，即无法使用牛顿法等优化方法对其进行求解，需要利用坐标下坡或者近端梯度法进行求解。

ElasticNet回归

ElasticNe回归在基本线性回归损失函数的基础上结合了$L_1$和$L_2$正则项：
$L_{EN} = \frac{1}{2n} L(\beta) + R_{EN} = \frac{1}{2n} \sum_{i=1}^n ||y_i - \hat{y}_i||_2^2 + \lambda p \sum_{j=1}^{k} ||\beta_j||_1 + \frac{\lambda (1-p)}{2} \sum_{j=1}^{k} ||\beta_j||_2^2 \tag{7.1.32}$
可以看到我们通过参数$p$控制$R_{L_1}$和$R_{L_2}$的组合，当$p=1$时，$R_{EN} \sim R_{L_1}$，当$p=0$时，$R_{EN} \sim \frac{1}{2} R_{L_2}$

在ElasticNe回归中，我们有2个超参数需要设置，如何选取一组合理的参数成为一个难题，这时，我们可以利用交叉验证的技术来判断参数的优劣，从而选取一组较为合适的调节因子$\lambda$。

python实现

现在，我们使用sklearn中的datasets.make_regression方法生成回归模型数据，并分别使用sklearn中的linear_model.Ridge类，linear_model.Lasso类和linear_model.ElasticNetCV类(ElasticNetCV可以通过交叉验证来设置参数$\lambda$和$p$)，建立以上三种回归模型。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNetCV
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression

#print(type(data_ron))
X, y = make_regression(n_samples = 300, n_features = 20, 
                       n_informative = 10, noise = 100, random_state = 1)

X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size = 0.2, random_state = 10)

#普通线性
model_lr = LinearRegression()
model_lr.fit(X_train, y_train)

#岭回归
model_r = Ridge(alpha = 0.5)
model_r.fit(X_train, y_train)

#Lasso回归
model_l = Lasso(alpha = 0.1)
model_l.fit(X_train, y_train)

#ElasticNet回归
model_en = ElasticNetCV(cv = 8)
model_en.fit(X_train, y_train)

print("lambda:{:.4f}, p:{}".format(model_en.alpha_, model_en.l1_ratio_))
#lambda:0.1727, p:0.5

可以看到，通过交叉验证得到的$\lambda$为0.1727，而$p$为默认值0.5。现在，我们查看这4个回归的得分。

def get_score(model):
    print("=====================================")
    print("train_score:{:.4f}, test_score:{:.4f}".format(model.score(X_train, y_train), 
          model.score(X_test, y_test)))


for item in [model_lr, model_r, model_l, model_en]:
    get_score(item)

=====================================
train_score:0.7862, test_score:0.6850
=====================================
train_score:0.7862, test_score:0.6858
=====================================
train_score:0.7862, test_score:0.6857
=====================================
train_score:0.7793, test_score:0.7040

由结果可知，虽然在训练集上前3个模型得分一致，但是在测试集上岭回归和Lasso回归的得分略高于普通线性回归，且虽然ElasticNet回归在训练集上的得分小于前三个回归模型，但是它在测试集上的得分是最高的。可以看到，我们通过给普通线性模型添加正则项的方式，提高了模型泛化能力的作用。

参考文献

[1]唐亘. 精通数据科学 从线性回归到深度学习[M]. 人民邮电出版社, 2018.

[2]诸葛越. 百面机器学习[M]. 人民邮电出版社, 2018.

[3]雷明. 机器学习与应用[M]. 清华大学出版社, 2018.

[4]李航. 统计学习方法[M]. 清华大学出版社, 2019.

[5] Aurélien Géron. Hands-On Machine Learning with Scikit-Learn and TensorFlow[M]. O’Reilly Media, 2017.

[6] James, G, Witten, D, Hastie, T, Tibshirani, R. An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics[M]. Springer, 2013.

[7] Peter Harrington. Machine Learning in Action[M].Manning Publications, 2012.

[8]https://scikit-learn.org/stable/user_guide.html