碧海苍梧
2023/03/10阅读：22主题：默认主题
DLT-03-多元线性回归
本文是深度学习入门（deep learning tutorial, DLT）系列的第三篇文章，主要介绍一下多元线性回归。想要学习深度学习或者想要了解机器学习的同学可以关注公众号GeodataAnalysis，我会逐步更新这一系列的文章。
多元线性回归和一元线性回归的模型表达基本一致，区别只在于输入的变量个数更多。因此，二者在假设函数和梯度下降方面有些区别，本文会在上一篇一元线性回归的基础上对其进行详细介绍。为了准确的衡量模型的精度，本文介绍了最小二乘法，使用该方法计算出的假设函数能最小化使损失函数的值。此外，多个变量的量纲不一定一致，因此本文还会介绍一下无量纲化在机器学习中的作用。
1 测试数据集
关注公众号GeodataAnalysis，回复20230302下载测试数据。测试数据为波士顿房价数据集，数据集的因变量为房价；共有十三个自变量，表示影响房价的十三个因素。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv('./housing.data', delim_whitespace=' ')
df.head()
| CRIM | ZN | INDUS | CHAS | NOX | RM | AGE | DIS | RAD | TAX | PTRATIO | B | LSTAT | MEDV | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.00632 | 18.0 | 2.31 | 0.0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1.0 | 296.0 | 15.3 | 396.90 | 4.98 | 24.0 |
| 1 | 0.02731 | 0.0 | 7.07 | 0.0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2.0 | 242.0 | 17.8 | 396.90 | 9.14 | 21.6 |
| 2 | 0.02729 | 0.0 | 7.07 | 0.0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2.0 | 242.0 | 17.8 | 392.83 | 4.03 | 34.7 |
| 3 | 0.03237 | 0.0 | 2.18 | 0.0 | 0.458 | 6.998 | 45.8 | 6.0622 | 3.0 | 222.0 | 18.7 | 394.63 | 2.94 | 33.4 |
| 4 | 0.06905 | 0.0 | 2.18 | 0.0 | 0.458 | 7.147 | 54.2 | 6.0622 | 3.0 | 222.0 | 18.7 | 396.90 | 5.33 | 36.2 |
x = df[df.columns[:-1]].to_numpy().T
y = df[df.columns[-1]].to_numpy()
x.shape, y.shape
((13, 506), (506,))
2 假设函数
多元线性回归的输入变量是多维变量，其表示方法为：
容纳这些多个特征的假设函数的多变量形式如下：
使用矩阵乘法的定义，我们的多变量假设函数可以简洁地表示为如下形式，其中 为全是1的矩阵，用于表示偏置项：
代码表示如下：
def hypothesis_fun(x, parameters):
return parameters @ x
3 损失函数
多元线性回归的损失函数与一元线性回归一致，也是平⽅误差函数，即：
代码表示如下：
def loss_fun(x, y, parameters):
y_predict = hypothesis_fun(x, parameters)
return np.sum(np.square(y_predict-y))/(2*y.size)
4 梯度下降
多元线性回归的梯度下降算法和一元线性回归基本一致，区别只在于参数的数量多了一些，具体如下：
代码表示如下：
def gradient_decent(x, y, parameters, learning_rate):
y_predict = hypothesis_fun(x, parameters)
return parameters - learning_rate*(y_predict-y) @ x.T/y.size
5 最小二乘法
最小二乘法就是要找到一组 使得 （残差平方和) 最小，即，求 。
最小二乘法的代数法解法就是对 求偏导数，令偏导数为0，再解方程组，得到 。矩阵法比代数法要简洁，下面主要讲解下矩阵法解法。
假设函数 的矩阵并表达形式为：
其中， 为 维的矩阵， 代表样本的特征数， 代表样本的个数。假设函数 为 的向量, 为 的向量，里面有 个代数法的模型参数。
损失函数定义为 ，其中 是样本的输出向量，维度为 。 在这主要是为了求导后系数为1，方便计算。
根据最小二乘法的原理，我们要对这个损失函数的 向量求导取0。结果如下式：
对上述求导等式整理后可得：
如果 不可逆，通常有以下两种原因：
-
冗余特征，其中两个特征非常密切相关（即它们是线性依赖的） -
特征太多（例如m≤n）。 在这种情况下，应删除某些特征或使用“正则化”。
使用最小二乘法计算假设函数的参数的代码如下：
x2 = np.vstack((np.ones(x.shape[1]), x))
parameters = np.linalg.inv(x2 @ x2.T) @ x2 @ y
6 无量纲化
本文介绍的无量纲化方法有两种，一是特征缩放，其计算公式为：
二是均值标准化，其计算公式为：
这两个无量纲化算法的代码表示如下：
def feature_scale(x):
x_min = np.min(x, axis=1)[..., np.newaxis]
x_max = np.max(x, axis=1)[..., np.newaxis]
return (x - x_min)/(x_max - x_min)
def mean_normalization(x):
x_min = np.min(x, axis=1)[..., np.newaxis]
x_max = np.max(x, axis=1)[..., np.newaxis]
x_mean = np.mean(x, axis=1)[..., np.newaxis]
return (x - x_mean)/(x_max - x_min)
7 模型训练及预测
多元线性回归的训练步骤和一元线性回归基本一致，唯一的区别在于在正式训练之前对输入变量进行了无量纲化。此外，在输入变量x矩阵的第一行前插入全是1的一行，用于表示偏置项。
parameters = np.random.rand(x.shape[0]+1)
learning_rate = 0.1
losses = []
x2 = feature_scale(x)
x2 = np.vstack((np.ones(x2.shape[1]), x2))
batch_size = 300
epoch_size = 100
for epoch in range(epoch_size):
for i in range(y.size//batch_size+1):
random_samples = np.random.choice(x2.shape[1], batch_size)
parameters = gradient_decent(x2[:, random_samples],
y[random_samples],
parameters,
learning_rate)
losses.append(loss_fun(x2[:, random_samples],
y[random_samples],
parameters))
plt.plot(losses);
现在我们的模型已经训练好了，下面我把模型预测的代码封装一下，具体如下：
def predict(x, parameters):
x = feature_scale(x)
x = np.vstack((np.ones(x.shape[1]), x))
return parameters @ x
可视化预测结果的代码如下，这里我们对比了最小二乘法的结算结果和我们模型的结果，可以发现二者相差不大，具体如下：
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 9))
ax1.plot(losses)
ax1.set_title('Loss', fontdict={'size': 30})
ax2.scatter(predict(x, parameters), y)
ax2.plot(y, y, 'r-')
ax2.set_title('Linear Regression', fontdict={'size': 30})
# Best theta
x2 = np.vstack((np.ones(x.shape[1]), x))
theta = np.linalg.inv(x2 @ x2.T) @ x2 @ y
ax3.scatter(theta @ x2, y);
ax3.plot(y, y, 'r-')
ax3.set_title('Least squares', fontdict={'size': 30});
8 对模型进行封装
这一节我们将上文构建的模型封装成了Python的一个类，方便对其调用进行训练、预测等任务。目前这个类较为简单，可以将其视为一个建议的线性神经网络模型，理解这个类的构造对理解后文的代码十分有帮助，建议好好看一下。
class MutiLinearRe():
def __init__(self, input_shape):
self.input_shape = input_shape
self.parameters = np.random.rand(self.input_shape[0]+1)
def hypothesis_fun(self, x, parameters):
return parameters @ x
def loss_fun(self, x, y, parameters):
y_predict = self.hypothesis_fun(x, parameters)
return np.sum(np.square(y_predict-y))/(2*y.size)
def gradient_decent(self, x, y, parameters, learning_rate):
y_predict = self.hypothesis_fun(x, parameters)
return parameters - learning_rate*(y_predict-y) @ x.T/y.size
def feature_scale(self, x):
x_min = np.min(x, axis=1)[..., np.newaxis]
x_max = np.max(x, axis=1)[..., np.newaxis]
return (x - x_min)/(x_max - x_min)
def mean_normalization(self, x):
x_min = np.min(x, axis=1)[..., np.newaxis]
x_max = np.max(x, axis=1)[..., np.newaxis]
x_mean = np.mean(x, axis=1)[..., np.newaxis]
return (x - x_mean)/(x_max - x_min)
def _normalization(self, x, normalization):
if normalization:
if 'feature' == normalization:
x = self.feature_scale(x)
elif 'mean' == normalization:
x = self.mean_normalization(x)
return x
def fit(self, x, y, epoch_size, batch_size, learning_rate, normalization="feature"):
self.learning_rate = learning_rate
self.loss = []
self.normalization = normalization
self.x = self._normalization(x, self.normalization)
self.x = np.vstack((np.ones(self.input_shape[1]), self.x))
self.y = y.copy()
for epoch in range(epoch_size):
for i in range(y.size//batch_size+1):
random_samples = np.random.choice(x2.shape[1], batch_size)
self.parameters = self.gradient_decent(self.x[:, random_samples],
self.y[random_samples],
self.parameters,
self.learning_rate)
loss_ = self.loss_fun(self.x[:, random_samples],
self.y[random_samples],
self.parameters)
if np.isinf(loss_):
raise ValueError("Overflow, please change the learning_rate")
self.loss.append(loss_)
def predict(self, x):
x = self._normalization(x, self.normalization)
x = np.vstack((np.ones(x.shape[1]), x))
return self.parameters @ x
model = MutiLinearRe(input_shape=x.shape)
model.fit(x, y, epoch_size=100, batch_size=300,
learning_rate=0.1, normalization='feature')
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 9))
ax1.plot(model.loss[700:])
ax1.set_title('Loss', fontdict={'size': 30})
ax2.scatter(model.predict(x), y)
ax2.plot(y, y, 'r-')
ax2.set_title('Linear Regression', fontdict={'size': 30})
# Least squares: Best theta
x2 = np.vstack((np.ones(x.shape[1]), x))
theta = np.linalg.inv(x2 @ x2.T) @ x2 @ y
ax3.scatter(theta @ x2, y);
ax3.plot(y, y, 'r-')
ax3.set_title('Least squares', fontdict={'size': 30});
作者介绍