多元线性回归

【假设形式】

在回归分析中，如果有两个及以上的自变量，且因变量和自变量之间是线性关系，则称为多元线性回归（Multivariate Linear Regression）

其假设形式如下：

对于给定的容量为 $n$ 的训练集 $D=\{(\mathbf{x_1},y_1),(\mathbf{x_2},y_2),…,(\mathbf{x_n},y_n)\}$，第 $i$ 组样本中的输入 $\mathbf{x_i}$ 具有 $m$ 个特征值，即：$\mathbf{x_i}=(x_i^{(1)},x_i^{(2)},…,x_i^{(m)})\in \mathbb{R}^m$，输出为 $y_i$，多元线性回归学习到的模型为 $f(\mathbf{x_i};\boldsymbol{\theta})$，使得 $f(\mathbf{x_i};\boldsymbol{\theta})\simeq y_i$

用假设函数 $f(\mathbf{x_i};\boldsymbol{\theta})$ 来表示对第 $i$ 组数据的预测结果：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

其中，特征参数 $\boldsymbol{\theta}$ 为 $(m+1)\times 1$ 的列向量，即：

$$\boldsymbol{\theta}=[\theta^{(0)},\theta^{(1)},...,\theta^{(m)}]^T\in \mathbb{R}^{m+1}$$

为了表述方便，对假设函数进行简化，定义一个额外的第 $0$ 个特征量，这个特征量对所有样本的取值全部为 $1$，这使得特征量从过去的 $m$ 个变为 $m+1$ 个，即设：$x_i^{(0)}=1$

那么假设函数就可以写为：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} x_i^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

与一元线性回归模型相同，需要一个损失函数损失函数 $J(\boldsymbol{\theta})$ 来作为衡量预测结果的指标，对于损失函数最小化时的 $\boldsymbol{\theta}$，通常使用最小二乘法或以梯度下降法为代表的迭代法来求解

无论使用何种方法，最终的目标，都是要令这个损失函数的值最小化，即

$$\boldsymbol{\theta}^* = \arg \min \limits_{\boldsymbol{\theta}} \:J(\boldsymbol{\theta})$$

【梯度下降法求解】

对于给定的容量为 $n$ 的样本集 $T=\{(\mathbf{x_1},y_1),(\mathbf{x_2},y_2),…,(\mathbf{x_n},y_n)\}$，第 $i$ 组样本中的输入 $\mathbf{x_i}$ 具有 $m$ 个特征值，即：$\mathbf{x_i}=(x_i^{(1)},x_i^{(2)},…,x_i^{(m)})\in \mathbb{R}^m$

设习得的模型的假设函数为：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

其中，特征参数 $\boldsymbol{\theta}$ 为 $(m+1)\times 1$ 的列向量，即：

$$\boldsymbol{\theta}=[\theta^{(0)},\theta^{(1)},...,\theta^{(m)}]^T\in \mathbb{R}^{m+1}$$

为了表述方便，对假设函数进行简化，定义一个额外的第 $0$ 个特征量，这个特征量对所有样本的取值全部为 $1$，这使得特征量从过去的 $m$ 个变为 $m+1$ 个，即设：$x_i^{(0)}=1$

那么假设函数就可以写为：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} x_i^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

设损失函数为：

$$J(\boldsymbol{\theta})=\frac{1}{2}\big[f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big]^2$$

之所以要乘以 $\frac{1}{2}$，是因为在求导后会带来 $\times 2$，不利于表达与计算，当乘以 $\frac{1}{2}$ 后，求导带来的 $\times 2$ 就与 $\frac{1}{2}$ 抵消，从而简化计算

目标是通过最小化代价函数 $J(\boldsymbol{\theta})$ 来在参数空间 $\Theta$ 中找到合适的 $\boldsymbol{\theta}$ 参数，即：

$$\boldsymbol{\theta}^*=\arg \min_{\Theta}\:J(\boldsymbol{\theta})$$

