1.3 .波士顿房价预测任务1.3 波士顿房价预测任务波士顿房价预测任务上一节我们初步认识了神经网络的基本概念（如神

1.3 波士顿房价预测任务

样例如下

 0.00632  18.00   2.310  0  0.5380  6.5750  65.20  4.0900   1  296.0  15.30 396.90   4.98  24.00
 0.02731   0.00   7.070  0  0.4690  6.4210  78.90  4.9671   2  242.0  17.80 396.90   9.14  21.60
 0.02729   0.00   7.070  0  0.4690  7.1850  61.10  4.9671   2  242.0  17.80 392.83   4.03  34.70
 0.03237   0.00   2.180  0  0.4580  6.9980  45.80  6.0622   3  222.0  18.70 394.63   2.94  33.40
 0.06905   0.00   2.180  0  0.4580  7.1470  54.20  6.0622   3  222.0  18.70 396.90   5.33  36.20

上一节我们初步认识了神经网络的基本概念（如神经元、多层连接、前向计算、计算图）和模型结构三要素（模型假设、评价函数和优化算法）。本节将以“波士顿房价预测”任务为例，向读者介绍使用Python语言和Numpy库来构建神经网络模型的思考过程和操作方法。

波士顿房价预测是一个经典的机器学习任务，类似于程序员世界的“Hello World”。和大家对房价的普遍认知相同，波士顿地区的房价受诸多因素影响。该数据集统计了13种可能影响房价的因素和该类型房屋的均价，期望构建一个基于13个因素进行房价预测的模型，如图1 所示。

图1：波士顿房价影响因素示意图

对于预测问题，可以根据预测输出的类型是连续的实数值，还是离散的标签，区分为回归任务和分类任务。因为房价是一个连续值，所以房价预测显然是一个回归任务。下面我们尝试用最简单的线性回归模型解决这个问题，并用神经网络来实现这个模型。

1.3.1 线性回归模型

假设房价和各影响因素之间能够用线性关系来描述：

$y$ = $\sum_{j=1}^{M} x_j w_j +b$

模型的求解即是通过数据拟合出每个wjwj和bb。其中，wjwj和bb分别表示该线性模型的权重和偏置。一维情况下，wjwj 和 bb 是直线的斜率和截距。

线性回归模型使用均方误差作为（Mean Squared Error，MSE）损失函数（Loss），用以衡量预测房价和真实房价的差异，公式如下：

$MSE$ = $\frac{1}{n} \sum_{i=1}^{n} {y_i - \hat{y_i}}$

思考：

为什么要以均方误差作为损失函数？即将模型在每个训练样本上的预测误差加和，来衡量整体样本的准确性。这是因为损失函数的设计不仅仅要考虑“合理性”（有物理意义），同样需要考虑“易解性”（易于求解），这个问题在后面的内容中会详细阐述。

1.3.3 基于线性回归实现波士顿房价预测任务

深度学习不仅实现了模型的端到端学习，还推动了人工智能进入工业大生产阶段，产生了标准化、自动化和模块化的通用框架。不同场景的深度学习模型具备一定的通用性，五个步骤即可完成模型的构建和训练，如图3 所示。

图3：构建神经网络/深度学习模型的基本步骤

正是由于深度学习的建模和训练的过程存在通用性，即在构建不同的模型时，只有模型三要素不同，其它步骤基本一致，才产生了深度学习框架来加速建模。

1.3.3.1 数据处理

数据处理包含五个部分：数据导入、数据形状变换、数据集划分、数据归一化处理和封装load data函数。数据预处理后，才能被模型调用。

说明：

本教程中的代码都可以在AI Studio上直接运行，Print结果都是基于程序真实运行的结果。
由于是真实案例，代码之间存在依赖关系，因此需要读者逐条、全部运行，否则会导致命令执行报错。

（1）数据读取

通过如下代码读入数据，了解下波士顿房价的数据集结构，数据存放在本地目录下housing.data文件中。

# 导入需要用到的package
import numpy as np
import json
# 读入训练数据
datafile = './work/housing.data'
data = np.fromfile(datafile, sep=' ')
data

array([6.320e-03, 1.800e+01, 2.310e+00, ..., 3.969e+02, 7.880e+00,
       1.190e+01])

（2）数据形状变换

