最优化方法Python计算：无约束优化应用——线性回归预测器

一、线性回归预测器

用于预测的回归模型，指的是标签数据取连续实数值，回归函数是连续的。回归模型经训练得到最优模型参数${\mathbf{w}}_0$后，可用一组带有标签的数据集$({\mathbf{x}}_i^{\top },y_i),i=1,2,\cdots ,m$测试训练效果。记测试数据$\mathbf{y}=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right)$，$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}_1^{\top },1\\ {\mathbf{x}}_2^{\top },1\\ ⋮\\ {\mathbf{x}}_m^{\top },1\end{array}\right)$，预测值
$$\left(\begin{array}{c}{y}_{t_1}\\ y_{t_2}\\ ⋮\\ y_{t_m}\end{array}\right)={\mathbf{y}}_t=F({\mathbf{w}}_0;\mathbf{X})$$
与真实标签值$\mathbf{y}$之间一般说来是有差别的，测试的目的就是衡量这一差别的大小。记$y_{\text{mean}}=\frac{1}{m}\sum _{i=1}^m{y}_i$，$u=\sum _{i=1}^m(y_i-y_{t_i}{)}^2=\Vert \mathbf{y}-{\mathbf{y}}_t{\Vert }^2$，$v=\sum _{i=1}^m(y_i-y_{\text{mean}}{)}^2=\Vert \mathbf{y}-y_{\text{mean}}{\Vert }^2$，考虑比值$\frac{u}{v}$。

• 若模型的预测值$y_{t_i}$，$i=1,2,\cdots ,m$远离$y_{\text{mean}}$，即构建的模型不理想，则$\frac{u}{v}>1$。此时，$1-\frac{u}{v}<0$；
• 若预测值$y_{t_i}$，$i=1,2,\cdots ,m$聚集在$y_{\text{mean}}$附近，则$0\le 1-\frac{u}{v}\le 1$；
• 在情形(2)下，若若预测值$y_{t_i}=y_i$，$i=1,2,\cdots ,m$，即每个预测值与对应的标签真值一致。此时$u=0$，故$1-\frac{u}{v}=1$；
• 在情形(2)下，$y_{t_i}=y_{\text{mean}}$，$i=1,2,\cdots ,m$，即构建的是一个取值为$y_{\text{mean}}$的恒等模型。此时，$\frac{u}{v}=1$，故$1-\frac{u}{v}=0$。
  由此可见，$1-\frac{u}{v}$越接近1，模型训练得越成功。我们就以称为$R^2$测定系数的
  $$\text{score}=1-\frac{u}{v}$$
  作为指标来评估模型训练优劣。

import numpy as np
identity = lambda x: x							#恒等函数
class Regression():								#预测模型
    def score(self, x, y):						#测试函数
        yp = self.predict(x)					#预测值
        u = np.linalg.norm(y - yp) ** 2
        v = np.linalg.norm(y - y.mean()) ** 2
        return 1.0 - u / v
class LinearRegressor(Regression, LineModel):	#线性回归预测器
    def __init__(self):
        self.tagVal = identity

程序中

• 第2行定义了恒等函数identity。该函数的参数x表示输入数据，返回值为x。
• 第3~4行定义了Regression类。该类中实现了用于预测的测试函数score。该函数的两个参数x和y表示测试数据，x为样本特征数据，y为样本标签。第5行调用预测函数predict，计算对应x的预测值yp。第6、7行分别计算$\Vert \mathbf{y}-{\mathbf{y}}_t{\Vert }^2$和$\Vert \mathbf{y}-y_{\text{mean}}{\Vert }^2$赋予u和v。第8行计算评估分score并返回。
• 第9~11行定义了LinearRegressor类。该类继承了Regression和LineModel两个类，因此拥有这两个类的所有属性与函数。包括拟合函数F、训练函数fit、预测函数predict和测试函数score，以及在fit函数中产生的诸如A、y、scalar、w0等各项属性数据。类定义体中第10~11行声明的构造函数__init__将LinearRegressor类对象的标签值函数设置为恒等函数identity（第2行定义）。因为作为线性回归，对新的特征数据$\mathbf{x}$对应的标签预测值直接由式
  $$y=y\cdot (\max y-\min y)+\min y$$
  给出，因此执行predict函数时，返回
  $$\text{self.tagVal(yp * (self.ymax - self.ymin) + self.ymin)}$$
  （见博文《最优化方法Python计算：无约束优化应用——线性回归模型》中predict函数的定义）即可满足要求。

