推荐阅读

前言

​ 多元线性回归总体思路和一元线性回归相同，都是代价函数+梯度下降，所以中心思想参考一元线性回归。

​ 在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。

​ 在开始之前，你可能需要先了解以下知识：

前置知识

向量化

​ 当我们想自己实现向量的点乘时，通常会想到利用for循环来完成，例如有\(\vec{a}\cdot\vec{b}\)，可以写为：

1
2
3
4
x=0
for i in range(a.shape[0]):
    x = x + a[i] * b[i]
return x

​ 然而，Python的Numpy库提供了一个dot()函数可以帮助我们进行向量化计算，作用是加快向量的运算速度，在数据量较大时效率会明显提高。其原理是Numpy利用了计算机底层硬件单指令多数据(SIMD)管道，这在数据集非常大的机器学习中至关重要。所以，向量化是机器学习中一个非常重要和实用的方法。

​ 下图是使用Numpy的dot函数与自己利用for循环实现的向量点乘分别对长度各为10000000的向量a、b点乘运行时间比较：

特征缩放

为什么要特征缩放

​ 当有多个特征时，在有的算法中，我们需要使这些特征具有相似的尺度（无量纲化）。

​ 这里介绍特征缩放在多元线性回归中的作用。

​ 拿下面”问题引入“里得数据来说，各个特征的范围差距太大，我们将每个特征对价格的影响可视化，可以看出哪些因素对价格影响更大。会得到以下图像：

​ 由于各个特征的数量差距过大，代价函数的等高线将会是扁长的，在梯度下降时也会是曲折的，而且计算时长相对会很长（因为学习率是通用的，为了照顾尺度大的特征，学习率必须设置的很小，学习率越小，下降速度就越慢）：

​ 特征缩放将每个特征的范围统一，可以使梯度下降变”平滑“，并且大大提高计算速度（因为可以调大学习率）。

特征缩放的方法

​ 特征缩放的方法有许多种，这里介绍两种：

均值归一化

​ 公式：\(x_i = \frac{x_i - \mu}{max-min}\)

​ 其中，\(\mu\)为样本中该特征的均值

1
2
3
4
# 均值归一化
def MeanNormalization(x):
   '''x为列表'''
    return [(float(i)-np.mean(x))/float(max(x)-min(x)) for i in x]

Z-score标准化(推荐)

​ 公式：\(x_j^{(i)} = \frac{x_j^{(i)}-\mu_j}{\delta_j}\)

​ 其中，\(j\)为\(X\)矩阵中的特征（或列),\((i)\)为样本序号。\(u_j\)为特征\(j\)的均值，\(\delta_j\)为特征\(j\)的标准差。

​ \(\mu_j = \frac{1}{m}\sum_{i=0}^{m-1}x_j^{(i)}\)

​ \(\delta_j^2 = \frac{1}{m}\sum_{i=0}^{m-1}(x_j^{(i)}-\mu_j)^2\)

1
2
3
4
5
6
7
8
9
10
# Z-score标准化
def Zscore(X):
    '''x是(m,n)的矩阵，m为样本个数，n为特征数目'''
    # 找每列(特征)均值
    mu = np.mean(X,axis=0)
    # 找每列(特征)标准差
    sigma = np.std(X,axis=0)
    X_norm = (X - mu) / sigma

    return X_norm

​ 或者使用sklearn:

1
2
3
4
from sklearn import preprocessing

# X为同上矩阵
X_norm = preprocessing.scale(X)

问题引入

​ 示例：现在你有以下数据，请利用这些值构建一个线性回归模型，并预测一栋1200平米，3间卧室，1层，年龄为40年的房屋的价格。

面积(平方)卧室数量楼层数房屋年龄价格（1000s dollars）
9522165271.5
12443264300
19473217509.8
....................

多元线性回归模型

多元线性回归函数

​ 对于上面提到的数据，一共有四种特征（面积，卧室数量，楼层，房屋面积），记为\(x_1,x_2,x_3,x_4\),每个特征分别需要一个\(w\),所以对应的线性回归函数为\(f(x) = w_1x_1+w_2x_2+w_3x_3+w_4x_4+b\).

​ 推广到一般多元线性回归函数，即：