由于读入的原始数据是1维的，所有数据都连在一起。因此需要我们将数据的形状进行变换，形成一个2维的矩阵，每行为一个数据样本（14个值），每个数据样本包含13个 $x$ （影响房价的特征）和一个 $y$ （该类型房屋的均价）。

# 读入之后的数据被转化成1维array，其中array的第0-13项是第一条数据，第14-27项是第二条数据，以此类推.... 
# 这里对原始数据做reshape，变成N x 14的形式
feature_names = [ 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE','DIS', 
                 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV' ]
feature_num = len(feature_names)
data = data.reshape([data.shape[0] // feature_num, feature_num])

（3）数据集划分

将数据集划分成训练集和测试集，其中训练集用于确定模型的参数，测试集用于评判模型的效果。为什么要对数据集进行拆分，而不能直接应用于模型训练呢？这与学生时代的授课和考试关系比较类似，如图4 所示。

图4：训练集和测试集拆分的意义

上学时总有一些自作聪明的同学，平时不认真学习，考试前临阵抱佛脚，将习题死记硬背下来，但是成绩往往并不好。因为学校期望学生掌握的是知识，而不仅仅是习题本身。另出新的考题，才能鼓励学生努力去掌握习题背后的原理。同样我们期望模型学习的是任务的本质规律，而不是训练数据本身，模型训练未使用的数据，才能更真实的评估模型的效果。

在本案例中，我们将80%的数据用作训练集，20%用作测试集，实现代码如下。通过打印训练集的形状，可以发现共有404个样本，每个样本含有13个特征和1个预测值。

ratio = 0.8
offset = int(data.shape[0] * ratio)
training_data = data[:offset]
training_data.shape

(404, 14)

（4）数据归一化处理

对每个特征进行归一化处理，使得每个特征的取值缩放到0~1之间。这样做有两个好处：一是模型训练更高效；二是特征前的权重大小可以代表该变量对预测结果的贡献度（因为每个特征值本身的范围相同）。

# 计算train数据集的最大值，最小值
maximums, minimums = \
                     training_data.max(axis=0), \
                     training_data.min(axis=0), 
# 对数据进行归一化处理
for i in range(feature_num):
    data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])

（5）封装成load data函数

将上述几个数据处理操作封装成load data函数，以便下一步模型的调用，实现方法如下。


def load_data():
    # 从文件导入数据
    datafile = './work/housing.data'
    data = np.fromfile(datafile, sep=' ')

    # 每条数据包括14项，其中前面13项是影响因素，第14项是相应的房屋价格中位数
    feature_names = [ 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', \
                      'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV' ]
    feature_num = len(feature_names)

    # 将原始数据进行Reshape，变成[N, 14]这样的形状
    data = data.reshape([data.shape[0] // feature_num, feature_num])

    # 将原数据集拆分成训练集和测试集
    # 这里使用80%的数据做训练，20%的数据做测试
    # 测试集和训练集必须是没有交集的
    ratio = 0.8
    offset = int(data.shape[0] * ratio)
    training_data = data[:offset]

    # 计算训练集的最大值，最小值
    maximums, minimums = training_data.max(axis=0), \
                            training_data.min(axis=0)

    # 对数据进行归一化处理
    for i in range(feature_num):
        data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])

    # 训练集和测试集的划分比例
    training_data = data[:offset]
    test_data = data[offset:]
    return training_data, test_data

# 获取数据
training_data, test_data = load_data()
x = training_data[:, :-1]
y = training_data[:, -1:]

# 查看数据
print(x[0])
print(y[0])

[0.         0.18       0.07344184 0.         0.31481481 0.57750527
 0.64160659 0.26920314 0.         0.22755741 0.28723404 1.
 0.08967991]
[0.42222222]

1.3.3.2 模型设计¶

模型设计是深度学习模型关键要素之一，也称为网络结构设计，相当于模型的假设空间，即实现模型“前向计算”（从输入到输出）的过程。

如果将输入特征和输出预测值均以向量表示，输入特征𝑥x有13个向量，𝑦y有1个向量，那么参数权重的形状是13×113×1。假设我们以如下任意数字赋值参数做初始化：

𝑤=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,−0.1,−0.2,−0.3,−0.4,0.0]


w = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.1, -0.2, -0.3, -0.4, 0.0]
w = np.array(w).reshape([13, 1])