二、综合案例

文件SeoulBikeData.csv来自UC Irvine Machine Learning Repository，包含了某城市一共享单车运营公司从2017年12月1日0时至2018年11月30日23时共8760条记录（每行表示一条记录）。数据集包含日期信息（第0列）、每小时租赁的自行车数量（第1列）、时间（第2列）和气象信息：温度、湿度、风速、能见度、体感温度、太阳辐射、降雪量、降雨量（第3~10列）以及第11~13列的季节、假日和运营三个字段，其数据需数字化：约定

• 季节：Spring:0，Summer:1，Autumn:2，Winter:3；
• 假日：No Holiday:0，Holiday:1；
• 运营：No:0，Yes:1。

日期	租赁数	时间	温度	湿度	风速	能见度	体感温度	太阳辐射	降雨量	降雪量	季节	假日	运营
01/12/2017	254	0	-5.2	37	2.2	2000	-17.6	0	0	0	Winter	No Holiday	Yes
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
01/01/2018	206	0	-3.2	40	0.5	1358	-14.9	0	0	0	Winter	Holiday	Yes
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
13/05/2018	579	9	14.8	95	1	986	14	0.53	0	0	Spring	No Holiday	Yes
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
09/11/2018	0	0	12	96	3.1	1185	11.3	0	18	0	Autumn	No Holiday	No
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
30/11/2018	584	23	1.9	43	1.3	1909	-9.3	0	0	0	Autumn	No Holiday	Yes

下列代码读取数据，并按以上规定作数据类型的转换。此外，数据表中第0行是表示各项数据名称的表头，应予剔除。表中表示一小时单车租赁数的第1列为标签数据，应单列。表示日期的第0列数据在预测中并不适用，也应剔除。考虑到体感温度与气温和湿度相关，可忽略这一特征数据，故应将第7列剔除。

import numpy as np												#导入numpy
data = np.loadtxt('SeoulBikeData.csv', delimiter=',', dtype=str)#读取数据文件
X = np.array(data)												#转换为数组
X = X[1:,:]														#去掉表头
Y = X[:,1].astype(float)										#读取标签数据
X = X[:,[2,3,4,5,6,8,9,10,11,12,13]]							#去掉日期、租赁数、体感温度
X1 = X[:,8]														#读取季节列
X1 = np.array([0 if x == 'Spring' else							#串转换为数值
               1 if x == 'Summer' else
               2 if x == 'Autumn' else
               3 for x in X1])
X[:,8] = X1														#回填季节列
X1 = X[:,9]														#读取假日列
X1 = np.array([0 if x == 'No Holiday' else						#串转换为数值
               1 for x in X1])
X[:,9] = X1														#回填假期列
X1 = X[:,10]													#读取运营列
X1 = np.array([0 if x == 'No' else								#串转换为数值
               1 for x in X1])
X[:,10] = X1													#回填运营列
X = X.astype(float)												#所有数据转换为实数型
print('共享单车特征数据：')
print(X)
print('共享单车标签数：')
print(Y)
m=X.shape[0]													#读取样本个数
print('共有%d个数据样本'%m)

程序的第2行以串类型读取数据集文件SeoulBikeData.csv为data。第3行将data转换为Numpy数组X。第4行将X中表示表头的第0行去掉。第5行将表中第1列（租赁数）转存为表示标签数据的数组Y。第6行去掉原表中的第0列（日期）和第1列（租赁数，已转存），第7列体感温度，保留其他特征数据列。第7~12行将表中表示季节数据的第8列数据转换为数值数据‘Spring’转换为0，‘Summer’转换为1，‘Autumn’转换为2，‘Winter’转换为3。相仿地，第13~16行转换假日列数据，第17~20行转换运营列数据。第21行将表示样本特征数据的X中所有数据转为浮点型。运行程序，将输出

