03. 线性回归
一、线性回归
线性回归（Linear Regression）是一种用于建模两个或多个变量之间线性关系的统计方法。它通过拟合一条直线（或超平面）来描述自变量（输入特征）与因变量（输出目标）之前的关联，并可用于预测或分析变量间的影响关系。假设因变量 y 与自变量 \(x_{1}, x_{2}, ..., x_{n}\) 之间的关系可以用如下线性方程表示：
\[y = \beta_{0} + \beta_{1}x_{1} + ... + \beta_{n}x_{n}
\]
其中，\(\beta_{0}\) 代表 截距，模型在自变量全为 0 时的基准值。\(\beta_{1}, \beta_{2}, ... ,\beta_{n}\) 为 自变量的系数，表示每个自变量对因变量的影响程度。
我们可以在终端中使用 pip 安装 sklearn 机器学习库和 matplotlib 绘图库。默认是从国外的主站上下载，因此，我们可能会遇到网络不好的情况导致下载失败。我们可以在 pip 指令后通过 -i 指定国内镜像源下载。
pip install scikit-learn pandas matplotlib -i https://mirrors.aliyun.com/pypi/simple
国内常用的 pip 下载源列表：
- 阿里云 https://mirrors.aliyun.com/pypi/simple
- 中国科技大学 https://pypi.mirrors.ustc.edu.cn/simple
- 清华大学 https://pypi.tuna.tsinghua.edu.cn/simple
- 中国科学技术大学 http://pypi.mirrors.ustc.edu.cn/simple
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
X = [[5], [8], [10], [12], [15], [3], [7], [9], [14], [6]] # 自变量，每周学习时长
y = [55, 65, 70, 75, 85, 50, 60, 72, 80, 50] # 因变量，每周学习时长对应的考试成绩
model = LinearRegression() # 线性回归模型
model.fit(X, y) # 训练模型
print(model.coef_) # 输出模型的系数
print(model.intercept_) # 输出模型的截距
# 预测
X_new = pd.DataFrame([[10], [12], [15]])
y_pred = model.predict(X_new)
print(y_pred)
# 可视化
x_line = np.arange(0, 15, 0.1).reshape(-1, 1) # 生成0到15的等差数列，用于绘制回归线
y_line = model.predict(x_line) # 计算回归线的y值
plt.plot(x_line, y_line, color='red', linewidth=2) # 绘制回归线
plt.scatter(X, y, color='blue') # 绘制原始数据点
plt.scatter(X_new, y_pred, color='green') # 绘制预测数据点
plt.xlabel('Hours of Study') # 设置x轴标签
plt.ylabel('Test Score') # 设置y轴标签
plt.title('Linear Regression') # 设置标题
plt.show() # 显示图像
二、损失函数
模型的预测值与真实值之间存在误差。损失函数用于衡量模型与真实值之间的误差，并通过最小化损失函数来优化模型参数。当损失函数最小时，得到的自变量系数就是最优解。
均方误差 是回归任务中最常用的损失函数。均方误差对应了欧氏距离。基于均方误差最小化来进行模型求解的方法称为 最小二乘法。在线性回归中，最小二乘法就是试图找到一条直线（或超平面），使所有样本到直线（或超平面）上的欧式距离之和最小。
\[MSE = \frac{1}{n}\sum_{i=1}^{n}{(f(x_{i}) - y_{i})}^{2}
\]
其中，n 为 样本个数，\(y_{i}\) 为 第 i 个样本的真实值，\(f(x_{i})\) 为 第 i 个样本的预测值。
均方误差具有以下特点：
- 均方误差对大误差比较敏感，因为平方项会放大较大的误差。
- 均方误差是凸函数，存在全局的唯一最小值，平方项又使损失函数处处可导，便于求解最优参数。
- 最小二乘法的解析解可通过矩阵运算直接求出。
- 如均方误差服从正态分布，则均方误差对应极大似然估计，是最优的损失函数。
这里，我们以一元线性回归为例，我们将 \(f(x_{i}) = \beta_{0} + \beta_{1}x_{i}\) 代入到 MSE 中，可以得到：
\[MSE = \frac{1}{n}\sum_{i=1}^{n}(\beta_{0} + \beta_{1}x_{i} - y_{i})^{2}
\]
然后，我们分别对 \(\beta_{0}\) 求偏导，可以得到：
\[\begin{align}
\frac{\partial MSE}{\partial \beta_{0}}
& = \frac{\partial}{\partial \beta_{0}}[\frac{1}{n}\sum_{i=1}^{n}{(\beta_{0} + \beta_{1}x_{i} - y_{i})^{2}}] \\
& = \frac{2}{n} \sum_{i=1}^{n}{(\beta_{0} + \beta_{1}x_{i} - y_{i})}
\end{align}
\]
我们令 \(\beta_{0}\) 的偏导等于 0，可以得出：
\[\begin{align}
& \frac{2}{n} \sum_{i=1}^{n}{(\beta_{0} + \beta_{1}x_{i} - y_{i})} = 0 \\
& \sum_{i=1}^{n}{\beta_{0}} = \sum_{i=1}^{n}{y_{i}} - \sum_{i=1}^{n}{\beta_{1}x_{i}} \\
& n{\beta_{0}} = \sum_{i=1}^{n}{y_{i}} - \sum_{i=1}^{n}{\beta_{1}x_{i}} \\
& \beta_{0} = \bar{y} - \beta_{1}\bar{x}
\end{align}
\]
我们对 \(\beta_{1}\) 求偏导，可以得出：
\[\begin{align}
\frac{\partial MSE}{\partial \beta_{1}}
&= \frac{\partial}{\partial \beta_{1}}[\frac{1}{n}\sum_{i=1}^{n}{(\beta_{0} + \beta_{1}x_{i} - y_{i})^{2}}] \\
&= \frac{2}{n}\sum_{i=1}^{n}{(\beta_{0} + \beta_{1}x_{i} - y_{i}).x_{i}} \\
&= \frac{2}{n}\sum_{i=1}^{n}{(\beta_{0}x_{i} + \beta_{1}x_{i}^{2} - x_{i}y_{i})}
\end{align}
\]
我们将 \(\beta_{0} = \bar{y} - \beta_{1}\bar{x}\) 代入 \(\beta_{1}\) 的偏导，可以得出：
\[\begin{align}
\frac{\partial MSE}{\partial \beta_{1}}
&= \frac{2}{n}\sum_{i=1}^{n}{((\bar{y} - \beta_{1}\bar{x})x_{i} + \beta_{1}x_{i}^{2} - x_{i}y_{i})} \\
&= \frac{2}{n}\sum_{i=1}^{n}{(x_{i}\bar{y} - \beta_{1}\bar{x}x_{i} + \beta_{1}x_{i}^{2} - x_{i}y_{i})} \\
&= \frac{2}{n}\sum_{i=1}^{n}{(\beta_{1}(x_{i}^{2} - \bar{x}x_{i}) + x_{i}(\bar{y}-y_{i}))} \\
\end{align}
\]
我们令 \(\beta_{1}\) 的偏导为 0，可以得出：
\[\begin{align}
& \frac{2}{n}\sum_{i=1}^{n}{(\beta_{1}(x_{i}^{2} - \bar{x}x_{i}) + x_{i}(\bar{y}-y_{i}))} = 0 \\
& \sum_{i=1}^{n}{\beta_{1}(x_i^{2} - \bar{x}x_{i})} = \sum_{i=1}^{n}{x_{i}(y_{i} - \bar{y})} \\
& \beta_{1} \sum_{i=1}^{n}{(x_i^{2} - \bar{x}x_{i})} = \sum_{i=1}^{n}{x_{i}(y_{i} - \bar{y})} \\
\end{align}
\]
我们解得 \(\beta_{1}\) 如下：
\[\begin{align}
\beta_{1}
& = \frac{\sum_{i=1}^{n}{x_{i}(y_{i} - \bar{y})}}{\sum_{i=1}^{n}{(x_i^{2} - \bar{x}x_{i})}} \\
&= \frac{\sum_{i=1}^{n}{x_{i}(y_{i} - \bar{y})}}{\sum_{i=1}^{n}{x_{i}((x_i - \bar{x}))}} \\
&= \frac{\sum_{i=1}^{n}{(x_{i} - \bar{x} + \bar{x})(y_{i} - \bar{y})}}{\sum_{i=1}^{n}{(x_{i} - \bar{x} + \bar{x})((x_i - \bar{x}))}} \\
&= \frac{\sum_{i=1}^{n}{(x_{i} - \bar{x})(y_{i} - \bar{y}) + \bar{x}(y_{i} - \bar{y})}}{\sum_{i=1}^{n}{(x_{i} - \bar{x})(x_i - \bar{x}) + \bar{x}(x_i - \bar{x})}} \\
&= \frac{\sum_{i=1}^{n}{(x_{i} - \bar{x})(y_{i} - \bar{y}) + \sum_{i=1}^{n}{(\bar{x}(y_{i} - \bar{y}))}}}{\sum_{i=1}^{n}{(x_{i} - \bar{x})(x_i - \bar{x}) + \sum_{i=1}^{n}{(\bar{x}(x_i - \bar{x}))}}} \\
&= \frac{\sum_{i=1}^{n}{(x_{i} - \bar{x})(y_{i} - \bar{y}) + \sum_{i=1}^{n}{(\bar{x}(y_{i} - \bar{y}))}}}{\sum_{i=1}^{n}{(x_{i} - \bar{x})(x_i - \bar{x}) + \sum_{i=1}^{n}{(\bar{x}(x_i - \bar{x}))}}} \\
\end{align}
\]
其中，
\[\begin{align}
& \sum_{i=1}^{n}{(y_{i} - \bar{y})} = \sum_{i=1}^{n}{y_{i}} - n\bar{y} = 0 \\
& \sum_{i=1}^{n}{(x_{i} - \bar{x})} = \sum_{i=1}^{n}{x_{i}} - n\bar{x} = 0
\end{align}
\]
所以，
\[\beta_{1} = \frac{\sum_{i=1}^{n}{(x_{i} - \bar{x})(y_{i} - \bar{y})}}{\sum_{i=1}^{n}{(x_{i} - \bar{x})^{2}}}
\]
三、正规方程法
如果我们使用的线性回归模型有多元未知数时，我们可以正规方程法求解线性回归的解析解。它是基于最小二乘法，通过求解矩阵方程来直接获取参数值。
将损失函数转换为矩阵形式如下：
\[\begin{align}
MSE &= \frac{1}{n}\sum_{i=1}^{n}{(f(x_{i}) - y_{i})^2} \\
& = \frac{1}{n}\sum_{i=1}^{n}{((\beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + ... + \beta_{m}x_{im}) - y_{i})^{2}} \\
& = \frac{1}{n}\sum_{i=1}^{n}{(\vec{\beta}.\vec{x_{i}} - y_{i})} ^ 2 \\
& = \frac{1}{n}\sum_{i=1}^{n}{(\vec{\beta}^{T} \vec{x_{i}} - y_{i})} ^ 2
\end{align}
\]
其中，向量的 L2 范围如下：
\[\lvert\lvert {X} \rvert\rvert_{2} = \sqrt{\sum_{i=1}^{n}x_{i}^{2}}
\]
其中，\(\vec{\beta_{i}}^{T}\vec{x_{i}}\) 整体可表示为 \(X\beta\)，我们根据向量的 L2 范式展开损失函数可得：
\[\begin{align}
MSE & = \frac{1}{n}\sum_{i=1}^{n}{(\vec{\beta}^{T}\vec{x_{i}} - y_{i})} ^ 2 \\
& = \frac{1}{n}\lvert\lvert {\vec{\beta}^{T}\vec{x} - \vec{y}} \rvert\vert_{2}^{2} \\
& = \frac{1}{n}\lvert\lvert {X\beta - \vec{y}} \rvert\vert_{2}^{2} \\
& = \frac{1}{n}(X\beta - \vec{y})^{T}(X\beta - \vec{y})
\end{align}
\]
其中，X 是 n * (m + 1) 的矩阵，包含一个全 1 的列。
\[X =
\begin{matrix}
1 & x_{11} & x_{12} & ... & x_{1m} \\
1 & x_{21} & x_{22} & ... & x_{2m} \\
... & ... & ...& ...& ...\\
1 & x_{31} & x_{32} & ... & x_{3m}
\end{matrix}
\]
\(\vec{\beta}\) 为 (m + 1) * 1 的参数矩阵，包含截距项 \(\beta{0}\)。
\[\vec{\beta} =
\begin{matrix}
\beta_{0} \\
\beta_{1} \\
... \\
\beta_{m}
\end{matrix}
\]
\(\vec{y}\) 为 n * 1 的因变量矩阵。
我们对 \(\beta\) 求偏导，可以得出：
\[\begin{align}
\frac{\partial MSE}{\partial \beta}
& = \frac{\partial(\frac{1}{n}(\vec{x}\vec{\beta} - \vec{y})^{T}.(\vec{x}\vec{\beta} - \vec{y}))}{\partial \beta} \\
& = \frac{2}{n}X^{T}(X \vec{\beta}- \vec{y})
\end{align}
\]
当 \(X^{T}X\) 为满秩矩阵或正规矩阵时，令偏导等于 0。
\[\frac{2}{n}X^{T}X\vec{\beta} - \frac{2}{n}X^{T}\vec{y} = 0 \\
X^{T}X\vec{\beta} = X^{T}\vec{y} \\
\vec{\beta} = (X^{T}X)^{-1}X^{T}Y
\]
四、梯度下降法
默认情况下，线性回归模型中使用正规方程法求解，它适用于特征数量较少的情况。当特征数量较大时，计算逆矩阵的复杂度会显著增加，此时，更加适合使用梯度下降法。
梯度下降法是一种用于最小化目标函数的迭代优化算法，核心是沿着目标函数（如损失函数）的负梯度方向逐步调整参数，从而逼近函数的最小值。梯度方向指示了函数增长最快的方向，因此负梯度方向是函数下降最快的方向。
假设目标函数为 \(J(\beta)\)，其中 \(\beta\) 是 模型参数，梯度下降的更新公式如下：
\[\beta_{t+1} = \beta_{t} - \alpha.\nabla J(\beta_{t})
\]
其中，\(\alpha\) 是 学习率，用来控制步长，\(\nabla(\beta_{t})\) 是 目标函数在 \(\beta_{t}\) 处的梯度（偏导组成的向量）。
我们求 \(\vec{\beta} = \begin{bmatrix} \beta_{0} & \beta_{1} \end{bmatrix}\) 使得损失函数 \(J(\vec{\beta}) = \frac{1}{n}\sum_{i=1}^{n}{(f(x_{i}) - y_{i})^{2}}\) 最小，为此需要求出损失函数相对 \(\vec{\beta}\) 的梯度。我们为 \(\vec{x}\) 左侧添加一列 1，则求得损失函数的梯度如下：
\[\nabla J(\vec{\beta}) = \frac{2}{n}X^{T}(X\vec{\beta} - \vec{y})
\]
我们使用梯度下降法求解 \(\vec{\beta} = \begin{bmatrix} \beta_{0} & \beta_{1} \end{bmatrix}\) ，其中 \(\vec{\beta}\) 初始值取 \(\begin{bmatrix} 1 & 1 \end{bmatrix}^{T}\)，学习率 \(\alpha\) 取 0.01。
import numpy as np
import matplotlib.pyplot as plt
# 损失函数
def J(beta):
return np.sum((X @ beta - y) ** 2) / n
# 梯度函数
def gradient(beta):
return X.T @ (X @ beta - y) / n * 2
X = np.array([[5], [8], [10], [12], [15], [3], [7], [9], [14], [6]]) # 自变量，每周学习时长
y = np.array([[55], [65], [70], [75], [85], [50], [60], [72], [80], [50]]) # 因变量，每周学习时长对应的考试成绩
n = X.shape[0] # 样本数量，矩阵最外围的元素数量
X = np.hstack((np.ones((n, 1)), X)) # 在X矩阵左侧添加一列1，用于计算截距项
# 超参数
alpha = 0.01 # 学习率
iterations = 10000 # 迭代次数
beta = np.array([[1], [1]]) # 初始参数值
iter = 0
# 定义列表，保存参数变化轨迹
beta0 = []
bete1 = []
# 重复迭代
for i in range(iterations):
beta0.append(beta[0, 0])
bete1.append(beta[1, 0])
grad = gradient(beta) # 计算梯度
beta = beta - alpha * grad # 更新参数值
if (np.abs(grad) < 1e-10).all():
iter = i + 1
break
print(f"Final Iteration: {iter}")
plt.plot(beta0, bete1)
plt.show()
上述的代码，我们手动写的损失函数和梯度函数，在实际开发过程中，我们可以使用机器学习库提供的 SGDRegressor 类直接使用。
SGDRegressor(
loss="squared_error", # 损失函数，均方误差
penalty="l2", # 正则化项，L2正则化
alpha=0.0001, # 正则化系数，控制正则化项的强度
l1_ratio=0.15, # L1正则化系数，控制L1正则化项的强度
fit_intercept=True, # 是否拟合截距项
max_iter=1000, # 最大迭代次数
tol=1e-3, # 停止条件，当损失函数变化小于tol时，停止迭代
shuffle=True, # 是否在每次迭代时打乱数据
verbose=0, # Verbosity mode.
epsilon=DEFAULT_EPSILON, # 用于计算学习率的参数
random_state=None, # 随机种子，用于复现结果
learning_rate="invscaling", # 学习率调度策略，逆缩放学习率
eta0=0.01, # 初始学习率
power_t=0.25, # 学习率调度参数，用于逆缩放学习率
early_stopping=False, # 是否启用早停机制
validation_fraction=0.1, # 验证集比例，用于早停机制
n_iter_no_change=5, # 无改进迭代次数，用于早停机制
warm_start=False, # 是否启用热启动，即从上次训练中恢复
average=False, # 是否启用平均模型，即计算每个迭代的模型参数平均值
)
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import SGDRegressor
X = np.array([[5], [8], [10], [12], [15], [3], [7], [9], [14], [6]]) # 自变量，每周学习时长
y = np.array([55, 65, 70, 75, 85, 50, 60, 72, 80, 50]) # 因变量，每周学习时长对应的考试成绩
# 定义模型
model = SGDRegressor(
loss='squared_error', # 损失函数，均方误差
max_iter=10 ** 6, # 最大迭代次数
eta0=1e-5, # 学习率
learning_rate='constant', # 学习率调度策略，constant，即学习率保持不变
tol=1e-8, # 停止条件，当损失函数变化小于tol时，停止迭代
)
model.fit(X, y) # 训练模型
print(model.coef_) # 输出模型的系数
print(model.intercept_) # 输出模型的截距
学习率过大可能导致跳过最优解，甚至发散，学习率过小可能导致收敛速度满，容易陷入局部最小。
五、广告投放效果预测
Advertising 数据集：https://www.kaggle.com/datasets/tawfikelmetwally/advertising-dataset。
| 字段名 | 值 |
|---|---|
| ID | 序号 |
| TV | 电视广告投放金额，单位千元 |
| Radio | 广播广告投放金额，单位千元 |
| Newspaper | 报纸广告投放金额，单位千元 |
| Sales | 销售额，单位百万元 |
import pandas as pd
import sklearn.model_selection
import sklearn.metrics
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, SGDRegressor
# 读取数据集
dataset = pd.read_csv("assets/data/advertising.csv") # 读取数据集
dataset.dropna(inplace=True) # 删除缺失值
dataset.drop(columns=dataset.columns[0], axis=1, inplace=True) # 删除第一列，即索引列
# 划分数据集与测试集
X = dataset.drop(columns="Sales", axis=1) # 自变量，广告花费
y = dataset["Sales"] # 因变量，销售额
x_train, x_test, y_train, y_test = sklearn.model_selection.train_test_split(
X, y, test_size=0.3, random_state=42
)
# 特征工程
scaler = StandardScaler() # 标准化特征
x_train = scaler.fit_transform(x_train) # 拟合并转换训练集特征
x_test = scaler.transform(x_test) # 转换测试集特征
# 模型训练
# 正规方程法
lr_model = LinearRegression() # 线性回归模型
lr_model.fit(x_train, y_train) # 拟合模型
print("lr_model.coef_:", lr_model.coef_) # 打印模型系数
print("lr_model.intercept_:", lr_model.intercept_) # 打印模型截距
# 梯度下降法
sgd_model = SGDRegressor()
sgd_model.fit(x_train, y_train)
print("sgd_model.coef_:", sgd_model.coef_) # 打印模型系数
print("sgd_model.intercept_:", sgd_model.intercept_) # 打印模型截距
# 预测
y_pred1 = lr_model.predict(x_test) # 预测测试集
y_pred2 = sgd_model.predict(x_test) # 预测测试集
# 使用均方误差评价模型
mse1 = sklearn.metrics.mean_squared_error(y_test, y_pred1) # 计算均方误差
mse2 = sklearn.metrics.mean_squared_error(y_test, y_pred2) # 计算均方误差
print(f"lr_model MSE: {mse1:.4f}")
print(f"sgd_model MSE: {mse2:.4f}")
浙公网安备 33010602011771号