从零开始的Python机器学习指南（二）——监督学习之OLS回归

介绍

本博客将结合样例介绍监督学习/Supervised Learning/SL下的第一大分支：回归/Regression。

开始前的准备

开始前，请先确保你的python环境中有以下包：
pandas，numpy，sklearn。

本文的所有代码可在Anaconda的Jupyter Lab里运行。

正文

首先来理解下：为什么回归/Regression是一种监督学习？这个问题的本质是什么？

我们首先要理解回归问题的本质。简单来说，回归问题的本质是，对于一个映射$f$,
$$f:{\mathbb{R}}^n\mapsto \mathbb{R}$$
，我们有它定义域/Domain和对应域/Codomain的一些数组。我们想要通过这些已有的数据来找到一个最符合这些数据的模型/回归函数。

换个角度看，对于每一条数据，该映射的定义域可以理解为特征集/Feature Set，其对应域可理解为标签集/Label Set。我们要找的就是未知的映射函数，对应着机器学习里的黑箱模型/Blackbox Model。

但是要把它变成一个完全的监督学习问题，我们还需要定义下我们要优化的目标是什么。既然我们要找到最符合数据的映射，我们可以考虑最小化误差平方和/Sum of Squared Errors/SSE。对于一个有$m$条数据的数据集和一个对应的映射$f$，该映射对于该数据集的误差平方和定义为

$$\ell =\sum _{i=1}^m(y^{(i)}-f({\mathbf{x}}^{(i)}){)}^2$$

，也就是把每一条数据的真实标签和映射的预测标签的差进行平方后求和。如果我们把误差平方和作为监督学习的损失函数/Loss Function，则该回归问题可被称为普通最小二乘法线性回归问题/OLS Linear Regression。我们要做的就是找到映射$f$，使得误差平方和最小。

我们把这个问题更加严谨地下个定义。让$\mathbf{x}$作为特征集，$\mathbf{w}$作为每个特征的权重集，则对于线性映射
$$f(\mathbf{x};\mathbf{w})={\mathbf{w}}^T\cdot \mathbf{x}=w_0{x}_0+w_1{x}_1+...+w_n{x}_n$$
，我们要最小化的损失函数$\ell$可定义为：
$$\ell (\mathbf{w})=\sum _{i=1}^m(y^{(i)}-f({\mathbf{x}}^{(i)}){)}^2=\sum _{i=1}^n(y^{(i)}-{\mathbf{w}}^T\cdot {\mathbf{x}}^{(i)}{)}^2=\sum _{i=1}^n(y^{(i)}-{\hat{y}}^{(i)}{)}^2.$$
也就是说，线性回归在监督学习下的本质可以理解为一个优化问题：
$$\mathbf{w}=\arg \underset{\mathbf{w}}{\min }\ell (\mathbf{w}).$$

那我们怎么能找到最合适的线性回归函数？
首先，这个问题是有唯一且最优解的，属于P/NP问题中的P。单从这一点来说，线性回归在机器学习界是相对简单和容易解决的问题了。

找到普通最小二乘法线性回归问题/OLS Linear Regression最优解的过程只需要依赖一些线性代数就能解决。推导过程较为繁琐和复杂，博主不在此展示，感兴趣的同学们可以自己搜索了解一下。但结论是，最符合数据（或者最小化误差平方和）的特征权重集${\mathbf{w}}^{\ast }$，可以通过以下矩阵乘法得出：
$${\mathbf{w}}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{Y}$$
，其中$\mathbf{X}$是特征集（大小为特征种类乘以数据行数），$\mathbf{Y}$是有$m$个元素的向量（标签集）。对于计算机来说，这种矩阵乘法是小case啦。

但是，没有学习的过程，还算机器学习吗？
很多同学，包括博主自己，在学习到这里的时候都有这个问题。模型似乎没有迭代和进步的过程，而是一步从零跳跃到最优解。但严格意义上来说，学习的过程是有的。模型获取了特征集和标签集，并且用这些数据提高了自己对数据的泛化能力。虽然没有一步一步走向更好的性能，但是电脑仍旧以有效的方式学习到了数据的特征。所以总而言之，线性回归也是一种机器学习。

对于一个有未知回归函数的数据集，如果我们通过上述方法取得了最优回归函数（即最小化误差平方和），那么我们对于新的数据就有一定的预测能力。

那线性回归有什么应用场景吗？
很多啊，比如：

1. 房价预测。如果你有过去房市的数据（比如年平均气温，日平均光照时间，离最近便利店的距离，等等）以及每年的房价，并且你如果认为房价和这些房间特征有（多维）线性关系，那么你可以用OLS来预测房产的价格。
2. 股票分析。类似房价，股票的一些特征（比如公司净利润，开盘时间，涨停次数等等）也是对于股价的重要指标。但要注意的是，由于时间和一些无法获得的数据（比如公司老板打算摆烂的概率）对股价的影响非常大，简单的线性回归模型可能不太适合预测股票这如此复杂的问题。
3. 顾客终身价值预测（CLV）。一个顾客给一个企业带来的收益总和很大程度上与顾客的收入水平、年龄、日均消费等信息挂钩。如果有顾客信息作为特征集和对应的对企业的总收益作为标签集，那么就可以用这些信息预测它们的CLV。