​ \(f(x) = w_1x_1+...+w_nx_n+b\) ,    其中，n为特征数量。

​ 观察\(f(x)\),我们发现可以将\(w\)看作一列，\(x\)看作一列。于是\(f(x)\)又可以写为：\(f_{\vec{w},b}(\vec{x}) = \vec{w}\cdot\vec{x}+b\),   （注意为点乘）

​ 这样我们的目标便是利用已知数据通过梯度下降算法找到最合适的\(\vec{w}\)和b。

转化为矩阵

​ 我们可以将训练集转化为\((m,n)\)的矩阵，m表示示例，n表示特征，于是训练集X可以写为：

\[ \begin{pmatrix} x_{0}^{(0)}&x_{1}^{(0)}&...& x_{n-1}^{(0)} \\ x_{0}^{(1)}&x_{1}^{(1)}&...& x_{n-1}^{(1)} \\...\\x_{0}^{(m-1)}&x_{1}^{(m-1)}&...& x_{n-1}^{(m-1)}\end{pmatrix} \quad \] ​ 注：\(\vec{x}^{(i)}\)表示含有第i个示例的向量，\(x_j^{(i)}\)表示第i个示例的第j个特征。因为每种特征对应一个\(w_i\),所以有向量：

\[ \vec{w} = \begin{pmatrix}w_0\\w_1\\...\\w_{n-1}\end{pmatrix} \quad \]

多元线性回归模型的代价函数

​ 因为\(\hat{y}^{(i)} = f_{\vec{w},b}(\vec{x}^{(i)})=\vec{w}\cdot\vec{x}^{(i)}+b\),    其中\(\hat{y}^{(i)}\)为将第i个示例的向量带入回归函数得出的预测值。那么根据一元线性回归的代价函数的定义，可得多元线性回归模型的代价函数为：

​ \(J(\vec{w},b) = \frac{1}{2m}\sum_{i=0}^{m-1}(\hat{y}^{(i)}-y^{(i)})^2= \frac{1}{2m}\sum_{i=0}^{m-1}(f_{\vec{w},b}(\vec{x}^{(i)})-y^{(i)})^2\)

多元线性回归模型梯度下降函数

​ 重复同时进行下列行为直到收敛：

​ \(w_j = w_j - a\frac{\partial J(\vec w,b)}{\partial w_j},\quad j = 0,1,...,n-1\)

​ \(b = b - a\frac{\partial J(\vec w,b)}{\partial b}\)

​ 其中，a为学习率。化简偏导得：

​ \(\frac{\partial J(\vec w,b)}{\partial w_j} = \frac{1}{m}\sum_{i=0}^{m-1}(f_{\vec w,b}(\vec x^{(i)})-y^{(i)})x_j^{(i)}\)

​ \(\frac{\partial J(\vec w,b)}{\partial b} = \frac{1}{m}\sum_{i=0}^{m-1}(f_{\vec w,b}(\vec x^{(i)})-y^{(i)})\)

​ 其中，n是特征数量，m是训练集示例数量。

​ 收敛时的\(\vec w,b\)即为所求。

问题解析（含代码）

导入并标准化训练集

​ 数据就不上传了，如果想动手测验，直接搞个矩阵把上面那三行当成训练集就行。

注意标准化时需要把均值和标准差也返回过去，便于对测试数据进行缩放。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import copy
import numpy as np

# Z-score标准化
def Zscore(X):
    '''X是(m,n)的矩阵，m为样本个数，n为特征数目'''

    # 找每列(特征)均值
    mu = np.mean(X,axis=0)
    # 找每列(特征)标准差
    sigma = np.std(X,axis=0)
    X_norm = (X - mu) / sigma

    return X_norm,mu,sigma

if __name__ == '__main__':
    data = np.loadtxt('houses.txt', dtype=np.float32, delimiter=',')
    X_train = data[:,:4]
    Y_train = np.ravel(data[:,4:5]) # 把Y从列转为一维数组
    # 将X训练集拿去标准化
    X_train = Zscore(X_train)
    # print(X_train)
    # print(Y_train)

多元线性代价函数

​ 这里只是把代码写出来，实际上计算回归函数时用不上。

