HBU_神经网络与深度学习 实验3 线性回归

---

写在前面的一些内容

1. 本文为HBU_神经网络与深度学习实验（2022年秋）实验3的实验报告，此文的基本内容参照 [1]神经网络与深度学习：案例与实践 - 第2章（上）：线性回归理论解读 和 [2]NNDL 实验2（下），检索时请按对应题号进行检索。
2. 本实验编程语言为Python 3.10，使用Pycharm进行编程。
3. 本实验报告目录标题级别顺序：一、1. （1）
4. 水平有限，难免有误，如有错漏之处敬请指正。

---

一、实现一个简单的线性回归模型

进行实验前，我们需要先导入torch包。

import torch

---

1. 数据集构建

首先，我们构造一个小的回归数据集。假设输入特征和输出标签的维度都为1，那么我们来定义一个被拟合的函数。

def linefunc(x, w=1.8, b=0.5):
    y = w * x + b
    return y

接下来的内容都将以$y=1.8x+0.5$ 这个被拟合函数为基准进行操作。

# create_data.py
def create_data(func, interval, sample_num, noise, add_outlier=False, outlier_ratio=0.001):
    """
    根据给定的函数，生成样本
    输入：
       - func：函数
       - interval： x的取值范围
       - sample_num： 样本数目
       - noise： 噪声均方差
       - add_outlier：是否生成异常值
       - outlier_ratio：异常值占比
    输出：
       - X: 特征数据，shape=[n_samples,1]
       - y: 标签数据，shape=[n_samples,1]
    """
    # 均匀采样
    # 使用paddle.rand在生成sample_num个随机数
    X = torch.rand(sample_num) * (interval[1] - interval[0]) + interval[0]
    y = func(X)
    # 生成高斯分布的标签噪声
    # 使用paddle.normal生成概率分布为N~(0,noise^2)的数据
    epsilon = torch.normal(0, noise, y.shape)
    y = y + epsilon
    if add_outlier:  # 生成额外的异常点
        outlier_num = int(len(y) * outlier_ratio)
        if outlier_num != 0:
            # 使用paddle.randint生成服从均匀分布的、范围在[0, len(y))的随机张量
            outlier_idx = torch.randint(len(y), [outlier_num])
            y[outlier_idx] = y[outlier_idx] * 5
    return X, y

利用上面的生成样本函数，生成150个带噪音的样本，其中100个训练样本，50个测试样本，并打印出训练数据的可视化分布。

from matplotlib import pyplot as plt
from create_data import create_data

def data_processing(func, interval, train_num, test_num, step, noise):
    X_train, y_train = create_data(func=func, interval=interval, sample_num=train_num, noise=noise, add_outlier=False)
    X_test, y_test = create_data(func=func, interval=interval, sample_num=test_num, noise=noise, add_outlier=False)
    # paddle.linspace返回一个张量，张量的值为在区间start和stop上均匀间隔的num个值，输出张量的长度为num
    X_underlying = torch.linspace(interval[0], interval[1], step)
    y_underlying = func(X_underlying)
    return X_train, y_train, X_test, y_test, X_underlying, y_underlying

func = linefunc
interval = (-10, 10)  # x的取值范围
train_num, test_num = 100, 50  # 训练样本数目和测试样本数目
noise = 3  # 标准差
step = 100
X_train, y_train, X_test, y_test, X_underlying, y_underlying = data_processing(func, interval, train_num, test_num, step, noise)

# 绘制数据
plt.scatter(X_train, y_train, marker='*', facecolor="none", edgecolor='#9932CC', s=50, label="train data") # 训练数据标点颜色为暗紫色 (darkorchid)
plt.scatter(X_test, y_test, facecolor="none", edgecolor='#00BFFF', s=50, label="test data") # 测试数据标点颜色为深天蓝色 (deepskyblue)
plt.plot(X_underlying, y_underlying, c='#000000', label=r"underlying distribution")
plt.legend(fontsize='x-large')  # 给图像加图例
plt.savefig('ml-vis.pdf') # 保存图像到PDF文件中
plt.show()

代码执行结果如下图所示：

这里的data_processing中的内容在本次实验后续内容还会再次使用，所以将其写成函数。

---

2. 模型构建

在线性回归中，自变量为样本的特征向量$\mathbf{x}\in {\mathbb{R}}^D$（每一维对应一个自变量），因变量是连续值的标签$y\in R$。
线性模型定义为：
$$\begin{array}{r}f(\mathbf{x};\mathbf{w},b)={\mathbf{w}}^T\mathbf{x}+b\end{array}$$其中权重向量$\mathbf{w}\in {\mathbb{R}}^D$和偏置$b\in \mathbb{R}$都是可学习的参数。[1]

在实践中，为了提高预测样本的效率，我们通常会将$N$ 样本归为一组进行成批地预测，这样可以更好地利用GPU设备的并行计算能力。
$$\begin{array}{r}\mathbf{y}=\mathbf{X}\mathbf{w}+b\end{array}$$其中$\mathbf{X}\in {\mathbb{R}}^{N\times D}$ 为$N$ 个样本的特征矩阵，$\mathbf{y}\in {\mathbb{R}}^N$ 为$N$ 个预测值组成的列向量。[1]