在最小化代价函数 $J(\boldsymbol{\theta})$ 时，其核心是损失函数对应的梯度函数，即将下列公式重复直到收敛为止：

$$\boldsymbol{\theta_{k+1}} = \boldsymbol{\theta_k} - \alpha \frac{\partial}{\partial \boldsymbol{\theta_k}}J(\boldsymbol{\theta_k})$$

根据链式法则：若 $z=f[g(x)]$，则 $z’=f[g(x)]’+g(x)’$

则对于偏导数项 $\frac{\partial}{\partial \boldsymbol{\theta}}J( \boldsymbol{\theta})$ 有：

$$\begin{align} \frac{\partial}{\partial \boldsymbol{\theta}}J( \boldsymbol{\theta}) &= \frac{\partial}{\partial \boldsymbol{\theta}} \frac{1}{2}\big[f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big]^2 \notag \\ &= \frac{1}{2} \cdot \frac{\partial}{\partial \boldsymbol{\theta}} \Big[\big(f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big)^2\Big] \cdot \Big[ \frac{\partial}{\partial \boldsymbol{\theta}} \big(f(\mathbf{x_i};\boldsymbol{\theta}) - y_i\big)\Big] \notag \\ &= \big(f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big) \cdot \Big[ \frac{\partial}{\partial \boldsymbol{\theta}} \big(f(\mathbf{x_i};\boldsymbol{\theta}) - y_i\big)\Big] \notag \\ &= \big(f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big) \cdot \frac{\partial}{\partial \boldsymbol{\theta}} \bigg( \sum_{j=0}^m \theta^{(j)} x_i^{(j)} - y_i \bigg) \notag \\ &= \big(f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big) \cdot \bigg[ \frac{\partial}{\partial \boldsymbol{\theta}} \sum_{j=0}^m \theta^{(j)} x_i^{(j)} - 0 \bigg] \notag \\ &= \big(f(\mathbf{x_i};\boldsymbol{\theta})-y_i\big) \cdot \frac{\partial}{\partial \boldsymbol{\theta}} \big(\theta^{(0)}x_i^{(0)}+\theta^{(1)}x_i^{(1)}+...+\theta^{(m)}x_i^{(m)}\big) \notag \end{align}$$

对于求和函数展开的部分，当对其中第 $j$ 项求导时，其余各项实际上就是一个常数，它们在求导这一刻是固定不变的，因此有：

$$\frac{\partial}{\partial \boldsymbol{\theta}} f(\mathbf{x_i};\boldsymbol{\theta}) = 0+...+ \frac{\partial}{\partial \boldsymbol{\theta}} \theta^{(j)}x_i^{(j)}+...+0$$

对于得到的乘积函数 $\frac{\partial}{\partial \boldsymbol{\theta}} \theta^{(j)}x_i^{(j)}$，有：

$$\frac{\partial}{\partial \boldsymbol{\theta}} \theta^{(j)}x_i^{(j)}= x_i^{(j)}$$

将得到的结果组合回去，即有：

$$\frac{\partial}{\partial \boldsymbol{\theta}} J(\boldsymbol{\theta})= \big[f(\mathbf{x_i};\boldsymbol{\theta}) -y_i\big] \cdot x_i^{(j)}$$

故而，只需将下列公式重复至收敛即可：

$$\boldsymbol{\theta_{k+1}} = \boldsymbol{\theta_k} - \alpha \big[f(\mathbf{x_i};\boldsymbol{\theta}) -y_i \big] \cdot x_i^{(j)},j = 0,1,...,m$$

其中，$x_i^{(0)}=1$

关于梯度下降法的具体介绍，详见：梯度下降法

【最小二乘法求解】

将数据集 $D$ 写为 $(m+1)\times n$ 的矩阵，即：

$$X=\begin{bmatrix} x_{1}^{(0)} & x_{2}^{(0)} & ... & x_{n}^{(0)} \\ x_{1}^{(1)} & x_{2}^{(1)} & ... & x_{n}^{(1)} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1}^{(m)} & x_{2}^{(m)} & ... & x_{n}^{(m)} \end{bmatrix}$$