取出第1条样本数据，观察样本的特征向量与参数向量相乘的结果。


x1=x[0]
t = np.dot(x1, w)
print(t)

[0.69474855]

完整的线性回归公式，还需要初始化偏移量 $b$ ，同样随意赋初值-0.2。那么，线性回归模型的完整输出是 $z=t+b$ ，这个从特征和参数计算输出值的过程称为“前向计算”。


b = -0.2
z = t + b
print(z)

将上述计算预测输出的过程以“类和对象”的方式来描述，类成员变量有参数 $w$ 和 $b$ 。通过写一个forward函数（代表“前向计算”）完成上述从特征和参数到输出预测值的计算过程，代码实现如下。


class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，
        # 此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z

基于Network类的定义，模型的计算过程如下所示。

net = Network(13)
x1 = x[0]
y1 = y[0]
z = net.forward(x1)
print(z)

[2.39362982]

1.3.3.3 训练配置

模型设计完成后，需要通过训练配置寻找模型的最优值，即通过损失函数来衡量模型的好坏。训练配置也是深度学习模型关键要素之一。

通过模型计算 $x_1$ 表示的影响因素所对应的房价应该是 $z$ , 但实际数据告诉我们房价是 $y$ 。这时我们需要有某种指标来衡量预测值 $z$ 跟真实值 $y$ 之间的差距。对于回归问题，最常采用的衡量方法是使用均方误差作为评价模型好坏的指标，公式为

$Loss = (y - z)^2$

上式中的 $Loss$ 通常也被称作损失函数，它是衡量模型好坏的指标。在回归问题中常用均方误差作为损失函数，而在分类问题中常用采用交叉熵（Cross-Entropy）作为损失函数，在后续的章节中会更详细的介绍。对其中任意一个样本计算损失函数值的代码实现如下：

Loss = (y1 - z)*(y1 - z)
print(Loss)

因为计算损失函数时需要把每个样本的损失函数值都考虑到，所以我们需要对单个样本的损失函数进行求和，并除以样本总数 $N$ 。公式为

$Loss= \frac{1}{N}\sum_{i=1}^N{(y_i - z_i)^2}$

在Network类下面添加损失函数的代码实现如下：

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    def loss(self, z, y):
        error = z - y
        cost = error * error
        cost = np.mean(cost)
        return cost

使用定义的Network类，可以方便的计算预测值和损失函数。需要注意的是，类中的变量 $\mathbf{x}$ , $w$ ， $b$ , $z$ , $error$ 等均是向量。以变量 $\mathbf{x}$ 为例，共有两个维度，一个代表特征数量（值为13），一个代表样本数量，代码如下所示。

net = Network(13)
# 此处可以一次性计算多个样本的预测值和损失函数
x1 = x[0:3]
y1 = y[0:3]
z = net.forward(x1)
print('predict: ', z)
loss = net.loss(z, y1)
print('loss:', loss)

predict:  [[2.39362982]
 [2.46752393]
 [2.02483479]]
loss: 3.384496992612791

1.3.3.4 训练过程

上述计算过程描述了如何构建神经网络，通过神经网络完成预测值和损失函数的计算。接下来介绍如何求解参数 $w$ 和 $b$ 的数值，这个过程也称为模型训练过程。训练过程是深度学习模型的关键要素之一，其目标是让定义的损失函数尽可能的小，也就是说找到一个参数解 $w$ 和 $b$ ，使得损失函数取得极小值。

我们先做一个小测试：如图5 所示，基于微积分知识，求一条曲线在某个点的斜率等于函数在该点的导数值。那么大家思考下，当处于曲线的极值点时，该点的斜率是多少？

图5：曲线斜率等于导数值

这个问题并不难回答，处于曲线极值点时的斜率为0，即函数在极值点的导数为0。那么，让损失函数取极小值的 $w$ 和 $b$ 应该是下述方程组的解： $\frac{\partial{L}}{\partial{w}}=0$ $\frac{\partial{L}}{\partial{b}}=0$

将样本数据 $(x, y)$ 带入上面的方程组中即可求解出 $w$ 和 $b$ 的值，但是这种方法只对线性回归这样简单的任务有效。如果模型中含有非线性变换，或者损失函数不是均方差这种简单的形式，则很难通过上式求解。为了解决这个问题，下面我们将引入更加普适的数值求解方法：梯度下降法。