好了，接下来我们先定义一个算子（op.py）：

# op.py
import torch
torch.manual_seed(10)  # 设置随机种子

class Op(object):
    def __init__(self):
        pass

    def __call__(self, inputs):
        return self.forward(inputs)

    def forward(self, inputs):
        raise NotImplementedError

    def backward(self, inputs):
        raise NotImplementedError

# 线性算子
class Linear(Op):
    def __init__(self, input_size):
        """
        输入：
           - input_size:模型要处理的数据特征向量长度
        """
        self.input_size = input_size
        # 模型参数
        self.params = {}
        self.params['w'] = torch.randn(self.input_size, 1)
        self.params['b'] = torch.zeros([1])

    def __call__(self, X):
        return self.forward(X)

    # 前向函数
    def forward(self, X):
        """
        输入：
           - X: tensor, shape=[N,D]
           注意这里的X矩阵是由N个x向量的转置拼接成的，与原教材行向量表示方式不一致
        输出：
           - y_pred： tensor, shape=[N]
        """
        N, D = X.shape
        if self.input_size == 0:
            return torch.full(size=[N, 1], fill_value=self.params['b'])
        assert D == self.input_size  # 输入数据维度合法性验证

        # 使用torch.matmul计算两个tensor的乘积
        y_pred = torch.matmul(X, self.params['w']) + self.params['b']
        return y_pred

然后我们来构建一个线性回归模型：

import op

# 注意这里我们为了和后面章节统一，这里的X矩阵是由N个x向量的转置拼接成的，与原教材行向量表示方式不一致
input_size = 3
N = 2
X = torch.randn(N, input_size)  # 生成2个维度为3的数据
model = op.Linear(input_size)
y_pred = model(X)
print("y_pred:", y_pred)  # 输出结果的个数也是2个

代码执行结果：

y_pred: tensor([[1.8529],
        [0.6011]])

---

3. 损失函数

回归任务是对连续值的预测，希望模型能根据数据的特征输出一个连续值作为预测值。因此回归任务中常用的评估指标是均方误差。
令$\mathbf{y}\in {\mathbb{R}}^N$，$\hat{\mathbf{y}}\in {\mathbb{R}}^N$分别为$N$ 个样本的真实标签和预测标签，均方误差的定义为：
$$\begin{array}{r}\mathcal{L}(\mathbf{y},\hat{\mathbf{y}})=\frac{1}{2N}\Vert \mathbf{y}-\hat{\mathbf{y}}{\Vert }^2=\frac{1}{2N}\Vert \mathbf{X}\mathbf{w}+\mathbf{b}-\mathbf{y}{\Vert }^2\end{array}$$其中$\mathbf{b}$ 为$N$ 维向量，所有元素取值都为$b$。[1]

首先我们需要写一个均方误差函数的代码：

# mean_squared_error.py
import torch

def mean_squared_error(y_true, y_pred):
    assert y_true.shape[0] == y_pred.shape[0]
    # torch.square 计算输入的平方值
    # torch.mean 沿 axis 计算 x 的平均值，默认axis是None, 则对输入的全部元素计算平均值。
    error = torch.mean(torch.square(y_true - y_pred))
    return error

然后我们构造一个样例来测试一下：

from mean_squared_error import mean_squared_error

# 构造一个简单的样例进行测试:[N,1], N=2
y_true = torch.tensor([[-0.2], [4.9]], dtype=torch.float32)  # tensor, 样本真实标签
y_pred = torch.tensor([[1.3], [2.5]], dtype=torch.float32)  # tensor, 样本预测标签
error = mean_squared_error(y_true=y_true, y_pred=y_pred).item()  # error: float, 误差值
print("error:", error)

代码执行结果：

error: 4.005000114440918

---

4. 模型优化

参数学习的过程通过优化器完成。由于我们可以基于最小二乘方法可以直接得到线性回归的解析解，此处的训练是求解析解的过程，代码实现如下：

# optimizer.py
import torch

def optimizer_lsm(model, X, y, reg_lambda=0):
    """
        输入：
           - model: 模型
           - X: tensor, 特征数据，shape=[N,D]
           - y: tensor,标签数据，shape=[N]
           - reg_lambda: float, 正则化系数，默认为0
        输出：
           - model: 优化好的模型
    """
    N, D = X.shape

    # 对输入特征数据所有特征向量求平均
    x_bar_tran = torch.mean(X, axis=0).T

    # 求标签的均值,shape=[1]
    y_bar = torch.mean(y)

    # torch.subtract通过广播的方式实现矩阵减向量
    x_sub = torch.subtract(X, x_bar_tran)

    # 使用torch.all判断输入tensor是否全0
    if torch.all(x_sub == 0):
        model.params['w'] = torch.zeros(size=[D])
        model.params['b'] = y_bar
        return model

    # torch.inverse求方阵的逆
    tmp = torch.inverse(torch.matmul(x_sub.T, x_sub) + reg_lambda * torch.eye(D))

    w = torch.matmul(torch.matmul(tmp, x_sub.T), (y - y_bar))
    b = y_bar - torch.matmul(x_bar_tran, w)

    model.params['w'] = torch.squeeze(w, axis=-1)
    model.params['b'] = b

    return model