​ 对应公式：\(J(\vec{w},b) = \frac{1}{2m}\sum_{i=0}^{m-1}(\hat{y}^{(i)}-y^{(i)})^2= \frac{1}{2m}\sum_{i=0}^{m-1}(f_{\vec{w},b}(\vec{x}^{(i)})-y^{(i)})^2\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 代价函数
def compute_cost(X, y, w, b):
    '''
    :param X: 矩阵，每一行代表一个示例向量
    :param y: 列表
    :param w: 向量
    :param b: 标量
    :return: cost(标量)
    '''
    m = X.shape[0]
    cost = 0.0
    for i in range(m):
        f_x_i = np.dot(X[i],w)
        cost = cost + (f_x_i - y[i])**2
    cost = cost/(2*m)
    return cost

梯度计算函数

​ 对应公式：

​ \(\frac{\partial J(\vec w,b)}{\partial w_j} = \frac{1}{m}\sum_{i=0}^{m-1}(f_{\vec w,b}(\vec x^{(i)})-y^{(i)})x_j^{(i)}\)

​ \(\frac{\partial J(\vec w,b)}{\partial b} = \frac{1}{m}\sum_{i=0}^{m-1}(f_{\vec w,b}(\vec x^{(i)})-y^{(i)})\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# 计算梯度
def compute_gradient(X, y, w, b):
    '''

    :param X: 矩阵，每一行代表一个示例向量
    :param y: 列表
    :param w: 向量
    :param b: 标量
    :return:
        sum_dw(标量，wj对于代价函数的偏导)
        sum_db(标量，b对于代价函数的偏导)
    '''
    m,n = X.shape
    sum_dw = np.zeros((n,))
    sum_db = 0

    for i in range(m):
        err = (np.dot(X[i],w) + b) - y[i]
        for j in range(n):
            sum_dw[j] = sum_dw[j] + err*X[i,j]
        sum_db = sum_db + err
    sum_dw = sum_dw/m
    sum_db = sum_db/m

    return sum_dw,sum_db

梯度迭代函数

​ 对应公式：

     重复同时进行下列行为直到收敛。

     \(w_j = w_j - a\frac{\partial J(\vec w,b)}{\partial w_j},\quad j = 0,1,...,n-1\)

     \(b = b - a\frac{\partial J(\vec w,b)}{\partial b}\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 梯度下降
def gradient_descent(X, y, w_init, b_init,alpha, num_iters):
    '''
    :param X: 矩阵，每一行代表一个示例向量
    :param y: 列表
    :param w_init: w初始值，向量
    :param b_init: b初始值，标量
    :param alpha: 学习率
    :param num_iters: 迭代次数
    :return: 训练出的w和b
    '''
    w = copy.deepcopy(w_init)
    b = b_init

    for i in range(num_iters):
        dw_j, db = compute_gradient(X, y, w, b)
        w = w - alpha * dw_j
        b = b - alpha * db
    return w,b

设置学习率，计算回归函数和预测

​ 预测时必须先使用之前计算出的均值和标准差把测试数据缩放再进行预测！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
if __name__ == '__main__':
    data = np.loadtxt('houses.txt', dtype=np.float32, delimiter=',')
    X_train = data[:,:4]
    Y_train = np.ravel(data[:,4:5]) # 把Y转为一维数组
    # 将X训练集拿去标准化
    X_train,mu,sigma = Zscore(X_train)
    # print(X_train)
    # print(Y_train)
    # 设置学习率等
    w_init = np.array([0,0,0,0])
    b_init = 0
    alpha = 0.5
    num_iters = 10000
    w,b = gradient_descent(X_train,Y_train,w_init,b_init,alpha,num_iters)
    print(f"预测出的w:{w}")
    print(f"预测出的b:{b}")
    # 预测
    X_forecast = [1200,3,1,40]
    X_forecast = (X_forecast-mu)/sigma
    print(f"压缩后的测试数据{X_forecast}")
    X_predice_price = np.dot(X_forecast,w)+b
    print(f"一个1200平米，3个卧室，1层楼，存在40年的房屋 = ${X_predice_price*1000:0.0f}")

​ 结果：

​ 就当那个世界的房子都这么便宜吧。

​ 关于线性回归模型的评估：

​ 常用的有均方误差和均方根误差等方法，这里留个空以后出篇文章，暂时放一放。

补充：学习率a的评估

​ 一般可以通过两个方面来评估：

​ 1.每次计算出的梯度，如果学习率合适，那么这个值应该在不断下降。

​ 2.梯度下降在整个代价函数上的”跳跃“变化，如果学习率合适，它应该是不断往下跳的。

​ 学习率过大：如果学习率过大，每次计算出的梯度会发生类似下面左图的情况。右图是学习率过大的情况下每次梯度下降在总成本函数上的变化，注意它是从下往上跳的，最后会导致结果发散。

​ 学习率过小：它虽然可以收敛，但是过程比较慢。

​ 学习率合适：

补充：浅谈正规方程

​ 求解线性回归不止有梯度下降一种方法，其实还有另一种方法求解，即正规方程。不过在这里不讨论它的详细公式，只拿它与梯度下降做一个对比。

梯度下降正规方程
需要选择学习率不需要选择学习率
需要多次迭代一次计算得出
当特征数量较多时也能较好使用时间复杂度为\(O(n^3)\),当特征数量较多时太慢
也适用于除线性回归的其他模型只适用于线性回归模型

补充：利用sklearn完成线性回归预测

关于StandardScaler()函数：sklearn中StandardScaler()

关于SGDRegressor：随机梯度线性回归，随机梯度下降是不将所有样本都进行计算后再更新参数，而是选取一个样本，计算后就更新参数。

关于LinearRegression：也是线性回归模型，这里没用，可以自己查。

​ 绘图模块的代码可能要自己改改。其余详细说明见代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import numpy as np
np.set_printoptions(precision=2) # 设置numpy可见小数点
import matplotlib.pyplot as plt
dlblue = '#0096ff'; dlorange = '#FF9300'; dldarkred='#C00000'; dlmagenta='#FF40FF'; dlpurple='#7030A0';
plt.style.use('../deeplearning.mplstyle')

from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.preprocessing import StandardScaler

if __name__ == '__main__':
    data = np.loadtxt('houses.txt', dtype=np.float32, delimiter=',')
    X_train = data[:, :4]
    X_features = ['size(sqft)','bedrooms','floors','age']
    y_train = np.ravel(data[:, 4:5])  # 把Y转为一维数组

    # 使用sklearn进行Z-score标准化
    scaler = StandardScaler()
    X_norm = scaler.fit_transform(X_train) # 标准化训练集X
    # print(X_norm)

    # 创建回归拟合模型
    sgdr = SGDRegressor(max_iter=1000) # 设置最大迭代次数
    sgdr.fit(X_norm, y_train)
    print(sgdr)
    print(f"完成的迭代次数: {sgdr.n_iter_}, 权重更新数: {sgdr.t_}")

    # 查看参数
    b_norm = sgdr.intercept_
    w_norm = sgdr.coef_
    print(f"SGDRegressor拟合参数:       w: {w_norm}, b:{b_norm}")
    print(f"之前自己编写的线性回归拟合参数: w: [110.61 -21.47 -32.66 -37.78], b: 362.23")

    # 使用sgdr.predict()进行预测
    y_pred_sgd = sgdr.predict(X_norm)
    # 使用w和b进行预测
    y_pred = np.dot(X_norm, w_norm) + b_norm
    print(f"使用 np.dot() 预测值是否与sgdr.predict预测值相同: {(y_pred == y_pred_sgd).all()}")
    print(f"训练集预测:\n{y_pred[:4]}")
    print(f"训练集真实值\n{y_train[:4]}")

    # 绘制训练集与预测值匹配情况
    fig, ax = plt.subplots(1, 4, figsize=(12, 3), sharey=True)
    for i in range(len(ax)):
        ax[i].scatter(X_train[:, i], y_train, label='target')
        ax[i].set_xlabel(X_features[i])
        ax[i].scatter(X_train[:, i], y_pred, color='orange', label='predict')
    ax[0].set_ylabel("Price")
    ax[0].legend()
    fig.suptitle("target versus prediction using z-score normalized model")
    plt.show()

​ 结果：

​ 训练集预测值与真实值结果对比