线性模型

概述

线性模型是自然界最简单的模型之一,它描述了一个(或多个)自变量对另一个因变量的影响是呈简单的比例、线性关系.例如:

  • 住房每平米单价为1万元,100平米住房价格为100万元,120平米住房为120万元;
  • 一台挖掘机每小时挖100m3100m^3沙土,工作4小时可以挖掘400m3400m^3沙土.

线性模型在二维空间内表现为一条直线,在三维空间内表现为一个平面,更高维度下的线性模型很难用几何图形来表示(称为超平面).如下图所示:

二维空间下线性模型表现为一条直线

三维空间下线性模型表现为一个平面
线性回归是要根据一组输入值和输出值(称为样本),寻找一个线性模型,能最佳程度上拟合于给定的数值分布,从而再给定新的输入时预测输出.样本如下表所示:
输入(x) 输出(y)
0.5 5.0
0.6 5.5
0.8 6.0
1.1 6.8
1.4 6.8

根据样本拟合的线性模型如下图所示:

linear_3

线性模型定义

设给定一组属性x,x=(x1;x2;...;xn)x, x=(x_1;x_2;...;x_n),线性方程的一般表达形式为:

y=w1x1+w2x2+w3x3+...+wnxn+by = w_1x_1 + w_2x_2 + w_3x_3 + ... + w_nx_n + b

写成向量形式为:

y=wTx+by = w^Tx + b

其中,w=(w1;w2;...;wn),x=(x1;x2;...;xn)w=(w_1;w_2;...;w_n), x=(x_1;x_2;...;x_n),w和b经过学习后,模型就可以确定. 当自变量数量为1时,上述线性模型即为平面下的直线方程:

y=wx+by = wx + b

线性模型形式简单、易于建模,却蕴含着机器学习中一些重要的基本思想. 许多功能强大的非线性模型可以在线性模型基础上引入层级结构或高维映射而得. 此外,由于ww直观表达了各属性在预测中的重要性,因此线性模型具有很好的可解释性.例如,判断一个西瓜是否为好瓜,可以用如下表达式来判断:

f好瓜(x)=0.2x色泽+0.5x根蒂+0.3x敲声+1f_{\text{好瓜}}(x) = 0.2x_{\text{色泽}} + 0.5x_{\text{根蒂}} + 0.3x_{\text{敲声}} + 1

上述公式可以解释为,一个西瓜是否为好瓜,可以通过色泽、根蒂、敲声等因素共同判断,其中根蒂最重要(权重最高),其次是敲声和色泽.

模型训练

在二维平面中,给定两点可以确定一条直线.但在实际工程中,可能有很多个样本点,无法找到一条直线精确穿过所有样本点,只能找到一条与样本”足够接近“或”距离足够小“的直线,近似拟合给定的样本.如下图所示:

linear_4

如何确定直线到所有样本足够近呢?可以使用损失函数来进行度量.

损失函数

损失函数用来度量真实值(由样本中给出)和预测值(由模型算出)之间的差异.损失函数值越小,表明模型预测值和真实值之间差异越小,模型性能越好;损失函数值越大,模型预测值和真实值之间差异越大,模型性能越差.在回归问题中,均方差是常用的损失函数,其表达式如下所示:

E=12i=1n(yy)2E = \frac{1}{2}\sum_{i=1}^{n}{(y - y')^2}

其中,y为模型预测值,y’为真实值. 均方差具有非常好的几何意义,对应着常用的欧几里得距离(简称欧式距离). 线性回归的任务是要寻找最优线性模型,是的损失函数值最小,即:

(w,b)=argmin12i=1n(yy)2=argmin12i=1n(ywxib)2(w^*, b^*) = arg min \frac{1}{2}\sum_{i=1}^{n}{(y - y')^2} \\ = arg min \frac{1}{2}\sum_{i=1}^{n}{(y' - wx_i - b)^2}

基于均方误差最小化来进行模型求解的方法称为“最小二乘法”. 线性回归中,最小二乘法就是试图找到一条直线,是所有样本到直线的欧式距离之和最小. 可以将损失函数对w和b分别求导,得到损失函数的导函数,并令导函数为0即可得到w和b的最优解.

梯度下降法

① 为什么使用梯度下降

在实际计算中,通过最小二乘法求解最优参数有一定的问题:

(1)最小二乘法需要计算逆矩阵,有可能逆矩阵不存在;

(2)当样本特征数量较多时,计算逆矩阵非常耗时甚至不可行.

所以,在实际计算中,通常采用梯度下降法来求解损失函数的极小值,从而找到模型的最优参数.

② 什么是梯度下降

梯度(gradient)是一个向量(矢量,有方向),表示某一函数在该点处的方向导数沿着该方向取得最大值,即函数在该点处沿着该方向(此梯度的方向)变化最快,变化率最大.损失函数沿梯度相反方向收敛最快(即能最快找到极值点).当梯度向量为零(或接近于零),说明到达一个极值点,这也是梯度下降算法迭代计算的终止条件.

这种按照负梯度不停地调整函数权值的过程就叫作“梯度下降法”.通过这样的方法,改变权重让损失函数的值下降得更快,进而将值收敛到损失函数的某个极小值.

通过损失函数,我们将“寻找最优参数”问题,转换为了“寻找损失函数最小值”问题.梯度下降法算法描述如下:

(1)损失是否足够小?如果不是,计算损失函数的梯度.
(2)按梯度的反方向走一小步,以缩小损失.
(3)循环到(1).

gradent_decent

梯度下降法中通过沿着梯度负方向不断调整参数,从而逐步接近损失函数极小值所在点. 如下图所示:

gradent_decent2

③ 参数更新法则

在直线方程中,有两个参数需要学习,w0w_0w1w_1,梯度下降过程中,分别对这两个参数单独进行调整,调整法则如下:

w0=w0+Δw0w1=w1+Δw1w_0 = w_0 + \Delta w_0\\ w_1 = w_1 + \Delta w_1

Δw0\Delta w_0Δw1\Delta w_1可表示为:

Δw0=ηΔlossΔw0Δw1=ηΔlossΔw1\Delta w_0 = -\eta \frac{\Delta loss}{\Delta w_0}\\ \Delta w_1 = -\eta \frac{\Delta loss}{\Delta w_1}\\

其中,η\eta称为学习率,ΔlossΔwi\frac{\Delta loss}{\Delta w_i}为梯度(即损失函数关于参数wiw_i的偏导数). 损失函数表达式为:

loss=12(yy)2=12((y(w0+w1x))2)loss =\frac{1}{2}\sum(y - y')^2 = \frac{1}{2}\sum((y-(w_0+w_1x))^2)

对损失函数求导(求导过程见补充知识),可得w0,w1w_0, w_1的偏导数为:

ΔlossΔw0=((yy)(1))=(yy)ΔlossΔw1=((yy)(x))=(x(yy))\frac{\Delta loss}{\Delta w_0} = \sum((y - y')(-1)) = -\sum(y - y')\\ \frac{\Delta loss}{\Delta w_1} = \sum((y - y')(-x)) = -\sum(x(y - y'))

实现线性回归

① 自己编码实现

以下是实现线性回归的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
# 线性回归示例
import numpy as np
import matplotlib.pyplot as mp
from mpl_toolkits.mplot3d import axes3d
import sklearn.preprocessing as sp

# 训练数据集
train_x = np.array([0.5, 0.6, 0.8, 1.1, 1.4]) # 输入集
train_y = np.array([5.0, 5.5, 6.0, 6.8, 7.0]) # 输出集

n_epochs = 1000 # 迭代次数
lrate = 0.01 # 学习率
epochs = [] # 记录迭代次数
losses = [] # 记录损失值

w0, w1 = [1], [1] # 模型初始值

for i in range(1, n_epochs + 1):
epochs.append(i) # 记录第几次迭代

y = w0[-1] + w1[-1] * train_x # 取出最新的w0,w1计算线性方程输出
# 损失函数(均方差)
loss = (((train_y - y) ** 2).sum()) / 2
losses.append(loss) # 记录每次迭代的损失值

print("%d: w0=%f, w1=%f, loss=%f" % (i, w0[-1], w1[-1], loss))

# 计算w0,w1的偏导数
d0 = -(train_y - y).sum()
d1 = -(train_x * (train_y - y)).sum()

# 更新w0,w1
w0.append(w0[-1] - (d0 * lrate))
w1.append(w1[-1] - (d1 * lrate))

程序执行结果:

1
2
3
4
5
6
7
8
9
10
11
1 w0=1.00000000 w1=1.00000000 loss=44.17500000
2 w0=1.20900000 w1=1.19060000 loss=36.53882794
3 w0=1.39916360 w1=1.36357948 loss=30.23168666
4 w0=1.57220792 w1=1.52054607 loss=25.02222743
5 w0=1.72969350 w1=1.66296078 loss=20.71937337
......
996 w0=4.06506160 w1=2.26409126 loss=0.08743506
997 w0=4.06518850 w1=2.26395572 loss=0.08743162
998 w0=4.06531502 w1=2.26382058 loss=0.08742820
999 w0=4.06544117 w1=2.26368585 loss=0.08742480
1000 w0=4.06556693 w1=2.26355153 loss=0.08742142

可以给数据加上可视化,让结果更直观.添加如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
###################### 训练过程可视化 ######################
# 训练过程可视化
## 损失函数收敛过程
w0 = np.array(w0[:-1])
w1 = np.array(w1[:-1])

mp.figure("Losses", facecolor="lightgray") # 创建一个窗体
mp.title("epoch", fontsize=20)
mp.ylabel("loss", fontsize=14)
mp.grid(linestyle=":") # 网格线:虚线
mp.plot(epochs, losses, c="blue", label="loss")
mp.legend() # 图例
mp.tight_layout() # 紧凑格式

## 显示模型直线
pred_y = w0[-1] + w1[-1] * train_x # 根据x预测y
mp.figure("Linear Regression", facecolor="lightgray")
mp.title("Linear Regression", fontsize=20)
mp.xlabel("x", fontsize=14)
mp.ylabel("y", fontsize=14)
mp.grid(linestyle=":")
mp.scatter(train_x, train_y, c="blue", label="Traing") # 绘制样本散点图
mp.plot(train_x, pred_y, c="red", label="Regression")
mp.legend()

# 显示梯度下降过程(复制粘贴即可,不需要编写)
# 计算损失函数曲面上的点 loss = f(w0, w1)
arr1 = np.linspace(0, 10, 500) # 0~9间产生500个元素的均匀列表
arr2 = np.linspace(0, 3.5, 500) # 0~3.5间产生500个元素的均匀列表

grid_w0, grid_w1 = np.meshgrid(arr1, arr2) # 产生二维矩阵

flat_w0, flat_w1 = grid_w0.ravel(), grid_w1.ravel() # 二维矩阵扁平化
loss_metrix = train_y.reshape(-1, 1) # 生成误差矩阵(-1,1)表示自动计算维度
outer = np.outer(train_x, flat_w1) # 求外积(train_x和flat_w1元素两两相乘的新矩阵)
# 计算损失:((w0 + w1*x - y)**2)/2
flat_loss = (((flat_w0 + outer - loss_metrix) ** 2).sum(axis=0)) / 2
grid_loss = flat_loss.reshape(grid_w0.shape)

mp.figure('Loss Function')
ax = mp.gca(projection='3d')
mp.title('Loss Function', fontsize=14)
ax.set_xlabel('w0', fontsize=14)
ax.set_ylabel('w1', fontsize=14)
ax.set_zlabel('loss', fontsize=14)
ax.plot_surface(grid_w0, grid_w1, grid_loss, rstride=10, cstride=10, cmap='jet')
ax.plot(w0, w1, losses, 'o-', c='orangered', label='BGD', zorder=5)
mp.legend(loc='lower left')

mp.show()

数据可视化结果如下图所示:

回归得到的线性模型

损失函数收敛过程

梯度下降过程

② 通过sklearn API实现

同样,可以使用sklearn库提供的API实现线性回归.代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# 利用LinearRegression实现线性回归
import numpy as np
import sklearn.linear_model as lm # 线性模型# 线性模型
import sklearn.metrics as sm # 模型性能评价模块
import matplotlib.pyplot as mp

train_x = np.array([[0.5], [0.6], [0.8], [1.1], [1.4]]) # 输入集
train_y = np.array([5.0, 5.5, 6.0, 6.8, 7.0]) # 输出集

# 创建线性回归器
model = lm.LinearRegression()
# 用已知输入、输出数据集训练回归器
model.fit(train_x, train_y)
# 根据训练模型预测输出
pred_y = model.predict(train_x)

print("coef_:", model.coef_) # 系数
print("intercept_:", model.intercept_) # 截距

# 可视化回归曲线
mp.figure('Linear Regression', facecolor='lightgray')
mp.title('Linear Regression', fontsize=20)
mp.xlabel('x', fontsize=14)
mp.ylabel('y', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')

# 绘制样本点
mp.scatter(train_x, train_y, c='blue', alpha=0.8, s=60, label='Sample')

# 绘制拟合直线
mp.plot(train_x, # x坐标数据
pred_y, # y坐标数据
c='orangered', label='Regression')

mp.legend()
mp.show()

执行结果:

linear_10

模型评价指标

(1)平均绝对误差(Mean Absolute Deviation):单个观测值与算术平均值的偏差的绝对值的平均;

(2)均方误差:单个样本到平均值差值的平方平均值;

(3)MAD(中位数绝对偏差):与数据中值绝对偏差的中值;

(4)R2决定系数:趋向于1,模型越好;趋向于0,模型越差.

多项式回归

什么是多项式回归

线性回归适用于数据呈线性分布的回归问题.如果数据样本呈明显非线性分布,线性回归模型就不再适用(下图左),而采用多项式回归可能更好(下图右).例如:

poly_1

多项式模型定义

与线性模型相比,多项式模型引入了高次项,自变量的指数大于1,例如一元二次方程:

y=w0+w1x+w2x2y = w_0 + w_1x + w_2x^2

一元三次方程:

y=w0+w1x+w2x2+w3x3y = w_0 + w_1x + w_2x^2 + w_3x ^ 3

推广到一元n次方程:

y=w0+w1x+w2x2+w3x3+...+wnxny = w_0 + w_1x + w_2x^2 + w_3x ^ 3 + ... + w_nx^n

上述表达式可以简化为:

y=i=1Nwixiy = \sum_{i=1}^N w_ix^i

与线性回归的关系

多项式回归可以理解为线性回归的扩展,在线性回归模型中添加了新的特征值.例如,要预测一栋房屋的价格,有x1,x2,x3x_1, x_2, x_3三个特征值,分别表示房子长、宽、高,则房屋价格可表示为以下线性模型:

y=w1x1+w2x2+w3x3+by = w_1 x_1 + w_2 x_2 + w_3 x_3 + b

对于房屋价格,也可以用房屋的体积,而不直接使用x1,x2,x3x_1, x_2, x_3三个特征:

y=w0+w1x+w2x2+w3x3y = w_0 + w_1x + w_2x^2 + w_3x ^ 3

相当于创造了新的特征x,xx, x = 长 * 宽 * 高. 以上两个模型可以解释为:

  • 房屋价格是关于长、宽、高三个特征的线性模型
  • 房屋价格是关于体积的多项式模型

因此,可以将一元n次多项式变换成n元一次线性模型.

多项式回归实现

对于一元n次多项式,同样可以利用梯度下降对损失值最小化的方法,寻找最优的模型参数w0,w1,w2,...,wnw_0, w_1, w_2, ..., w_n.可以将一元n次多项式,变换成n元一次多项式,求线性回归.以下是一个多项式回归的实现.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# 多项式回归示例
import numpy as np
# 线性模型
import sklearn.linear_model as lm
# 模型性能评价模块
import sklearn.metrics as sm
import matplotlib.pyplot as mp
# 管线模块
import sklearn.pipeline as pl
import sklearn.preprocessing as sp

train_x, train_y = [], [] # 输入、输出样本
with open("poly_sample.txt", "rt") as f:
for line in f.readlines():
data = [float(substr) for substr in line.split(",")]
train_x.append(data[:-1])
train_y.append(data[-1])

train_x = np.array(train_x) # 二维数据形式的输入矩阵,一行一样本,一列一特征
train_y = np.array(train_y) # 一维数组形式的输出序列,每个元素对应一个输入样本
# print(train_x)
# print(train_y)

# 将多项式特征扩展预处理,和一个线性回归器串联为一个管线
# 多项式特征扩展:对现有数据进行的一种转换,通过将数据映射到更高维度的空间中
# 进行多项式扩展后,我们就可以认为,模型由以前的直线变成了曲线
# 从而可以更灵活的去拟合数据
# pipeline连接两个模型
model = pl.make_pipeline(sp.PolynomialFeatures(3), # 多项式特征扩展,扩展最高次项为3
lm.LinearRegression())

# 用已知输入、输出数据集训练回归器
model.fit(train_x, train_y)
# print(model[1].coef_)
# print(model[1].intercept_)

# 根据训练模型预测输出
pred_train_y = model.predict(train_x)

# 评估指标
err4 = sm.r2_score(train_y, pred_train_y) # R2得分, 范围[0, 1], 分值越大越好
print(err4)

# 在训练集之外构建测试集
test_x = np.linspace(train_x.min(), train_x.max(), 1000)
pre_test_y = model.predict(test_x.reshape(-1, 1)) # 对新样本进行预测

# 可视化回归曲线
mp.figure('Polynomial Regression', facecolor='lightgray')
mp.title('Polynomial Regression', fontsize=20)
mp.xlabel('x', fontsize=14)
mp.ylabel('y', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.scatter(train_x, train_y, c='dodgerblue', alpha=0.8, s=60, label='Sample')

mp.plot(test_x, pre_test_y, c='orangered', label='Regression')

mp.legend()
mp.show()

打印输出:

1
0.9224401504764776

执行结果:

poly_2

过拟合与欠拟合

什么是欠拟合、过拟合

在上一小节多项式回归示例中,多项特征扩展器PolynomialFeatures()进行多项式扩展时,指定了最高次数为3,该参数为多项式扩展的重要参数,如果选取不当,则可能导致不同的拟合效果.下图显示了该参数分别设为1、20时模型的拟合图像:

poly_3

这两种其实都不是好的模型. 前者没有学习到数据分布规律,模型拟合程度不够,预测准确度过低,这种现象称为“欠拟合”;后者过于拟合更多样本,以致模型泛化能力(新样本的适应性)变差,这种现象称为“过拟合”. **欠拟合模型一般表现为训练集、测试集下准确度都比较低;过拟合模型一般表现为训练集下准确度较高、测试集下准确度较低. **一个好的模型,不论是对于训练数据还是测试数据,都有接近的预测精度,而且精度不能太低.

【思考1】以下哪种模型较好,哪种模型较差,较差的原因是什么?

训练集R2值 测试集R2值
0.6 0.5
0.9 0.6
0.9 0.88

【答案】第一个模型欠拟合;第二个模型过拟合;第三个模型适中,为可接受的模型.

【思考2】以下哪个曲线为欠拟合、过拟合,哪个模型拟合最好?

overfit

【答案】第一个模型欠拟合;第三个模型过拟合;第二个模型拟合较好.

如何处理欠拟合、过拟合

  • 欠拟合:提高模型复杂度,如增加特征、增加模型最高次幂等等;
  • 过拟合:降低模型复杂度,如减少特征、降低模型最高次幂等等.

线性回归模型变种

正则化

什么是正则化

过拟合还有一个常见的原因,就是模型参数值太大,所以可以通过抑制参数的方式来解决过拟合问题.如下图所示,右图产生了一定程度过拟合,可以通过弱化高次项的系数(但不删除)来降低过拟合.

overfit2

例如,可以通过在θ3,θ4\theta_3, \theta_4上添加一定的系数,来压制这两个高次项的系数,这种方法称为正则化。但在实际问题中,可能有更多的系数,我们并不知道应该压制哪些系数,所以,可以通过收缩所有系数来避免过拟合.

正则化的定义

正则化是指,在目标函数(如损失函数)后面加上一个范数,来防止过拟合的手段,这个范数定义为:

xp=(i=1Nxp)1p||x||_p = (\sum_{i=1}^N |x|^p)^{\frac{1}{p}}

当p=1时,称为L1范数(即所有系数绝对值之和):

x1=(i=1Nx)||x||_1 = (\sum_{i=1}^N |x|)

当p=2是,称为L2范数(即所有系数平方之和再开方):

x2=(i=1Nx2)12||x||_2 = (\sum_{i=1}^N |x|^2)^{\frac{1}{2}}

通过对目标函数添加正则项,整体上压缩了参数的大小,从而防止过拟合.

Lasso回归与岭回归

Lasso 回归和岭回归(Ridge Regression)都是在标准线性回归的基础上修改了损失函数的回归算法. Lasso回归全称为 Least absolute shrinkage and selection operator,又译“最小绝对值收敛和选择算子”、”套索算法”,其损失函数如下所示:

E=1n(i=1Nyiyi)2+λw1E = \frac{1}{n}(\sum_{i=1}^N y_i - y_i')^2 + \lambda ||w||_1

岭回归损失函数为:

E=1n(i=1Nyiyi)2+λw2E = \frac{1}{n}(\sum_{i=1}^N y_i - y_i')^2 + \lambda ||w||_2

从逻辑上说,Lasso回归和岭回归都可以理解为通过调整损失函数,减小函数的系数,从而避免过于拟合于样本,降低偏差较大的样本的权重和对模型的影响程度.

以下关于Lasso回归于岭回归的sklearn实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# Lasso回归和岭回归示例
import numpy as np
# 线性模型
import sklearn.linear_model as lm
# 模型性能评价模块
import sklearn.metrics as sm
import matplotlib.pyplot as mp

x, y = [], [] # 输入、输出样本
with open("abnormal.txt", "rt") as f:
for line in f.readlines():
data = [float(substr) for substr in line.split(",")]
x.append(data[:-1])
y.append(data[-1])

x = np.array(x) # 二维数据形式的输入矩阵,一行一样本,一列一特征
y = np.array(y) # 一维数组形式的输出序列,每个元素对应一个输入样本
# print(x)
# print(y)

# 创建线性回归器
model = lm.LinearRegression()
# 用已知输入、输出数据集训练回归器
model.fit(x, y)
# 根据训练模型预测输出
pred_y = model.predict(x)

# 创建岭回归器并进行训练
# Ridge: 第一个参数为正则强度,该值越大,异常样本权重就越小
model_2 = lm.Ridge(alpha=200, max_iter=1000) # 创建对象, max_iter为最大迭代次数
model_2.fit(x, y) # 训练
pred_y2 = model_2.predict(x) # 预测

# lasso回归
model_3 = lm.Lasso(alpha=0.5, # L1范数相乘的系数
max_iter=1000) # 最大迭代次数
model_3.fit(x, y) # 训练
pred_y3 = model_3.predict(x) # 预测

# 可视化回归曲线
mp.figure('Linear & Ridge & Lasso', facecolor='lightgray')
mp.title('Linear & Ridge & Lasso', fontsize=20)
mp.xlabel('x', fontsize=14)
mp.ylabel('y', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.scatter(x, y, c='dodgerblue', alpha=0.8, s=60, label='Sample')
sorted_idx = x.T[0].argsort()

mp.plot(x[sorted_idx], pred_y[sorted_idx], c='orangered', label='Linear') # 线性回归
mp.plot(x[sorted_idx], pred_y2[sorted_idx], c='limegreen', label='Ridge') # 岭回归
mp.plot(x[sorted_idx], pred_y3[sorted_idx], c='blue', label='Lasso') # Lasso回归

mp.legend()
mp.show()

以下是执行结果:

linear_9

模型保存与加载

可以使用Python提供的功能对模型对象进行保存.使用方法如下:

1
2
3
4
5
import pickle
# 保存模型
pickle.dump(模型对象, 文件对象)
# 加载模型
model_obj = pickle.load(文件对象)

保存训练模型应该在训练完成或评估完成之后,完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 模型保存示例
import numpy as np
import sklearn.linear_model as lm # 线性模型
import pickle

x = np.array([[0.5], [0.6], [0.8], [1.1], [1.4]]) # 输入集
y = np.array([5.0, 5.5, 6.0, 6.8, 7.0]) # 输出集

# 创建线性回归器
model = lm.LinearRegression()
# 用已知输入、输出数据集训练回归器
model.fit(x, y)

print("训练完成.")

# 保存训练后的模型
with open('linear_model.pkl', 'wb') as f:
pickle.dump(model, f)
print("保存模型完成.")

执行完成后,可以看到与源码相同目录下多了一个名称为linear_model.pkl的文件,这就是保存的训练模型.使用该模型代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 模型加载示例
import numpy as np
import sklearn.linear_model as lm # 线性模型
import sklearn.metrics as sm # 模型性能评价模块
import matplotlib.pyplot as mp
import pickle

x = np.array([[0.5], [0.6], [0.8], [1.1], [1.4]]) # 输入集
y = np.array([5.0, 5.5, 6.0, 6.8, 7.0]) # 输出集

# 加载模型
with open('linear_model.pkl', 'rb') as f:
model = pickle.load(f)
print("加载模型完成.")

# 根据加载的模型预测输出
pred_y = model.predict(x)

# 可视化回归曲线
mp.figure('Linear Regression', facecolor='lightgray')
mp.title('Linear Regression', fontsize=20)
mp.xlabel('x', fontsize=14)
mp.ylabel('y', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.scatter(x, y, c='blue', alpha=0.8, s=60, label='Sample')

mp.plot(x, pred_y, c='orangered', label='Regression')

mp.legend()
mp.show()

执行结果和训练模型预测结果一样.

线性回归总结

(1)什么是线性模型:线性模型是自然界最简单的模型之一,反映自变量、因变量之间的等比例增长关系

(2)什么时候使用线性回归:线性模型只能用于满足线性分布规律的数据中

(3)如何实现线性回归:给定一组样本,给定初始的w和b,通过梯度下降法求最优的w和b

决策树回归

决策树

基本算法原理

核心思想:相似的输入必会产生相似的输出。例如预测某人薪资:

年龄:1-青年,2-中年,3-老年
学历:1-本科,2-硕士,3-博士
经历:1-出道,2-一般,3-老手,4-骨灰
性别:1-男性,2-女性

年龄 学历 经历 性别 ==> 薪资
1 1 1 1 ==> 6000(低)
2 1 3 1 ==> 10000(中)
3 3 4 1 ==> 50000(高)
==>
1 3 2 2 ==> ?
1
2
3
4
5
6
7
8
样本数量非常庞大  100W个样本

换一种数据结构,来提高检索效率

树形结构

回归 : 均值
分类 : 投票(概率)

为了提高搜索效率,使用树形数据结构处理样本数据:

年龄=1{学历1学历2学历3年龄=2{学历1学历2学历3年龄=3{学历1学历2学历3\text{年龄}=1\left\{ \begin{aligned} \text{学历}1 \\ \text{学历}2 \\ \text{学历}3 \\ \end{aligned} \right. \quad\quad \text{年龄}=2\left\{ \begin{aligned} \text{学历}1 \\ \text{学历}2 \\ \text{学历}3 \\ \end{aligned} \right. \quad\quad \text{年龄}=3\left\{ \begin{aligned} \text{学历}1 \\ \text{学历}2 \\ \text{学历}3 \\ \end{aligned} \right.

首先从训练样本矩阵中选择一个特征进行子表划分,使每个子表中该特征的值全部相同,然后再在每个子表中选择下一个特征按照同样的规则继续划分更小的子表,不断重复直到所有的特征全部使用完为止,此时便得到叶级子表,其中所有样本的特征值全部相同。对于待预测样本,根据其每一个特征的值,选择对应的子表,逐一匹配,直到找到与之完全匹配的叶级子表,用该子表中样本的输出,通过平均(回归)或者投票(分类)为待预测样本提供输出。

首先选择哪一个特征进行子表划分决定了决策树的性能。这么多特征,使用哪个特征先进行子表划分?

sklearn提供的决策树底层为cart树(Classification and Regression Tree),cart回归树在解决回归问题时的步骤如下:

  1. 原始数据集S,此时树的深度depth=0;

  2. 针对集合S,遍历每一个特征的每一个value(遍历数据中的所有离散值 (12个) )

    用该value将原数据集S分裂成2个集合:左集合left(<=value的样本)、右集合right(>value的样本),

    分别计算这2个集合的mse(均方误差),找到使(left_mse+right_mse)最小的那个value,记录下此时的特征名称和value,这个就是最佳分割特征以及最佳分割值;

    ​ mse:均方误差

    ​ ((y1-y1’)^2 + (y2-y2’)^2 + (y3-y3’)^2 + (y4-y4’)^2 ) / 4 = mse

x1 y1 y1’

x2 y2 y2’

x3 y3 y3’

x4 y4 y4’

  1. 找到最佳分割特征以及最佳分割value之后,用该value将集合S分裂成2个集合,depth+=1;

  2. 针对集合left、right分别重复步骤2,3,直到达到终止条件。

决策树底层结构 为二叉树

1
2
3
4
5
终止条件有如下几种:
1、特征已经用完了:没有可供使用的特征再进行分裂了,则树停止分裂;
2、子节点中没有样本了:此时该结点已经没有样本可供划分,该结点停止分裂;
3、树达到了人为预先设定的最大深度:depth >= max_depth,树停止分裂。
4、节点的样本数量达到了人为设定的阈值:样本数量 < min_samples_split ,则该节点停止分裂;

决策树回归器模型相关API:

1
2
3
4
5
6
7
8
9
10
import sklearn.tree as st

# 创建决策树回归器模型 决策树的最大深度为4
model = st.DecisionTreeRegressor(max_depth=4)
# 训练模型
# train_x: 二维数组样本数据
# train_y: 训练集中对应每行样本的结果
model.fit(train_x, train_y)
# 测试模型
pred_test_y = model.predict(test_x)

案例:预测波士顿地区房屋价格。

  1. 读取数据,打断原始数据集。 划分训练集和测试集。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
import sklearn.datasets as sd
import sklearn.utils as su
# 加载波士顿地区房价数据集
boston = sd.load_boston()
print(boston.feature_names)
# |CRIM|ZN|INDUS|CHAS|NOX|RM|AGE|DIS|RAD|TAX|PTRATIO|B|LSTAT|
# 犯罪率|住宅用地比例|商业用地比例|是否靠河|空气质量|房间数|年限|距中心区距离|路网密度|房产税|师生比|黑人比例|低地位人口比例|
# 打乱原始数据集的输入和输出
x, y = su.shuffle(boston.data, boston.target, random_state=7)
# 划分训练集和测试集
train_size = int(len(x) * 0.8)
train_x, test_x, train_y, test_y = \
x[:train_size], x[train_size:], \
y[:train_size], y[train_size:]
  1. 创建决策树回归器模型,使用训练集训练模型。使用测试集测试模型。
1
2
3
4
5
6
7
8
9
10
import sklearn.tree as st
import sklearn.metrics as sm

# 创建决策树回归模型
model = st.DecisionTreeRegressor(max_depth=4)
# 训练模型
model.fit(train_x, train_y)
# 测试模型
pred_test_y = model.predict(test_x)
print(sm.r2_score(test_y, pred_test_y))

集合算法

​ 三个臭皮匠,顶个诸葛亮

单个模型得到的预测结果总是片面的,根据多个不同模型给出的预测结果,利用平均(回归)或者投票(分类)的方法,得出最终预测结果。

基于决策树的集合算法,就是按照某种规则,构建多棵彼此不同的决策树模型,分别给出针对未知样本的预测结果,最后通过平均或投票得到相对综合的结论。常用的集合模型包括Boosting类模型(AdaBoost、GBDT)与Bagging(自助聚合、随机森林)类模型。

AdaBoost模型(正向激励)

首先为样本矩阵中的样本随机分配初始权重,由此构建一棵带有权重的决策树,在由该决策树提供预测输出时,通过加权平均或者加权投票的方式产生预测值。

1
2
已经构建好一个决策树  通过1322 找到所有的女博士  一个4个  6000  8000 9000 10000
由于正向激励,对每个样本都分配了初始权重 权重为:1 1 1 3 预测: 加权均值

将训练样本代入模型,预测其输出,对那些预测值与实际值不同的样本,提高其权重,由此形成第二棵决策树。重复以上过程,构建出不同权重的若干棵决策树。

1
实际值:10000  但是你预测的为6000   构建第二个决策树,提高10000样本的权重

正向激励相关API:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import sklearn.tree as st
import sklearn.ensemble as se
# model: 决策树模型(一颗)
model = st.DecisionTreeRegressor(max_depth=4)
# 自适应增强决策树回归模型
#
model = se.AdaBoostRegressor(model, n_estimators=400, random_state=7)

正向激励 的基础模型 : 决策树
n_estimators:构建400棵不同权重的决策树,训练模型

# 训练模型
model.fit(train_x, train_y)
# 测试模型
pred_test_y = model.predict(test_x)

案例:基于正向激励训练预测波士顿地区房屋价格的模型。

1
2
3
4
5
6
7
8
# 创建基于决策树的正向激励回归器模型
model = se.AdaBoostRegressor(
st.DecisionTreeRegressor(max_depth=4), n_estimators=400, random_state=7)
# 训练模型
model.fit(train_x, train_y)
# 测试模型
pred_test_y = model.predict(test_x)
print(sm.r2_score(test_y, pred_test_y))

特征重要性

作为决策树模型训练过程的副产品,根据划分子表时选择特征的顺序标志了该特征的重要程度,此即为该特征重要性指标。训练得到的模型对象提供了属性:feature_importances_来存储每个特征的重要性。

获取样本矩阵特征重要性属性:

1
2
model.fit(train_x, train_y)
fi = model.feature_importances_

案例:获取普通决策树与正向激励决策树训练的两个模型的特征重要性值,按照从大到小顺序输出绘图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import matplotlib.pyplot as mp

model = st.DecisionTreeRegressor(max_depth=4)
model.fit(train_x, train_y)
# 决策树回归器给出的特征重要性
fi_dt = model.feature_importances_
model = se.AdaBoostRegressor(
st.DecisionTreeRegressor(max_depth=4), n_estimators=400, random_state=7)
model.fit(train_x, train_y)
# 基于决策树的正向激励回归器给出的特征重要性
fi_ab = model.feature_importances_

mp.figure('Feature Importance', facecolor='lightgray')
mp.subplot(211)
mp.title('Decision Tree', fontsize=16)
mp.ylabel('Importance', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(axis='y', linestyle=':')
sorted_indices = fi_dt.argsort()[::-1]
pos = np.arange(sorted_indices.size)
mp.bar(pos, fi_dt[sorted_indices], facecolor='deepskyblue', edgecolor='steelblue')
mp.xticks(pos, feature_names[sorted_indices], rotation=30)
mp.subplot(212)
mp.title('AdaBoost Decision Tree', fontsize=16)
mp.ylabel('Importance', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(axis='y', linestyle=':')
sorted_indices = fi_ab.argsort()[::-1]
pos = np.arange(sorted_indices.size)
mp.bar(pos, fi_ab[sorted_indices], facecolor='lightcoral', edgecolor='indianred')
mp.xticks(pos, feature_names[sorted_indices], rotation=30)
mp.tight_layout()
mp.show()

GBDT

GBDT(Gradient Boosting Decision Tree 梯度提升树)通过多轮迭代,每轮迭代产生一个弱分类器,每个分类器在上一轮分类器的残差**(残差在数理统计中是指实际观察值与估计值(拟合值)之间的差)**基础上进行训练。基于预测结果的残差设计损失函数。GBDT训练的过程即是求该损失函数最小值的过程。

案例

GBDT案例

GBDT原理

原理1:

GBDT原理1

原理2:

GBDT原理2

原理3:

GBDT原理3

1
2
3
4
5
6
7
8
9
10
import sklearn.tree as st
import sklearn.ensemble as se
# 自适应增强决策树回归模型
# n_estimators:构建400棵不同权重的决策树,训练模型
model = se.GridientBoostingRegressor(
max_depth=10, n_estimators=1000, min_samples_split=2)
# 训练模型
model.fit(train_x, train_y)
# 测试模型
pred_test_y = model.predict(test_x)
1
boosting : Adaboost  GBDT

自助聚合

每次从总样本矩阵中以有放回抽样的方式随机抽取部分样本构建决策树,这样形成多棵包含不同训练样本的决策树,以削弱某些强势样本对模型预测结果的影响,提高模型的泛化特性。

随机森林

在自助聚合的基础上,每次构建决策树模型时,不仅随机选择部分样本,而且还随机选择部分特征,这样的集合算法,不仅规避了强势样本对预测结果的影响,而且也削弱了强势特征的影响,使模型的预测能力更加泛化。

随机森林相关API:

1
2
3
4
5
6
7
import sklearn.ensemble as se
# 随机森林回归模型 (属于集合算法的一种)
# max_depth:决策树最大深度10
# n_estimators:构建1000棵决策树,训练模型
# min_samples_split: 子表中最小样本数 若小于这个数字,则不再继续向下拆分
model = se.RandomForestRegressor(
max_depth=10, n_estimators=1000, min_samples_split=2)

案例:分析共享单车的需求,从而判断如何进行共享单车的投放。

1
2
3
1。加载并整理数据集
2.特征分析
3.打乱数据集,划分训练集,测试集
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
import numpy as np
import sklearn.utils as su
import sklearn.ensemble as se
import sklearn.metrics as sm
import matplotlib.pyplot as mp

data = np.loadtxt('../data/bike_day.csv', unpack=False, dtype='U20', delimiter=',')
day_headers = data[0, 2:13]
x = np.array(data[1:, 2:13], dtype=float)
y = np.array(data[1:, -1], dtype=float)

x, y = su.shuffle(x, y, random_state=7)
print(x.shape, y.shape)
train_size = int(len(x) * 0.9)
train_x, test_x, train_y, test_y = \
x[:train_size], x[train_size:], y[:train_size], y[train_size:]
# 随机森林回归器
model = se.RandomForestRegressor( max_depth=10, n_estimators=1000, min_samples_split=2)
model.fit(train_x, train_y)
# 基于“天”数据集的特征重要性
fi_dy = model.feature_importances_
pred_test_y = model.predict(test_x)
print(sm.r2_score(test_y, pred_test_y))

data = np.loadtxt('../data/bike_hour.csv', unpack=False, dtype='U20', delimiter=',')
hour_headers = data[0, 2:13]
x = np.array(data[1:, 2:13], dtype=float)
y = np.array(data[1:, -1], dtype=float)
x, y = su.shuffle(x, y, random_state=7)
train_size = int(len(x) * 0.9)
train_x, test_x, train_y, test_y = \
x[:train_size], x[train_size:], \
y[:train_size], y[train_size:]
# 随机森林回归器
model = se.RandomForestRegressor(
max_depth=10, n_estimators=1000,
min_samples_split=2)
model.fit(train_x, train_y)
# 基于“小时”数据集的特征重要性
fi_hr = model.feature_importances_
pred_test_y = model.predict(test_x)
print(sm.r2_score(test_y, pred_test_y))

画图显示两组样本数据的特征重要性:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
mp.figure('Bike', facecolor='lightgray')
mp.subplot(211)
mp.title('Day', fontsize=16)
mp.ylabel('Importance', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(axis='y', linestyle=':')
sorted_indices = fi_dy.argsort()[::-1]
pos = np.arange(sorted_indices.size)
mp.bar(pos, fi_dy[sorted_indices], facecolor='deepskyblue', edgecolor='steelblue')
mp.xticks(pos, day_headers[sorted_indices], rotation=30)

mp.subplot(212)
mp.title('Hour', fontsize=16)
mp.ylabel('Importance', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(axis='y', linestyle=':')
sorted_indices = fi_hr.argsort()[::-1]
pos = np.arange(sorted_indices.size)
mp.bar(pos, fi_hr[sorted_indices], facecolor='lightcoral', edgecolor='indianred')
mp.xticks(pos, hour_headers[sorted_indices], rotation=30)
mp.tight_layout()
mp.show()

补充知识

R2系数详细计算

R2系数详细计算过程如下:

若用yiy_i表示真实的观测值,用yˉ\bar{y}表示真实观测值的平均值,用yi^\hat{y_i}表示预测值则,有以下评估指标:

  • 回归平方和(SSR)

SSR=i=1n(yi^yˉ)2SSR = \sum_{i=1}^{n}(\hat{y_i} - \bar{y})^2

估计值与平均值的误差,反映自变量与因变量之间的相关程度的偏差平方和.

  • 残差平方和(SSE)

SSE=i=1n(yiyi^)2SSE = \sum_{i=1}^{n}(y_i-\hat{y_i} )^2

即估计值与真实值的误差,反映模型拟合程度.

  • 总离差平方和(SST)

SST=SSR+SSE=i=1n(yiyˉ)2SST =SSR + SSE= \sum_{i=1}^{n}(y_i - \bar{y})^2

即平均值与真实值的误差,反映与数学期望的偏离程度.

  • R2_score计算公式

R2_score,即决定系数,反映因变量的全部变异能通过回归关系被自变量解释的比例.计算公式:

R2=1SSESSTR^2=1-\frac{SSE}{SST}

即:

R2=1i=1n(yiy^i)2i=1n(yiyˉ)2R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)2}{\sum_{i=1}{n} (y_i - \bar{y})^2}

进一步化简为:

R2=1i(yiyi)2/ni(yiy^)2/n=1RMSEVarR^2 = 1 - \frac{\sum\limits_i(y_i - y_i)^2 / n}{\sum\limits_i(y_i - \hat{y})^2 / n} = 1 - \frac{RMSE}{Var}

分子就变成了常用的评价指标均方误差MSE,分母就变成了方差,对于R2R^2可以通俗地理解为使用均值作为误差基准,看预测误差是否大于或者小于均值基准误差.

R2_score = 1,样本中预测值和真实值完全相等,没有任何误差,表示回归分析中自变量对因变量的解释越好.

R2_score = 0,此时分子等于分母,样本的每项预测值都等于均值.

线性回归损失函数求导过程

线性函数定义为:

y=w0+w0x1y = w_0 + w_0 x_1

采用均方差损失函数:

loss=12(yy)2loss = \frac{1}{2} (y - y')^2

其中,y为真实值,来自样本;y’为预测值,即线性方程表达式,带入损失函数得:

loss=12(y(w0+w1x1))2loss = \frac{1}{2} (y - (w_0 + w_1 x_1))^2

将该式子展开:

loss=12(y22y(w0+w1x1)+(w0+w1x1)2)12(y22yw02yw1x1+w02+2w0w1x1+w12x12)loss = \frac{1}{2} (y^2 - 2y(w_0 + w_1 x_1) + (w_0 + w_1 x_1)^2) \\ \frac{1}{2} (y^2 - 2y*w_0 - 2y*w_1x_1 + w_0^2 + 2w_0*w_1 x_1 + w_1^2x_1^2) \\

w0w_0求导:

lossw0=12(02y0+2w0+2w1x1+0)=12(2y+2w0+2w1x1)=122(y+(w0+w1x1))=(y+y)=(yy)\frac{\partial loss}{\partial w_0} = \frac{1}{2}(0-2y-0+2w_0 + 2w_1 x_1 +0) \\ =\frac{1}{2}(-2y + 2 w_0 + 2w_1 x_1) \\ = \frac{1}{2} * 2(-y + (w_0 + w_1 x_1)) \\ =(-y + y') = -(y - y')

w1w_1求导:

lossw1=12(002yx1+0+2w0x1+2w1x12)=12(2yx1+2w0x1+2w1x12)=122x1(y+w0+w1x1)=x1(y+y)=x1(yy)\frac{\partial loss}{\partial w_1} = \frac{1}{2}(0-0-2y*x_1+0+2 w_0 x_1 + 2 w_1 x_1^2) \\ = \frac{1}{2} (-2y x_1 + 2 w_0 x_1 + 2w_1 x_1^2) \\ = \frac{1}{2} * 2 x_1(-y + w_0 + w_1 x_1) \\ = x_1(-y + y') = - x_1(y - y')

推导完毕.