---

5. 模型训练

在准备了数据、模型、损失函数和参数学习的实现之后，我们开始模型的训练。在回归任务中，模型的评价指标和损失函数一致，都为均方误差。
通过上文实现的线性回归类来拟合训练数据，并输出模型在训练集上的损失。

from opitimizer import optimizer_lsm

input_size = 1
model = op.Linear(input_size)
model = optimizer_lsm(model, X_train.reshape([-1, 1]), y_train.reshape([-1, 1]))
print("w_pred:", model.params['w'].item(), "b_pred: ", model.params['b'].item())

y_train_pred = model(X_train.reshape([-1, 1])).squeeze()
train_error = mean_squared_error(y_true=y_train, y_pred=y_train_pred).item()
print("train error: ", train_error)

X_train_large, y_train_large = create_data(func=func, interval=interval, sample_num=5000, noise = noise, add_outlier = False)

model_large = op.Linear(input_size)
model_large = optimizer_lsm(model_large,X_train_large.reshape([-1,1]),y_train_large.reshape([-1,1]))
print("w_pred large:",model_large.params['w'].item(), "b_pred large: ", model_large.params['b'].item())

y_train_pred_large = model_large(X_train_large.reshape([-1,1])).squeeze()
train_error_large = mean_squared_error(y_true=y_train_large, y_pred=y_train_pred_large).item()
print("train error large: ",train_error_large)

代码执行结果：

w_pred: 1.7244662046432495 b_pred:  0.10675737261772156
train error:  9.735129356384277
w_pred large: 1.7934911251068115 b_pred large:  0.48613476753234863
train error large:  8.730573654174805

显然，$w$ 和$b$ 的预测值与真实值$(1.8，0.5)$ 有一定的差距。

---

6. 模型评估

下面用训练好的模型预测一下测试集的标签，并计算在测试集上的损失。

y_test_pred = model(X_test.reshape([-1, 1])).squeeze()
test_error = mean_squared_error(y_true=y_test, y_pred=y_test_pred).item()
print("test error: ", test_error)

y_test_pred_large = model_large(X_test.reshape([-1,1])).squeeze()
test_error_large = mean_squared_error(y_true=y_test, y_pred=y_test_pred_large).item()
print("test error large: ",test_error_large)

代码执行结果：

test error:  8.030173301696777
test error large:  6.8764328956604

---

还没完事，咱们来加一点训练数据的样本数量，再加点正则化系数的调整。

（1）调整训练数据的样本数量

修改的部分代码如下所示：

train_num, test_num = 5000, 50  # 训练样本数目和测试样本数目

代码执行结果：

w_pred: 1.8004437685012817 b_pred:  0.45358312129974365
train error:  8.924786567687988
w_pred large: 1.7934911251068115 b_pred large:  0.48613476753234863
train error large:  8.730573654174805
test error:  8.784930229187012
test error large:  8.785351753234863

相比于原训练样本数目，新的训练样本数目更多，计算得到的$w$ 和$b$ 和被拟合函数的更相近。

（2）调整正则化系数

修改的部分代码如下所示：

def optimizer_lsm(model, X, y, reg_lambda=5):

代码执行结果：

w_pred: 1.7161484956741333 b_pred:  0.4196764826774597
train error:  10.837532043457031
w_pred large: 1.7934348583221436 b_pred large:  0.4861411154270172
train error large:  8.730573654174805
test error:  9.686600685119629
test error large:  9.726728439331055

正则化系数调整后对训练没什么影响，这说明训练模型没有发生过拟合。

---

二、多项式回归

多项式回归是回归任务的一种形式，其中自变量和因变量之间的关系是$M$ 次多项式的一种线性回归形式，即：
$$\begin{array}{r}f(\mathbf{x};\mathbf{w})=w_{1}x+w_2{x}^2+...+w_M{x}^M+b={\mathbf{w}}^T\varphi (x)+b\end{array}$$其中$M$ 为多项式的阶数，$\mathbf{w}=[w_1,...,w_M{]}^T$ 为多项式的系数，$\varphi (x)=[x,x^2,\cdots ,x^M{]}^T$ 为多项式基函数，将原始特征$x$ 映射为$M$ 维的向量。当$M=0$ 时，$f(\mathbf{x};\mathbf{w})=b$。
公式$(4)$ 展示的是特征维度为1的多项式表达，当特征维度大于1时，存在不同特征之间交互的情况，这是线性回归无法实现。公式$(5)$ 展示的是当特征维度为2，多项式阶数为2时的多项式回归：
$$\begin{array}{r}f(\mathbf{x};\mathbf{w})=w_1{x}_1+w_2{x}_2+w_3{x}_1^2+w_4{x}_1{x}_2+w_5{x}_2^2+b\end{array}$$[1]