（1）梯度下降法

在现实中存在大量的函数正向求解容易，但反向求解较难，被称为单向函数，这种函数在密码学中有大量的应用。密码锁的特点是可以迅速判断一个密钥是否是正确的(已知 $x$ ，求 $y$ 很容易)，但是即使获取到密码锁系统，无法破解出正确的密钥是什么（已知 $y$ ，求 $x$ 很难）。神经网络模型的损失函数就是这样的单向函数，反向求解并不容易。

这种情况特别类似于一位想从山峰走到坡谷的盲人，他看不见坡谷在哪（无法逆向求解出 $Loss$ 导数为0时的参数值），但可以伸脚探索身边的坡度（当前点的导数值，也称为梯度）。那么，求解Loss函数最小值可以这样实现：从当前的参数取值，一步步的按照下坡的方向下降，直到走到最低点。这种方法笔者称它为“盲人下坡法”。哦不，有个更正式的说法“梯度下降法（Gradient Descent，GD）”。

训练的关键是找到一组 $(w, b)$ ，使得损失函数 $L$ 取极小值。我们先看一下损失函数 $L$ 只随两个参数 $w_5$ 、 $w_9$ 变化时的简单情形，启发下寻解的思路。 $L=L(w_5, w_9)$ 这里我们将 $w_0, w_1, ..., w_{12}$ 中除 $w_5, w_9$ 之外的参数和 $b$ 都固定下来，可以用图画出 $L(w_5, w_9)$ 的形式。

net = Network(13)
losses = []
#只画出参数w5和w9在区间[-160, 160]的曲线部分，以及包含损失函数的极值
w5 = np.arange(-160.0, 160.0, 1.0)
w9 = np.arange(-160.0, 160.0, 1.0)
losses = np.zeros([len(w5), len(w9)])

#计算设定区域内每个参数取值所对应的Loss
for i in range(len(w5)):
    for j in range(len(w9)):
        net.w[5] = w5[i]
        net.w[9] = w9[j]
        z = net.forward(x)
        loss = net.loss(z, y)
        losses[i, j] = loss

#使用matplotlib将两个变量和对应的Loss作3D图
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)

w5, w9 = np.meshgrid(w5, w9)

ax.plot_surface(w5, w9, losses, rstride=1, cstride=1, cmap='rainbow')
plt.show()

对于这种简单情形，我们利用上面的程序，可以在三维空间中画出损失函数随参数变化的曲面图。从图中可以看出有些区域的函数值明显比周围的点小。

需要说明的是：为什么这里我们选择 $w_5$ 和 $w_9$ 来画图？这是因为选择这两个参数的时候，可比较直观的从损失函数的曲面图上发现极值点的存在。其他参数组合，从图形上观测损失函数的极值点不够直观。

观察上述曲线呈现出“圆滑”的坡度，这正是我们选择以均方误差作为损失函数的原因之一。图6 呈现了只有一个参数维度时，均方误差和绝对值误差（只将每个样本的误差累加，不做平方处理）的损失函数曲线图。

图6：均方误差和绝对值误差损失函数曲线图

由此可见，均方误差表现的“圆滑”的坡度有两个好处：

曲线的最低点是可导的。
越接近最低点，曲线的坡度逐渐放缓，有助于通过当前的梯度来判断接近最低点的程度（是否逐渐减少步长，以免错过最低点）。

而绝对值误差是不具备这两个特性的，这也是损失函数的设计不仅仅要考虑“合理性”，还要追求“易解性”的原因。

现在我们要找出一组 $[w_5, w_9]$ 的值，使得损失函数最小，实现梯度下降法的方案如下：

步骤1：随机的选一组初始值，例如： $[w_5, w_9] = [-100.0, -100.0]$
步骤2：选取下一个点 $[w_5^{'} , w_9^{'}]$ ，使得 $L(w_5^{'} , w_9^{'}) < L(w_5, w_9)$
步骤3：重复步骤2，直到损失函数几乎不再下降。

如何选择 $[w_5^{'} , w_9^{'}]$ 是至关重要的，第一要保证 $L$ 是下降的，第二要使得下降的趋势尽可能的快。微积分的基础知识告诉我们，沿着梯度的反方向，是函数值下降最快的方向，如图7 所示。简单理解，函数在某一个点的梯度方向是曲线斜率最大的方向，但梯度方向是向上的，所以下降最快的是梯度的反方向。

