梯度下降法求解多元线性回归--Numpy实现

梯度下降法求解多元线性回归

我们来用编程实现一个多元线性回归问题。
这是样本数据

其中的属性包括房屋面积和房间数。可以看到面积和房间数的取值范围相差很大，如果直接使用这样的数据来训练模型，这个面积的贡献就会远远大于房间数的影响，在学习过程中占主导甚至是决定性的地位，这显然是不合理的，这时候应该首先对属性的值进行归一化处理。

归一化/标准化

归一化又称为标准化，就是将数据的大小限制在一定的范围之内。

在机器学习中对所有属性进行归一化处理，就是让他们处于同一个范围，同一个数量级下，这样才能更加的具有可比性。

使用归一化，不仅能够更快的收敛到最优解，而且还可能提高学习器的精度。

常用的归一化方法有

• 线性归一化
• 标准差归一化
• 非线性映射归一化

线性归一化：对原始数据的线性变换

这是它的转换函数

$x^{\ast }=\frac{x-min}{max-min}$

其中me和max分别是所有数据中的最小值和最大值。

• 线性归一化实现对原始数据的等比例缩放。
• 归一化之后，所有的数据都被映射到0和1之间

这种归一化方法适用于样本分布比较均匀，比较集中的情况。如果最大值或者最小值不稳定和绝大多数数据的差距非常大的时候，使用这种方法得到的结果也会不稳定。

为了避免这种情况，在实际应用中可以使用经验常量来代替最大值和最小值。

标准差归一化：将数据归一化为均值为0，方差为1的标准正态分布

这是它的转换函数

$x^8=\frac{x-\mu }{\sigma }$

其中的$\mu$和$\sigma$分别是均值和标准差。

标准差归一化适标准差规划适用于样本近似于正态分布或者最大最小值未知的情况。有时当最大或者最小值处于孤立点时，也可以使用标准差归一化。

非线性映射归一化：对原始数据的非线性变换

除此之外，还可以使用其他一些非线性的数学函数对原始数据进行映射，常用的函数有对数、指数、正切等。

非线性映射的方法，通常用于数据分化比较大的情况，也就是有的数据很大，有的数据很小，通过这种非线性映射的方法可以使数据变得尽量均匀或者有特点。

样本属性的归一化，需要根据实际数据的分布情况和特点进行定制化的设计，决定采用哪种函数。

线性归一化解决问题

在这节课的例子中，我们采用线性归一化的方法

import numpy as np

area = np.array([137.97, 104.50, 100.00, 124.32, 79.20, 99.00, 124.00, 114.00,
                 106.69, 138.05, 53.75, 46.91, 68.00, 63.02, 81.26, 86.21])
room = np.array([3, 2, 2, 3, 1, 2, 3, 2, 2, 3, 1, 1, 1, 1, 2, 2])

# 这是线性归一化
x1 = (area-area.min())/(area.max()-area.min())
x2 = (room-room.min())/(room.max()-room.min())

print(x1)
print(x2)

数组对象的min和max方法分别返回数组的最小值和最大值。
这是归一化之后的商品房面积可以看到最小的数被映射为了0，最大的数被映射为1，所有的数据都在0和1之间。

[0.99912223 0.63188501 0.58251042 0.84935264 0.3542901  0.57153829
 0.84584156 0.73612025 0.65591398 1.         0.07504937 0.
 0.23140224 0.17676103 0.37689269 0.43120474]

这是归一化之后的房间数线性归一化的结果，其实就是这个数在整个样本中所处的位置的比例。

[1.  0.5 0.5 1.  0.  0.5 1.  0.5 0.5 1.  0.  0.  0.  0.  0.5 0.5]

使用梯度下降法求解多元线性回归

下面我们看一下使用梯度下降法，实现多元线性回归的程序的主要步骤。