接下来我们基于特征维度为1的自变量介绍多项式回归实验。

---

1. 数据集构建

假设我们要拟合的非线性函数为一个缩放后的$sin$ 函数。

# sin函数: sin(2 * pi * x)
def sin(x):
    y = torch.sin(2 * math.pi * x)
    return y

构建训练和测试数据的方式和前文的方式相同，其中训练数样本15个，测试样本10个，高斯噪声标准差为0.5，自变量范围为 (0,1)。

# 生成数据
func = sin
interval = (0, 1)
train_num, test_num = 15, 10
noise = 0.5  # 0.1
step = 100
X_train, y_train, X_test, y_test, X_underlying, y_underlying = data_processing(func, interval, train_num, test_num, step, noise)

# 绘制图像
plt.rcParams['figure.figsize'] = (8.0, 6.0)
plt.scatter(X_train, y_train, marker='*', facecolor="none", edgecolor='#9932CC', s=50, label="train data")
# plt.scatter(X_test, y_test, facecolor="none", edgecolor="r", s=50, label="test data")
plt.plot(X_underlying, y_underlying, c='#00BFFF', label=r"$\sin(2\pi x)$")
plt.legend(fontsize='x-large')
plt.savefig('ml-vis2.pdf')
plt.show()

代码执行结果如下图所示：

在输出结果中，绿色的曲线是周期为1的$sin$ 函数曲线，蓝色的圆圈为生成的训练样本数据，红色的圆圈为生成的测试样本数据。

---

2. 模型构建

通过多项式的定义可以看出，多项式回归和线性回归一样，同样学习参数$\mathbf{w}$，只不过需要对输入特征$\varphi (x)$ 根据多项式阶数进行变换。因此，我们可以套用求解线性回归参数的方法来求解多项式回归参数。

def polynomial_basis_function(x, degree=2):
    """
    输入：
       - x: tensor, 输入的数据，shape=[N,1]
       - degree: int, 多项式的阶数
       example Input: [[2], [3], [4]], degree=2
       example Output: [[2^1, 2^2], [3^1, 3^2], [4^1, 4^2]]
       注意：本案例中,在degree>=1时不生成全为1的一列数据；degree为0时生成形状与输入相同，全1的Tensor
    输出：
       - x_result： tensor
    """
    if degree == 0:
        return torch.ones(size=x.shape, dtype=torch.float32)
    x_tmp = x
    x_result = x_tmp
    for i in range(2, degree + 1):
        x_tmp = torch.multiply(x_tmp, x)  # 逐元素相乘
        x_result = torch.concat((x_result, x_tmp), axis=-1)
    return x_result

# 简单测试
data = [[2], [3], [4]]
X = torch.tensor(data=data, dtype=torch.float32)
degree = 3
transformed_X = polynomial_basis_function(X, degree=degree)
print("转换前:", X)
print("阶数为", degree, "转换后:", transformed_X)

代码执行结果：

转换前: tensor([[2.],
        [3.],
        [4.]])
阶数为 3 转换后: tensor([[ 2.,  4.,  8.],
        [ 3.,  9., 27.],
        [ 4., 16., 64.]])

---

3. 模型训练

对于多项式回归，我们可以同样使用前面线性回归中定义的LinearRegression算子、训练函数train、均方误差函数mean_squared_error。拟合训练数据的目标是最小化损失函数，同线性回归一样，也可以通过矩阵运算直接求出$\mathbf{w}$ 的值。
我们设定不同的多项式阶，$M$ 的取值分别为0、1、3、8，之前构造的训练集上进行训练，观察样本数据对$sin$ 曲线的拟合结果。

from opitimizer import optimizer_lsm

plt.rcParams['figure.figsize'] = (12.0, 8.0)

for i, degree in enumerate([0, 1, 3, 8]):  # []中为多项式的阶数
    model = op.Linear(degree)
    X_train_transformed = polynomial_basis_function(X_train.reshape([-1, 1]), degree)
    X_underlying_transformed = polynomial_basis_function(X_underlying.reshape([-1, 1]), degree)
    model = optimizer_lsm(model, X_train_transformed, y_train.reshape([-1, 1]))  # 拟合得到参数
    y_underlying_pred = model(X_underlying_transformed).squeeze()
    print(model.params)
    # 绘制图像
    plt.subplot(2, 2, i + 1)
    plt.scatter(X_train, y_train, marker='*', facecolor="none", edgecolor='#9370DB', s=50, label="train data")
    plt.plot(X_underlying, y_underlying, c='#000000', label=r"$\sin(2\pi x)$")
    plt.plot(X_underlying, y_underlying_pred, c='#FF0000', label="predicted function")
    plt.ylim(-2, 1.5)
    plt.annotate("M={}".format(degree), xy=(0.95, -1.4))