共享单车特征数据：
[[ 0.  -5.2 37.  ...  3.   0.   1. ]
 [ 1.  -5.5 38.  ...  3.   0.   1. ]
 [ 2.  -6.  39.  ...  3.   0.   1. ]
 ...
 [21.   2.6 39.  ...  2.   0.   1. ]
 [22.   2.1 41.  ...  2.   0.   1. ]
 [23.   1.9 43.  ...  2.   0.   1. ]]
 共享单车标签数：
[254. 204. 173. ... 694. 712. 584.]
共有8760个数据样本

下面，我们将X、Y中的数据随机地分成两部分：X[train]、Y[train]为训练数据中的特征及标签，X[test]、Y[test]为测试数据中的特征与标签。然后用X[train]和Y[train]训练线性回归模型，并以X[test]、Y[test]作为测试数据，计算R² 评估分。在上述程序后添加下列代码完成训练和测试

……
a = np.arange(m)									#下标集
np.random.seed(2026)								#随机种子
train = np.random.choice(a, m//3, replace = False)	#随机取得三分之一作为训练集下标
test = np.setdiff1d(a,train)						#测试集下标
print('随机选取%d个样本训练线性回归模型。'%(m//3))
bikeSharing = LinearRegressor()						#构造线性回归模型
bikeSharing.fit(X[train], Y[train])					#训练模型
print('系数：', bikeSharingDemand.coef_)
print('截距：%.4f'%bikeSharingDemand.intercept_)
score = bikeSharingDemand.score(X[test], Y[test])	#测试模型
print('对其余%d个样本数据测试，模型评估分：%.4f'%(m-m//3, score))

程序第1行的省略号表示前段程序的代码。第2行构造序列 { 0 , 1 , ⋯   , 8759 } \{0,1,\cdots,8759\} {0,1,⋯,8759}为a。第4行从a中随机选取三分之一（m//3）数据，作为训练数据的下标train。第5行从a中除去train的数据，剩下部分作为测试数据下标test。6行创建LinearRegressor类对象bikeSharing。第7行调用fit函数，用X[train]，Y[train]训练bikeSharingDmand。第8~9行输出训练所得系数和截距。第10行调用bikeSharing的测试函数score对测试特征数据Xtest和Ytest进行测试，第11行输出模型的R² 评分。运行程序，输出

训练中...，稍候
20次迭代后完成训练。
系数：[ 28.45 37.69 -5.13 1.88 0.05 -7.59
       -80.30 -72.07 7.22 -23.41 -145.0 838.5]
截距：-551.1319
对其余5840个样本数据测试，模型评估分：0.5341

即模型经过20次迭代训练完毕，由所得模型参数算得系数：

• 忽略其他特征，时间每增加1个小时，共享单车租用量约增加28台；
• 忽略其他特征，温度每增加1\textcelsius，共享单车租用量约增加37台；
• 忽略其他特征，湿度每增加1$\%$，共享单车租用量约减少5台；
• 忽略其他特征，风速每增加1m/s，共享单车租用量约增加2台；
• 忽略其他特征，能见度每增加10m，共享单车租用量几乎无变化；
• 忽略其他特征，体感温度每增加1\textcelsius，共享单车租用量约减少8台；
• 忽略其他特征，太阳辐射每增加1MJ/m${}^2$，共享单车租用量约减少80台；
• 忽略其他特征，降雨量每增加1cm，共享单车租用量约减少72台；
• 忽略其他特征，降雪量每增加1cm，共享单车租用量约增加7台；
• 忽略其他特征，季节每变化1个（按春、夏、秋、冬顺序），共享单车租用量约减少23台；
• 忽略其他特征，假日共享单车租用量约减少145台； \item 忽略其他特征，公司运营单车租用量约增加838台。
  基本符合人们的预期。截距意为若所有特征归零，单车租赁数约减少551台。训练所得模型的评估分为0<0.5341<1，属正常模型（R² 分数越接近 1，模型拟合效果越好）。

写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！