图7：梯度下降方向示意图

**（2）计算梯度

上面我们讲过了损失函数的计算方法，这里稍微改写，为了使梯度计算更加简洁，引入因子 $\frac{1}{2}$ ，定义损失函数如下：

$L= \frac{1}{2N}\sum_{i=1}^N{(y_i - z_i)^2}$

其中 $z_i$ 是网络对第 $i$ 个样本的预测值：

$z_i = \sum_{j=0}^{12}{x_i^{j}\cdot w_j} + b$

梯度的定义：

$𝑔𝑟𝑎𝑑𝑖𝑒𝑛𝑡 = (\frac{\partial{L}}{\partial{w_0}},\frac{\partial{L}}{\partial{w_1}}, ... ,\frac{\partial{L}}{\partial{w_{12}}} ,\frac{\partial{L}}{\partial{b}})$

可以计算出 $L$ 对 $w$ 和 $b$ 的偏导数：

$\frac{\partial{L}}{\partial{w_j}} = \frac{1}{N}\sum_{i=1}^N{(z_i - y_i)\frac{\partial{z_i}}{\partial{w_j}}} = \frac{1}{N}\sum_{i=1}^N{(z_i - y_i)x_i^{j}}$

$\frac{\partial{L}}{\partial{b}} = \frac{1}{N}\sum_{i=1}^N{(z_i - y_i)\frac{\partial{z_i}}{\partial{b}}} = \frac{1}{N}\sum_{i=1}^N{(z_i - y_i)}$

从导数的计算过程可以看出，因子 $\frac{1}{2}$ 被消掉了，这是因为二次函数求导的时候会产生因子 $2$ ，这也是我们将损失函数改写的原因。

下面我们考虑只有一个样本的情况下，计算梯度：

$L= \frac{1}{2}{(y_i - z_i)^2}$

$z_1 = {x_1^{0}\cdot w_0} + {x_1^{1}\cdot w_1} + ... + {x_1^{12}\cdot w_{12}} + b$

可以计算出：

$L= \frac{1}{2}{({x_1^{0}\cdot w_0} + {x_1^{1}\cdot w_1} + ... + {x_1^{12}\cdot w_{12}} + b - y_1)^2}$

可以计算出 $L$ 对 $w$ 和 $b$ 的偏导数：

$\frac{\partial{L}}{\partial{w_0}} = ({x_1^{0}\cdot w_0} + {x_1^{1}\cdot w_1} + ... + {x_1^{12}\cdot w_12} + b - y_1)\cdot x_1^{0}=({z_1} - {y_1})\cdot x_1^{0}$

$\frac{\partial{L}}{\partial{b}} = ({x_1^{0}\cdot w_0} + {x_1^{1}\cdot w_1} + ... + {x_1^{12}\cdot w_{12}} + b - y_1)\cdot 1 = ({z_1} - {y_1})$

可以通过具体的程序查看每个变量的数据和维度。

x1 = x[0]
y1 = y[0]
z1 = net.forward(x1)
print('x1 {}, shape {}'.format(x1, x1.shape))
print('y1 {}, shape {}'.format(y1, y1.shape))
print('z1 {}, shape {}'.format(z1, z1.shape))

x1 [0.         0.18       0.07344184 0.         0.31481481 0.57750527 0.64160659 0.26920314 0.         0.22755741 0.28723404 1. 0.08967991], shape (13,)
y1 [0.42222222], shape (1,)
z1 [130.86954441], shape (1,)

按上面的公式，当只有一个样本时，可以计算某个 $w_j$ ，比如 $w_0$ 的梯度。

gradient_w0 = (z1 - y1) * x1[0]
print('gradient_w0 {}'.format(gradient_w0))

gradient_w0 [0.]

同样我们可以计算 $w_1$ 的梯度。


gradient_w1 = (z1 - y1) * x1[1]
print('gradient_w1 {}'.format(gradient_w1))

gradient_w1 [23.48051799]

依次计算 $w_2$ 的梯度。

gradient_w2= (z1 - y1) * x1[2]
print('gradient_w2 {}'.format(gradient_w2))

gradient_w2 [9.58029163]

1.3 .波士顿房价预测任务