# plt.legend(bbox_to_anchor=(1.05, 0.64), loc=2, borderaxespad=0.)
plt.legend(loc='lower left', fontsize='x-large')
plt.savefig('ml-vis3.pdf')
plt.show()

代码执行结果：

{'w': tensor([0.]), 'b': tensor(-0.0696)}
{'w': tensor([-2.3304]), 'b': tensor([1.1818])}
{'w': tensor([ 21.8905, -54.1687,  34.9021]), 'b': tensor([-1.7652])}
{'w': tensor([ 102.5296, -319.9999,  344.7184,   40.8362, -252.0008,   55.0874,
         -10.3919,   51.4386]), 'b': tensor([-10.1180])}

执行代码后得到下图：

观察可视化结果，红色的曲线表示不同阶多项式分布拟合数据的结果：
当$M=0$ 或$M=1$ 时，拟合曲线较简单，模型欠拟合；
当$M=8$ 时，拟合曲线较复杂，模型过拟合；
当$M=3$ 时，模型拟合最为合理。

---

4. 模型评估

下面通过均方误差来衡量训练误差、测试误差以及在没有噪音的加入下$sin$ 函数值与多项式回归值之间的误差，更加真实地反映拟合结果。多项式分布阶数从0到8进行遍历。

from mean_squared_error import mean_squared_error

# 训练误差和测试误差
training_errors = []
test_errors = []
distribution_errors = []

# 遍历多项式阶数
for i in range(9):
    model = op.Linear(i)

    X_train_transformed = polynomial_basis_function(X_train.reshape([-1, 1]), i)
    X_test_transformed = polynomial_basis_function(X_test.reshape([-1, 1]), i)
    X_underlying_transformed = polynomial_basis_function(X_underlying.reshape([-1, 1]), i)

    optimizer_lsm(model, X_train_transformed, y_train.reshape([-1, 1]))

    y_train_pred = model(X_train_transformed).squeeze()
    y_test_pred = model(X_test_transformed).squeeze()
    y_underlying_pred = model(X_underlying_transformed).squeeze()

    train_mse = mean_squared_error(y_true=y_train, y_pred=y_train_pred).item()
    training_errors.append(train_mse)

    test_mse = mean_squared_error(y_true=y_test, y_pred=y_test_pred).item()
    test_errors.append(test_mse)

    # distribution_mse = mean_squared_error(y_true=y_underlying, y_pred=y_underlying_pred).item()
    # distribution_errors.append(distribution_mse)

print("train errors: \n", training_errors)
print("test errors: \n", test_errors)
# print ("distribution errors: \n", distribution_errors)

# 绘制图片
plt.rcParams['figure.figsize'] = (8.0, 6.0)
plt.plot(training_errors, '-.', mfc="none", mec='#FF0000', ms=10, c='#FF0000', label="Training")
plt.plot(test_errors, '--', mfc="none", mec='#FF69B4', ms=10, c='#FF69B4', label="Test")
# plt.plot(distribution_errors, '-', mfc="none", mec="#3D3D3F", ms=10, c="#3D3D3F", label="Distribution")
plt.legend(fontsize='x-large')
plt.xlabel("degree")
plt.ylabel("MSE")
plt.savefig('ml-mse-error.pdf')
plt.show()

代码执行结果：

train errors: 
 [0.6441656351089478, 0.45997682213783264, 0.3104371428489685, 0.20104871690273285, 0.19899645447731018, 0.20385366678237915, 0.8440933227539062, 0.31182458996772766, 0.2052197903394699]
test errors: 
 [0.6824741363525391, 0.4967922568321228, 1.1924924850463867, 0.18429550528526306, 0.44991421699523926, 2.477367877960205, 4.155946254730225, 2.193247079849243, 1.981374740600586]

执行代码后得到下图：

观察可视化结果：

• 当阶数较低的时候，模型的表示能力有限，训练误差和测试误差都很高，代表模型欠拟合；
• 当阶数较高的时候，模型表示能力强，但将训练数据中的噪声也作为特征进行学习，一般情况下训练误差继续降低而测试误差显著升高，代表模型过拟合。

对于模型过拟合的情况，可以引入正则化方法，通过向误差函数中添加一个惩罚项来避免系数倾向于较大的取值。下面加入$l_2$ 正则化项，查看拟合结果。

degree = 8  # 多项式阶数
reg_lambda = 0.0001  # 正则化系数

X_train_transformed = polynomial_basis_function(X_train.reshape([-1, 1]), degree)
X_test_transformed = polynomial_basis_function(X_test.reshape([-1, 1]), degree)
X_underlying_transformed = polynomial_basis_function(X_underlying.reshape([-1, 1]), degree)

model = op.Linear(degree)
optimizer_lsm(model, X_train_transformed, y_train.reshape([-1, 1]))
y_test_pred = model(X_test_transformed).squeeze()
y_underlying_pred = model(X_underlying_transformed).squeeze()

