本文采用仅采用numpy进行房价预测,不包含任何深度学习框架,波士顿房价预测问题作为机器学习领域的“hello world”。本文有利于读者更好地理解并入门,为以后更深程度的学习打下良好的基础。
目录
前言
波士顿房价数据集源于美国某经济学杂志上,分析研究波士顿 房价的数据集。该数据集中的每一行数据都是对波士顿周边或者城镇 的房价情况的概述,统计的是房价的中位数和 13 个指标,试图能够 分析出这些指标和房价之间的关系。本文主要以数据预处理,损失函数,梯度计算,反向传播更新参数展开,具体介绍这几个步骤在代码中的体现。
一,数据处理:
读入数据集:
本文采用housing.data数据集。数据集的格式如图所示,该数据集共有506个样本,每个样本包含13个特征值和1个与之对应的标签。如图所示:
该数据集共有506行14列的数据,每一行数据代表一个样本(每个样本里包含14的数据),前13列分别对应13个特征值,最后一列为各个样本值对应的标签,即房价,所有数据均已空格分开。
datafile = 'housing.data'
data = np.fromfile(datafile, sep=' ')
其中'housing.data'为路径名,我的直接放在了本文件夹下。并赋值给datafile
np.fromfile表示读取该路径下的文件,sep = ''表示以空格分隔。此时data为一个1行7084列的向量。
数据维度转换:
由于读入的原始数据是1维的。因此需要我们将数据的形状进行变换,形成一个2维的矩阵,每行为一个数据样本(14个值),每个数据样本包含13个特征值(影响房价的特征)和一个对应的标签(该类型房屋的均价)。
# 读入之后的数据被转化成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])
feature_names为一个列表,内有14个值,前13个元素代表对应的特征向量,最后一个元素为标签。每个特征向量的含义如下图所示:
feature_number = len(feature_names),用len函数求出特征值的数量并赋值给feature_number,为14(包含标签),data.shape表示求data的维度为一维(7084,),data.shape[0]即为7084,所以data.shape[0]对feature_number(14)取余,得到506;所以data的维度为[506,14]。data为506行,14列的二维矩阵。
区分训练集和测试集,通常取总样本数的80%作为训练集,剩余的20%作为测试集。
代码如下:
ratio = 0.8
offset = int(data.shape[0] * ratio)
training_data = data[:offset]
test_data = data[offset:]
对data进行切片处理,offset为整形变量,代表506*0.8。data[:offset]是指对二维矩阵的第一维进行切片处理,第二维不变。得到的training_data和test_data分别为404行14列的矩阵和102行14列的矩阵。
数据归一化处理:
对每个特征进行归一化处理,使得每个特征的取值缩放到0~1之间。
这样做有两个好处:
一是模型训练更高效;
二是特征前的权重大小可以代表该变量对预测结果的贡献度(因为每个特征值本身的范围相同)
代码如下:
# 计算train数据集的最大值,最小值,平均值
maximums, minimums, avgs = \
training_data.max(axis=0), \
training_data.min(axis=0), \
training_data.sum(axis=0) / training_data.shape[0]
# 对数据进行归一化处理
for i in range(feature_num):
#print(maximums[i], minimums[i], avgs[i])
data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])
首先求出每一个特征值的最大值,最小值和平均值并分别赋值给maximums, minimums, avgs。即对所有的列分别求出每列的最大值,最小值和平均值。.max()函数可以求出最大值,但由于training_data是一个404行14列的矩阵,而我们分别求出每一列的最大值,最小值和平均值。所以要指定维度为0,即对列求。 (axis=0),若不指定,则对所有数据求最大值。.min()函数同理。 training_data.sum(axis=0)表示对每一列求和最后得到一个1行14列的矩阵,每一个数据为对应每一列所有行数之和。training_data.shape表示[404,14],其中0表示第一个元素。training_data.shape[0]表示404,二者之比自然为平均数。
之后的for循环表示依次遍历每一列,由于data是506行,14列的矩阵,所以循环里面用到切片操作data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])。其中data[:, i] 表示506行,第i列的矩阵,minimums[i]表示第i列的最小值,maximums[i]表示第i列的最大值。经过此处理后,便可使所有数据均在0-1之间。
代码封装成类:
import numpy as np
def load_data():
# 从文件导入数据
datafile = '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, avgs = training_data.max(axis=0), training_data.min(axis=0), \
training_data.sum(axis=0) / training_data.shape[0]
# 对数据进行归一化处理
for i in range(feature_num):
#print(maximums[i], minimums[i], avgs[i])
data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])
# 训练集和测试集的划分比例
training_data = data[:offset]
test_data = data[offset:]
return training_data, test_data
此操作仅为训练时直观有效,不封装成类也没有关系。但封装成类更简单清晰,我这是直接写在Tensor_3.py文件里。
Tensor_1为主程序,使用时直接从Tensor_3中调用load_data类即可。如图所示:
二,模型设计
线性回归:
模型设计是深度学习模型关键要素之一,也称为网络结构设计,相当于模型的假设空间,即实现模型“前向计算”(从输入到输出)的过程。本文只讨论线性回归,线性回归是利用数理统计中回归分析,来确定两种或两种以上 变量间相互依赖的定量关系的一种统计分析方法。它对一个或多个自 变量与因变量之间的关系进行建模,得出的函数是一个或多个称为回 归系数的模型参数的线性组合。用一个方程式来表示如下: 𝑦 = 𝜔𝑖𝑥𝑖 + 𝑏 𝑖 = 0,1, …
线性回归要解决的问题便是得出最佳的拟合线,建立自变量与 因变量之间的关系,获得该条最佳拟合线可以通过使用最小二乘法来 解决。对于观测数据,它通过最小化每一个数据点到线的垂直偏差平 方和来计算最佳拟合线。线性回归的目标函数如下: 𝑚𝑖𝑛 𝜔 ‖𝑋𝜔 − 𝑦‖ 2
前向传播:
将上述计算预测输出的过程以“类和对象”的方式来描述,类成员变量有参数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类,我们训练所用的方法都会集成在这个类里。
初始化num_of_weights权重的数量,即为特征值的数量,这里是13,因为我们的数据有14列,最后一列是标签不需要有权重w,然后随机初始化这些w,np.random.seed(0)这里采用随机数种子0
numpy.random.randn(d0,d1,…,dn)。randn函数返回一个或一组样本,具有标准正态分布。dn表格每个维度,返回值为指定维度的array,初始化b为0。
forword方法中,np.dot()表示矩阵相乘。x为任意行13列的矩阵,self.w为13行1列的矩阵。加上权重b,最后输出z为任意行1列的矩阵,实现前向传播。
损失函数:
模型设计完成后,需要通过训练配置寻找模型的最优值,即通过损失函数来衡量模型的好坏。训练配置也是深度学习模型关键要素之一。通过模型计算x1x_1x1表示的影响因素所对应的房价应该是zzz, 但实际数据告诉我们房价是yyy。这时我们需要有某种指标来衡量预测值z跟真实值y之间的差距。对于回归问题,最常采用的衡量方法是使用均方误差作为评价模型好坏的指标,具体定义如下:
Loss=(y−z)2Loss = (y - z)^2Loss=(y−z)2
上式中的LossLossLoss(简记为: L)通常也被称作损失函数,它是衡量模型好坏的指标。损失函数的定义并不唯一。
对一个样本计算损失函数值的实现如下:
def loss(self, z, y):
error = z - y
cost = error * error
cost = np.mean(cost)
return cost
z为前向传播求出的权重,为任意行1列的矩阵,为模型的预测值,由于还没有进行反向传播更新参数,w最开始的赋值的任意的,所以第一次得到的z值就是蒙的。相应的实际值y也为任意行1列的矩阵。error = z - y,将预测值矩阵z与实际值矩阵y做差,所得到的值赋给error。此时的error可能为正,也可能为负数。而cost为error的平方,对应数据相乘,并非矩阵的乘法,必定为正数,因此,我们通常会说,损失值cost越少,模型的质量越高,预测效果越好。np.mean()为求平均值函数,对任意行1列cost内的元素求均值,并赋值给cost,最后返回cost。
计算梯度:
所谓梯度,简单来说就是对应变量的导数。根据上面的结论,我们可以得到损失函数loss关于z,y的函数。而z是w和b函数,我们的最终目的是求出尽可能精确的w和b,从而使得预测值z无限逼近于y,使loss尽可能地减少。
由高等数学的知识我们可以知道,函数沿梯度方向上升最快。对应的,函数沿梯度反方向下降的最快,所以我们只要求出梯度取反就可以使loss下降最快,即loss分别对w和b求偏导。
以下为详细的数学推导部分,数学基础相对较好的读者可以仔细研读。数学基础相对较弱的读者可以直接跳到欢乐的代码环节。
代码实现如下:
def gradient(self, x, y):
z = self.forward(x)
gradient_w = (z-y)*x
gradient_w = np.mean(gradient_w, axis=0)
gradient_w = gradient_w[:, np.newaxis]
gradient_b = (z - y)
gradient_b = np.mean(gradient_b)
return gradient_w, gradient_b
gradient方法即为求梯度的方法,x为任意行13列的矩阵,代表任意组样本值,每一列代表对应的特征值,y为任意行1列的矩阵,代表任意行样本的实际值。gradient_w = (z-y)*x,其中 (z-y)*x为loss对w求偏导,gradient_w表示loss对w的梯度,为任意行13列。这是任意行样本值求出的偏导,用np.mean方法求任意组的均值axis=0表示对行求平均值,此时gradient_w为1维向量,长度为13
不符合模型输入要求,因此要单独加1个维度,gradient_w = gradient_w[:, np.newaxis]表示数据不变,单独加一个维度。
gradient_b同理。
反向传播:
反向传播更新参数,使之前的被随机初始化的w和b向更准确的权重靠近。
我们之前说过函数沿梯度的反方向下降最快,而梯度就是函数对变量的偏导数,我们既然已经求得对应变量的偏导数,函数的梯度自然已经得到。为使loss函数下降的很快,进而使模型更精确,我们需要更新这些“不太靠谱”的w和b。
代码如下所示:
def update(self, gradient_w, gradient_b, eta = 0.01):
self.w = self.w - eta * gradient_w
self.b = self.b - eta * gradient_b
eta表示学习率设置为0.01,也叫步长。简单来说,偏导数的方向只是一个正确的方向,但每次要向着这个方向迈出多大的步子,由学习率决定。
学习率设置过大,可能会错过正确值。
学习率设置过小,参数更新很慢。
封装函数:
目前,我们已经写好了训练需要的所有方法。通常我们会把这些方法连同训练过程一起封装到一个类里,利于操作。
代码如下:
def train(self, training_data, num_epochs, batch_size=10, eta=0.01):
n = len(training_data)
losses = []
for epoch_id in range(num_epochs):
# 在每轮迭代开始之前,将训练数据的顺序随机打乱
# 然后再按每次取batch_size条数据的方式取出
np.random.shuffle(training_data)
# 将训练数据进行拆分,每个mini_batch包含batch_size条的数据
mini_batches = [training_data[k:k + batch_size] for k in range(0, n, batch_size)]
for iter_id, mini_batch in enumerate(mini_batches):
# print(self.w.shape)
# print(self.b)
x = mini_batch[:, :-1]
y = mini_batch[:, -1:]
a = self.forward(x)
loss = self.loss(a, y)
gradient_w, gradient_b = self.gradient(x, y)
self.update(gradient_w, gradient_b, eta)
losses.append(loss)
print('Epoch {:3d} / iter {:3d}, loss = {:.4f}'.
format(epoch_id, iter_id, loss))
return losses
Epoch表示一轮迭代,一轮迭代表示所有的样本均被使用过一轮。
epoch_id表示第几轮迭代。
num_epochs表示迭代的轮数。
mini_batche表示一个batch所包含的最小样本数。
batch_size表示一个batch所包含样本数据的大小。
mini_batches表示batch的数量
x表示每个mini_batch的特征值
y表示每个mini_batch的标签
np.random.shuffle(training_data)表示将training_data乱序,进改变行的位置,样本数据不发生改变。
mini_batches = [training_data[k:k + batch_size] for k in range(0, n, batch_size)]
其中mini_batches为一个列表,列表里的每一个元素均为,batch_size行14列的矩阵。
外层循环表示迭代的轮数。
内层循环表示每次喂入网络batch_size行14列的矩阵,直到一个epoch结束。即所有数据均被使用一次。
关于为什么不只使用一次数据:
简单来说,就像我们人类学习一样,我们对后学的知识记忆程度相比较于之前学过的知识记得更清晰一些。在房价预测问题上,后出现的数据会对模型产生权重较大的影响。为了消除这一影响,通常会对训练集进行对次epoch,在此之前都要用shuffle进行乱序。
总结:
所有代码如图所示:
Tensor_3.py:
import numpy as np
def load_data():
# 从文件导入数据
datafile = '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, avgs = training_data.max(axis=0), training_data.min(axis=0), \
training_data.sum(axis=0) / training_data.shape[0]
# 对数据进行归一化处理
for i in range(feature_num):
#print(maximums[i], minimums[i], avgs[i])
data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])
# 训练集和测试集的划分比例
training_data = data[:offset]
test_data = data[offset:]
return training_data, test_data
Tensor_1.py:
from Tensor_3 import load_data
import numpy as np
import matplotlib.pyplot as plt
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
def gradient(self, x, y):
z = self.forward(x)
gradient_w = (z - y) * x
gradient_w = np.mean(gradient_w, axis=0)
gradient_w = gradient_w[:, np.newaxis]
gradient_b = (z - y)
gradient_b = np.mean(gradient_b)
return gradient_w, gradient_b
def update(self, gradient_w, gradient_b, eta=0.01):
self.w = self.w - eta * gradient_w
self.b = self.b - eta * gradient_b
def train(self, training_data, num_epochs, batch_size=10, eta=0.01):
n = len(training_data)
losses = []
for epoch_id in range(num_epochs):
# 在每轮迭代开始之前,将训练数据的顺序随机打乱
# 然后再按每次取batch_size条数据的方式取出
np.random.shuffle(training_data)
# 将训练数据进行拆分,每个mini_batch包含batch_size条的数据
mini_batches = [training_data[k:k + batch_size] for k in range(0, n, batch_size)]
for iter_id, mini_batch in enumerate(mini_batches):
# print(self.w.shape)
# print(self.b)
x = mini_batch[:, :-1]
y = mini_batch[:, -1:]
a = self.forward(x)
loss = self.loss(a, y)
gradient_w, gradient_b = self.gradient(x, y)
self.update(gradient_w, gradient_b, eta)
losses.append(loss)
print('Epoch {:3d} / iter {:3d}, loss = {:.4f}'.
format(epoch_id, iter_id, loss))
return losses
train_data, test_data = load_data()
# 创建网络
net = Network(13)
# 启动训练
losses = net.train(train_data, num_epochs=50, batch_size=100, eta=0.1)
# 画出损失函数的变化趋势
plot_x = np.arange(len(losses))
plot_y = np.array(losses)
plt.plot(plot_x, plot_y)
plt.show()
通过结果我们可以看到loss在逐渐减少。