• 第1步加载样本数据，area，room，price
• 第2步数据处理，首先对属性归一化，然后生成算法需要的多维数组
• 第3步设置超参数，学习率，迭代次数
• 第4步设置模型参数的初值，这里的$W_0(w_0,w_1,w_2)$是一个向量，其中包括三个模型参数
• 第4步使用权值更新算法训练模型W，
  • 首先计算梯度，
    $\frac{\partial Loss}{\partial W}=X^T(XW-Y)$
  • 然后把它带入迭代公式，就可以实现对权值的更新
    $W^{(k+1)}=W^{(k)}-\eta X^T(XW-Y)$
    $W^{(k+1)}=W^{(k)}-\eta \frac{\partial Loss}{\partial W}$
• 最后以可视化的形式输出结果

加载样本数据

下面我们来看一下程序的实现，首先导入需要的库，设置中文字体，加载数据，这分别是商品房面积，房间数和房价，num是样本数量。

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']

area = np.array([137.97, 104.50, 100.00, 124.32, 79.20, 99.00, 124.00, 114.00,
                 106.69, 138.05, 53.75, 46.91, 68.00, 63.02, 81.26, 86.21])
room = np.array([3, 2, 2, 3, 1, 2, 3, 2, 2, 3, 1, 1, 1, 1, 2, 2])
price = np.array([145.00, 110.00, 93.00, 116.00, 65.32, 104.00, 118.00, 91.00,
                  62.00, 133.00, 51.00, 45.00, 78.50, 69.65, 75.69, 95.30])

num = len(area)

数据处理

创建元素值全部为1的一维数组$x_0$， 对商品房面积和房间数进行归一化处理，将$x_0，x_1，x_2$堆叠为形状为(16, 3)，将房价转化为形状为(16, 1)的二维数组。

x0 = np.ones(num)

x1 = (area - area.min()) / (area.max() - area.min())
x2 = (room - room.min()) / (room.max() - room.min())

X = np.stack((x0, x1, x2), axis=1)
Y = price.reshape(-1, 1)

print(X.shape, Y.shape)
# ((16, 3), (16, 1))

设置超参数

设置学习率，迭代次数和显示间隔

learn_rate = 0.001
iter = 500

display_step = 50

设置模型参数初始值

这里采用随机函数，生成一个形状为(3, 1)的二维数组，作为向量W的初值。

np.random.seed(612)
W = np.random.randn(3, 1)

训练模型

这是计算偏导数

这是更新权值

这是计算房价的估计值这里的X形状为(16, 3)， W形状为(3, 1)相乘的结果，使形状为(16, 1)的二维数组

下面计算军方误差作为损失值，结果是一个数字。

把损失值追加到列表，mse中，每隔50次迭代，显示一下当前的循环次数和损失值

mse = []
for i in range(0, iter + 1):
    dl_dW = np.matmul(np.transpose(X), np.matmul(X, W) - Y)
    W = W - learn_rate * dl_dW

    PRED = np.matmul(X, W)
    Loss = np.mean(np.square(Y - PRED)) / 2

    mse.append(Loss)

    if i % display_step == 0:
        print("i: %i, Loss: %f" % (i, mse[i]))

这是运行的结果，可以看到随着迭代次数的增加，损失逐渐收敛。

结果可视化

最后将训练的过程和结果以可视化的形式输出。

首先设置画布尺寸，然后划分子图，在子图一中绘制损失的变化曲线，在子图二中绘制样本数据和模型的预测数据。

plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.plot(mse)
plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss", fontsize=14)

plt.subplot(1, 2, 2)
PRED = PRED.reshape(-1)
plt.plot(price, color="red", marker="o", label="销售记录")
plt.plot(PRED, color="blue", marker=".", label="预测房价")
plt.xlabel("Sample", fontsize=14)
plt.ylabel("Price", fontsize=14)

plt.legend()
plt.show()

这是运行的结果

现在这个模型就训练好了，可以使用它来预测未知的情况，大家也可以试着修改学习率迭代次数和模型参数的初值，看一下结果有什么不同好。