model_reg = op.Linear(degree)
optimizer_lsm(model_reg, X_train_transformed, y_train.reshape([-1, 1]), reg_lambda=reg_lambda)
y_test_pred_reg = model_reg(X_test_transformed).squeeze()
y_underlying_pred_reg = model_reg(X_underlying_transformed).squeeze()

mse = mean_squared_error(y_true=y_test, y_pred=y_test_pred).item()
print("mse:", mse)
mes_reg = mean_squared_error(y_true=y_test, y_pred=y_test_pred_reg).item()
print("mse_with_l2_reg:", mes_reg)

# 绘制图像
plt.scatter(X_train, y_train, marker='*', facecolor="none", edgecolor="#9370DB", s=50, label="train data")
plt.plot(X_underlying, y_underlying, c='#000000', label=r"$\sin(2\pi x)$")
plt.plot(X_underlying, y_underlying_pred, c='#FF0000', linestyle="--", label="$deg. = 8$")
plt.plot(X_underlying, y_underlying_pred_reg, c='#FF69B4', linestyle="-.", label="$deg. = 8, \ell_2 reg$")
plt.ylim(-1.5, 1.5)
plt.annotate("lambda={}".format(reg_lambda), xy=(0.82, -1.4))
plt.legend(fontsize='large')
plt.savefig('ml-vis4.pdf')
plt.show()

代码执行结果：

mse: 1.981374740600586
mse_with_l2_reg: 0.1406262069940567

执行代码后得到下图：

观察可视化结果，其中黄色曲线为加入$l_2$ 正则后多项式分布拟合结果，红色曲线为未加入$l_2$ 正则的拟合结果，黄色曲线的拟合效果明显好于红色曲线。

---

三、Runner类介绍

通过上面的实践，我们可以看到，在一个任务上应用机器学习方法的流程基本上包括：数据集构建、模型构建、损失函数定义、优化器、模型训练、模型评价、模型预测等环节。
为了更方便地将上述环节规范化，我们将机器学习模型的基本要素封装成一个Runner类。除上述提到的要素外，再加上模型保存、模型加载等功能。
Runner类的成员函数定义如下表所示：[2]

函数	定义
__ init __	实例化Runner类时默认调用，需要传入模型、损失函数、优化器和评价指标等
train	完成模型训练，指定模型训练需要的训练集和验证集
evaluate	通过对训练好的模型进行评价，在验证集或测试集上查看模型训练效果
predict	选取一条数据对训练好的模型进行预测
save_model	模型在训练过程和训练结束后需要进行保存
load_model	调用加载之前保存的模型

Runner类的框架定义如下：

# Runner.py
import os

class Runner(object):
    def __init__(self, model, optimizer, loss_fn, metric):
        # 优化器和损失函数为None,不再关注
        # 模型
        self.model = model
        # 评估指标
        self.metric = metric
        # 优化器
        self.optimizer = optimizer

    def train(self, dataset, reg_lambda, model_dir):
        X, y = dataset
        self.optimizer(self.model, X, y, reg_lambda)
        # 保存模型
        self.save_model(model_dir)

    def evaluate(self, dataset, **kwargs):
        X, y = dataset
        y_pred = self.model(X)
        result = self.metric(y_pred, y)
        return result

    def predict(self, X, **kwargs):
        return self.model(X)

    def save_model(self, model_dir):
        if not os.path.exists(model_dir):
            os.makedirs(model_dir)
        params_saved_path = os.path.join(model_dir, 'params.pdtensor')
        torch.save(self.model.params, params_saved_path)

    def load_model(self, model_dir):
        params_saved_path = os.path.join(model_dir, 'params.pdtensor')
        self.model.params = torch.load(params_saved_path)

Runner类的流程如上图所示，可以分为4个阶段：
初始化阶段：传入模型、损失函数、优化器和评价指标。
模型训练阶段：基于训练集调用$train()$ 函数训练模型，基于验证集通过$evaluate()$ 函数验证模型。通过$save$$\_model()$ 函数保存模型。
模型评价阶段：基于测试集通过$evaluate()$ 函数得到指标性能。
模型预测阶段：给定样本，通过$predict()$ 函数得到该样本标签。[2]

---

四、基于线性回归的波士顿房价预测

在本节中，我们使用线性回归来对马萨诸塞州波士顿郊区的房屋进行预测。实验流程主要包含如下5个步骤：

• 数据处理：包括数据清洗（缺失值和异常值处理）、数据集划分，以便数据可以被模型正常读取，并具有良好的泛化性;
• 模型构建：定义线性回归模型类；
• 训练配置:训练相关的一些配置，如：优化算法、评价指标等；
• 组装训练框架Runner:Runner用于管理模型训练和测试过程；
• 模型训练和测试:利用Runner进行模型训练和测试。
  进行实验前，我们需要先导入pandas包。

import pandas as pd  # 开源数据分析和操作工具

---

1. 数据处理

本实验使用波士顿房价预测数据集，共506条样本数据，每条样本包含了12种可能影响房价的因素和该类房屋价格的中位数。预览前5条数据，代码实现如下：