代码

理解了原理后，我们可以用python实现上述的回归学习算法。

首先，我们来生成一些简单的数据。

import numpy as np # 用来进行一些数学运算
import pandas as pd # 用来用数据框的方式储存数据
import matplotlib # 用来画图的
import matplotlib.pyplot as plt
from sklearn import linear_model # 我们需要的模型

# X和Y包含生成的数据。X可以理解为特征/Features，Y可以看作标签/Label
# 我们要做的就是训练一个模型，使它能够通过X准确预测Y的对应值
# coeff包含生成函数所使用的系数
X, Y, coeff = skd.make_regression(n_samples=5000, n_features=2, 
                                  noise=5.0, 
                                  coef=True,
                                  random_state=114514)

print(f"特征集大小为{X.shape}，标签集大小为{Y.shape}。")
# 特征集大小为(5000,2)，标签集大小为(5000,)。

在生成好数据后，我们可以检查下数据内容：

from mpl_toolkits.mplot3d import Axes3D # 一个用来画3D图形的库

fig = plt.figure()
plot3d = fig.gca(projection='3d') # 调用matplotlib的对象
plot3d.view_init(elev=15., azim=35) # 观察角度

x1 = X[:, 0] # 选中X里第0列的所有行
x2 = X[:, 1] # 选中X里第1列的所有行

# 画出数据散点。注意x1，x2，Y的行数要相同
plot3d.scatter3D(x1, x2, Y, c=Y, cmap='Greens', label='Noised Samples');

# 画出真实参考平面
x1min, x1max = int(np.floor(min(x1))), int(np.ceil(max(x1)))
x2min, x2max = int(np.floor(min(x2))), int(np.ceil(max(x2)))
x1plane = np.linspace(x1min, x1max, 2000)
x2plane = np.linspace(x2min, x2max, 2000)
xx1, xx2 = np.meshgrid(range(x1min, x1max), range(x2min, x2max))
x12_coord = np.array([xx1, xx2])
y = coeff[0] * x12_coord[0] + coeff[1] + x12_coord[1]
surf = plot3d.plot_surface(xx1, xx2, y, alpha=0.2, label="True model")

# 给3D图形加上标签
surf._facecolors2d = surf._facecolor3d
surf._edgecolors2d = surf._edgecolor3d 
plot3d.legend()

# 画出坐标轴标签
plot3d.set_xlabel('$x_1$')
plot3d.set_ylabel('$x_2$')
plot3d.set_zlabel('$y$')

plt.show()

我们得到了以下图形：

接下来就是用sklearn训练线性回归模型了。代码非常简单，只有两行。

# 创造线性回归对象
OLS = linear_model.LinearRegression(fit_intercept=True)

# 训练模型
OLS.fit(X, Y)

我们可以手动比较一下真实模型和训练模型的差距：

print('真实模型参数:')
print('\tw0: {:.2f}, w1: {:.2f}, w2: {:.2f}'.format(0, *coeff))
print('训练模型参数:')
print('\tw0: {:.2f}, w1: {:.2f}, w2: {:.2f}'.format(0, *OLS.coef_))
'''
真实模型参数:
	w0: 0.00, w1: 12.56, w2: 7.64
训练模型参数:
	w0: 0.00, w1: 12.58, w2: 7.77
'''

更加科学和全面的方法是用交叉验证/Cross-validation。这种方法的原理是：每次把数据随机分成训练集和验证集，并代入模型评估模型在训练和验证下的错误。数学上来说，这种方法缩小了随机取样验证的随机性，是对于模型性能更加全面的评估。代码如下：

from sklearn.model_selection import cross_validate

# 这里的cv指的是把数据分成几份。比如，cv=10就是把数据分成十份，取一份作为验证集，九份作为训练集
cv_results = cross_validate(OLS, X, Y, cv=10, scoring="r2",
                            return_train_score=True)

print('Mean test score: {:.3f} (std: {:.3f})'
      '\nMean train score: {:.3f} (std: {:.3f})'.format(
                                              np.mean(cv_results['test_score']),
                                              np.std(cv_results['test_score']),
                                              np.mean(cv_results['train_score']),
                                              np.std(cv_results['train_score'])))
'''
Mean test score: 0.894 (std: 0.009)
Mean train score: 0.895 (std: 0.001)
'''

我们的模型有89.4%的验证准确率和89.5%的训练准确率，总体来说是一个较为准确、没有过拟合也没有欠拟合的模型。

我们还可以用训练好的模型进行预测，比如：

print(OLS.predict([[42,24]])[0]) # 模型预测结果
print(coeff[0]*42+coeff[1]*24) # 真实值
'''
714.8098105630355
710.6836496552116
'''

可以观察到，虽然保证了最优解，但是模型的预测仍不是百分百准确的。这个现象的成因来源于数据集的随机噪声/Random Noise，也就是数据本身也没有百分百符合其潜在模型。有些时候我们可以减小噪声对算法的影响。比如，对于图像处理来说，我们可以使用高斯滤波/Gaussian Filter来减少图像噪点对学习算法的干扰。

结语

在下一篇博客博主会介绍如何用ML方法解决简单的分类/Classification问题。有任何问题和建议请随时评论或私信。码字不易，喜欢博主内容的话请点赞支持！