同时，将样本中的 $y_i$ 也写为矩阵形式，即输出变量 $Y$ 为 $n\times 1$ 的列向量：

$$Y=[y_1,y_2,...,y_n]^T \in \mathbb{R}^{n}$$

选用残差平方和 RSS 作为损失函数，则有：

$$\begin{align} J(\boldsymbol{\theta}) &= \sum_{i=1}^n (f(\mathbf{x_i})-y_i)^2 \notag \\ &= \sum_{i=1}^n(\sum_{j=0}^{m+1}\theta^{(j)}x_{i}^{(j)}-y_i)^2 \notag \\ &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 \notag \\ &= (X^T\boldsymbol{\theta}-Y)^T(X^T\boldsymbol{\theta}-Y) \notag \\ &= (\boldsymbol{\theta}^TX -Y^T)(X^T\boldsymbol{\theta}-Y) \notag \\ &= \boldsymbol{\theta}^TXX^T\boldsymbol{\theta} -\boldsymbol{\theta}^TXY -Y^TX^T\boldsymbol{\theta} +Y^TY \notag \end{align}$$

要令目标函数最小，显然要令 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})=0$

首先求 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})$，有：

$$\frac{\partial}{\partial \boldsymbol{\theta}}J(\boldsymbol{\theta})=2XX^T \boldsymbol{\theta} -2XY$$

解得：

$$\boldsymbol{\theta}=(XX^T)^{-1}XY$$

其中，$XX^T$ 为满秩矩阵，$(XX^T)^{-1}$ 为对应的逆矩阵

因此，只要根据样本给出的输入 $X$ 与输出 $Y$，若 $(XX^T)^{-1}$ 存在，即可计算出 $\boldsymbol{\theta}$ 的解析解

关于最小二乘法的具体介绍，详见：最小二乘法

【sklearn 实现】

以 sklearn 中的波士顿房价数据集为例，实现多元线性回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score

# 特征提取
def deal_data():
    boston = load_boston()  # sklearn的波士顿房价数据集
    df = pd.DataFrame(boston.data, columns=boston.feature_names)
    df['result'] = boston.target
    data = np.array(df)
    return data[:, :-1], data[:, -1]

# 模型训练
def train_model(features, labels):  
    # 建立线性回归模型
    model = LinearRegression()
    
    # 训练
    model.fit(features, labels)
    return model

# 模型评估
def estimate_model(y_true, y_pred):
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(MSE)
    MAE = mean_absolute_error(y_true, y_pred)
    R2 = r2_score(y_true, y_pred)
    indicators = {"MSE": MSE, "RMSE":RMSE, "MAE":MAE, "R2":R2}
    return indicators

# 可视化
def visualization(y_true, y_pred, model):
    # 绘图
    plt.plot(range(y_true.shape[0]), y_true, "b-") 
    plt.plot(range(y_true.shape[0]), y_pred, "r-.")
    plt.legend(["original value", "predicted value"])
    plt.xlabel("samples", fontsize="15")
    plt.ylabel("y", fontsize="15")
    
    plt.show()

if __name__ == "__main__":
    # 特征提取
    x, y = deal_data()
    
    # 简单交叉验证
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3,random_state=0)
    
    # 模型训练
    model = train_model(x_train, y_train)
    
    # 预测结果
    y_pred = model.predict(x_test) # predict()输入输出均为二维
    print("y test:", y_test[:10]) # 测试集y值
    print("y pred:", y_pred[:10]) # 预测y值
    
    # 模型评估
    indicators = estimate_model(y_test, y_pred)
    print("MSE:", indicators["MSE"])
    print("RMSE:", indicators["RMSE"])
    print("MAE:", indicators["MAE"])
    print("R2:", indicators["R2"])
    
    # 可视化
    visualization(y_test, y_pred, model)