（1）数据预览

# 利用pandas加载波士顿房价的数据集
data = pd.read_csv("boston_house_prices.csv")
# 预览前5行数据
print(data.head())

代码执行结果：

      CRIM    ZN  INDUS  CHAS    NOX  ...  RAD  TAX  PTRATIO  LSTAT  MEDV
0  0.00632  18.0   2.31     0  0.538  ...    1  296     15.3   4.98  24.0
1  0.02731   0.0   7.07     0  0.469  ...    2  242     17.8   9.14  21.6
2  0.02729   0.0   7.07     0  0.469  ...    2  242     17.8   4.03  34.7
3  0.03237   0.0   2.18     0  0.458  ...    3  222     18.7   2.94  33.4
4  0.06905   0.0   2.18     0  0.458  ...    3  222     18.7   5.33  36.2

（2）数据清洗

• 缺失值分析

对数据集中的缺失值或异常值等情况进行分析和处理，保证数据可以被模型正常读取。
通过isna()方法判断数据中各元素是否缺失，然后通过sum()方法统计每个字段缺失情况。

# 查看各字段缺失值统计情况
print(data.isna().sum())

代码执行结果：

[5 rows x 13 columns]
CRIM       0
ZN         0
INDUS      0
CHAS       0
NOX        0
RM         0
AGE        0
DIS        0
RAD        0
TAX        0
PTRATIO    0
LSTAT      0
MEDV       0
dtype: int64

从输出结果看，波士顿房价预测数据集中不存在缺失值的情况。

• 异常值处理

通过箱线图直观的显示数据分布，并观测数据中的异常值。
箱线图一般由五个统计值组成：最大值、上四分位、中位数、下四分位和最小值。
一般来说，观测到的数据大于最大估计值或者小于最小估计值则判断为异常值，其中

最大估计值=上四分位+1.5∗(上四分位−下四分位) 最小估计值=下四分位−1.5∗(上四分位−下四分位)

def boxplot(data, fig_name):
    # 绘制每个属性的箱线图
    data_col = list(data.columns)

    # 连续画几个图片
    plt.figure(figsize=(5, 5), dpi=300)
    # 子图调整
    plt.subplots_adjust(wspace=0.6)
    # 每个特征画一个箱线图
    for i, col_name in enumerate(data_col):
        plt.subplot(3, 5, i + 1)
        # 画箱线图
        plt.boxplot(data[col_name],
                    showmeans=True,
                    meanprops={"markersize": 1, "marker": "D", "markeredgecolor": '#f19ec2'},  # 均值的属性
                    medianprops={"color": '#e4007f'},  # 中位数线的属性
                    whiskerprops={"color": '#e4007f', "linewidth": 0.4, 'linestyle': "--"},
                    flierprops={"markersize": 0.4},
                    )
        # 图名
        plt.title(col_name, fontdict={"size": 5}, pad=2)
        # y方向刻度
        plt.yticks(fontsize=4, rotation=90)
        plt.tick_params(pad=0.5)
        # x方向刻度
        plt.xticks([])
    plt.savefig(fig_name)
    plt.show()
    
boxplot(data, 'ml-vis5.pdf')

代码执行结果如下图所示：

箱线图示例如下图所示。

从输出结果看，数据中存在较多的异常值（图中上下边缘以外的空心小圆圈）。

使用四分位值筛选出箱线图中分布的异常值，并将这些数据视为噪声，其将被临界值取代，代码实现如下：

# 四分位处理异常值
num_features = data.select_dtypes(exclude=['object', 'bool']).columns.tolist()

for feature in num_features:
    if feature == 'CHAS':
        continue

    Q1 = data[feature].quantile(q=0.25)  # 下四分位
    Q3 = data[feature].quantile(q=0.75)  # 上四分位

    IQR = Q3 - Q1
    top = Q3 + 1.5 * IQR  # 最大估计值
    bot = Q1 - 1.5 * IQR  # 最小估计值
    values = data[feature].values
    values[values > top] = top  # 临界值取代噪声
    values[values < bot] = bot  # 临界值取代噪声
    data[feature] = values.astype(data[feature].dtypes)

# 再次查看箱线图，异常值已被临界值替换（数据量较多或本身异常值较少时，箱线图展示会不容易体现出来）
boxplot(data, 'ml-vis6.pdf')

代码执行结果如下图所示：

从输出结果看，经过异常值处理后，箱线图中异常值得到了改善。

（3）数据集划分

import torch
torch.manual_seed(10)

# 划分训练集和测试集
def train_test_split(X, y, train_percent=0.8):
    n = len(X)
    shuffled_indices = torch.randperm(n)  # 返回一个数值在0到n-1、随机排列的1-D Tensor
    train_set_size = int(n * train_percent)
    train_indices = shuffled_indices[:train_set_size]
    test_indices = shuffled_indices[train_set_size:]

    X = X.values
    y = y.values

    X_train = X[train_indices]
    y_train = y[train_indices]

    X_test = X[test_indices]
    y_test = y[test_indices]

    return X_train, X_test, y_train, y_test

X = data.drop(['MEDV'], axis=1)
y = data['MEDV']

X_train, X_test, y_train, y_test = train_test_split(X,y)# X_train每一行是个样本，shape[N,D]

（4）特征工程

为了消除纲量对数据特征之间影响，在模型训练前，需要对特征数据进行归一化处理，将数据缩放到[0, 1]区间内，使得不同特征之间具有可比性。

# 特征工程
X_train = torch.tensor(X_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.float32)
y_test = torch.tensor(y_test, dtype=torch.float32)

X_min = torch.min(X_train, dim=0)[0]
X_max = torch.max(X_train, dim=0)[0]

X_train = (X_train - X_min) / (X_max - X_min)
X_test = (X_test - X_min) / (X_max - X_min)

# 训练集构造
train_dataset = (X_train, y_train)
# 测试集构造
test_dataset = (X_test, y_test)

进行X_train-X_min时，为避免报错，我们需要min和max的值（Value），加入索引[0]即可。

---

2. 模型构建

实例化一个线性回归模型，特征维度为 12:

from nndl.op import Linear

# 模型实例化
input_size = 12
model=Linear(input_size)

---

3. 完善Runner类

---

4. 模型训练

在组装完成Runner之后，我们将开始进行模型训练、评估和测试。首先，我们先实例化Runner，然后开始进行装配训练环境，接下来就可以开始训练了，相关代码如下：

# 模型保存文件夹
saved_dir = './models'

from Runner import Runner
from opitimizer import optimizer_lsm

optimizer = optimizer_lsm
# 实例化Runner
runner = Runner(model, optimizer=optimizer, loss_fn=None, metric=mse_loss)
# 启动训练
runner.train(train_dataset, reg_lambda=0, model_dir=saved_dir)

打印出训练得到的权重：

columns_list = data.columns.to_list()
weights = runner.model.params['w'].tolist()
b = runner.model.params['b'].item()

for i in range(len(weights)):
    print(columns_list[i], "weight:", weights[i])

print("b:", b)

代码执行结果：

CRIM weight: -5.2610182762146
ZN weight: 1.3627111911773682
INDUS weight: -0.024850845336914062
CHAS weight: 1.800213098526001
NOX weight: -7.556697368621826
RM weight: 9.557083129882812
AGE weight: -1.3511630296707153
DIS weight: -9.967939376831055
RAD weight: 7.528402328491211
TAX weight: -5.082475662231445
PTRATIO weight: -6.9966325759887695
LSTAT weight: -13.183658599853516
b: 32.62158966064453

从输出结果看，CRIM、PTRATIO等的权重为负数，表示该镇的人均犯罪率与房价负相关，学生与教师比例越大，房价越低。RAD和CHAS等为正，表示到径向公路的可达性指数越高，房价越高；临近Charles River房价高。

---

5. 模型测试

加载训练好的模型参数，在测试集上得到模型的MSE指标。

# 加载模型权重
runner.load_model(saved_dir)

mse = runner.evaluate(test_dataset)
print('MSE:', mse.item())

代码执行结果：

MSE: 11.210769653320312

---

6. 模型预测

使用Runner中load_model函数加载保存好的模型，使用predict进行模型预测，代码实现如下：

runner.load_model(saved_dir)
pred = runner.predict(X_test[:1])
print("真实房价：", y_test[:1].item())
print("预测的房价：", pred.item())

代码执行结果：

真实房价： 18.899999618530273
预测的房价： 21.529155731201172

从输出结果看，预测房价接近真实房价。

---

五、实验Q&A

问题1：使用类实现机器学习模型的基本要素有什么优点？

便于扩写和复用。

问题2：算子op、优化器opitimizer放在单独的文件中，主程序在使用时调用该文件。这样做有什么优点？

能够精简代码行数，使程序更简洁。

问题3：线性回归通常使用平方损失函数，能否使用交叉熵损失函数？为什么？

不能。平方损失函数还和错误的分类有关，但是这种调整在线性回归中没有必要。

---

六、实验总结

问题不小，错误不少，只能说之前欠下来的终究还是要还的。paddle和torch终归不能随意通用，还是要注意函数的使用啊。先把paddle-pytorch API对应表 混个眼熟吧。
马上就要考六级力（悲）
那就让暴风雨来的再猛烈一些吧！

下面贴一个normal()的用法：
TypeError: normal() received an invalid combination of arguments - got (float, int, Tensor), but expected one of:
* (Tensor mean, Tensor std, *, torch.Generator generator, Tensor out)
* (Tensor mean, float std, *, torch.Generator generator, Tensor out)
* (float mean, Tensor std, *, torch.Generator generator, Tensor out)
* (float mean, float std, tuple of ints size, *, torch.Generator generator, Tensor out, torch.dtype dtype, torch.layout layout, torch.device device, bool pin_memory, bool requires_grad)