文章目录
- 前言
- 一、逐步回归
- 1. 前进法（Forward Selection）
- 2. 后退法（Backward Elimination）
- 3. 逐步回归法（Stepwise Regression）
- 二、示例
- 三、代码实现----python
前言
- Matlab中逐步回归的实现可以使用 Matlab 的
stepwise函数，本文主要讨论逐步回归如何在 python 中使用。 - 思路参考视频：
https://www.bilibili.com/video/BV1kU4y1R7o2/?spm_id_from=333.1007.top_right_bar_window_history.content.click&vd_source=67471d3a1b4f517b7a7964093e62f7e6
一、逐步回归
- 逐步回归（Stepwise Regression）是一种选择统计模型的技术，用于找到最优模型，即通过添加或移除变量来选择合适的特征。
- 逐步回归主要有三种方法：前进法（Forward Selection）、后退法（Backward Elimination）和逐步回归法（Stepwise Regression）。下面是对这三种方法的简单介绍：
1. 前进法（Forward Selection）
概念：
- 前进法从一个空模型开始，即最初没有任何预测变量。
- 然后逐步添加预测变量，每次添加一个变量，使得新模型的评价指标（如AIC、BIC、R^2等）最优。
- 继续这个过程，直到添加任何更多的变量都不能显著提高模型的性能。
步骤：
- 从空模型开始，不包含任何预测变量。
- 评估每个未加入模型的变量，将使模型性能最优的变量加入模型。
- 重复步骤2，直到添加任何变量都不能显著改善模型。
2. 后退法（Backward Elimination）
概念：
- 后退法从包含所有预测变量的模型开始。
- 然后逐步移除预测变量，每次移除一个变量，使得新模型的评价指标最优。
- 继续这个过程，直到移除任何更多的变量都不能显著提高模型的性能。
步骤：
- 从包含所有可能的预测变量的全模型开始。
- 评估每个变量的显著性，移除最不显著的变量（即对模型贡献最小的变量）。
- 重复步骤2，直到移除任何变量都不能显著改善模型。
3. 逐步回归法（Stepwise Regression）
概念：
- 逐步回归法结合了前进法和后退法，既可以添加变量也可以移除变量。
- 每次步骤既可以是添加一个新变量，也可以是移除一个现有变量，以达到模型性能的最优。
步骤：
- 从空模型开始或包含所有预测变量的模型开始（具体取决于实现方式）。
- 在每一步中，评估所有可能的添加或移除变量的操作。
- 选择对模型性能最优的操作（添加或移除一个变量）。
- 重复步骤2和步骤3，直到添加或移除任何变量都不能显著改善模型。
二、示例
- 水泥凝固时放出的热量
y
y
y 与水泥中 4 种化学成分
x
1
,
x
2
,
x
3
,
x
4
x_1,x_2,x_3,x_4
x1,x2,x3,x4 有关，今测得一组数据如下,试用逐步回归确定一个线性模型，并找出影响水泥凝固时放出热量的必要因素
根据此示例，本文选用后退法选择出影响水泥凝固时放出热量的必要因素。
三、代码实现----python
1. 输入数据
import pandas as pd
import numpy as np
# 数据
x1 = np.array([7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10])
x2 = np.array([26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68])
x3 = np.array([6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8])
x4 = np.array([60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12])
y = np.array([78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 115.9, 83.8, 113.3, 109.4])
# 自变量矩阵
X = pd.DataFrame({'x1': x1, 'x2': x2, 'x3': x3, 'x4': x4})
2. 初始化
本文选用的后退法，所以被选择的因素初始化为包含所有因素，被排出的元素列表为空。
# 初始化未被选中的因素
excluded = list(initial_list)
# 初始化被选中的因素
included = list(set(X.columns) - set(excluded))
3. 评估每个变量的显著性
本文中使用 P P P 值的大小评估每个变量的显著性。
步骤：
- 拟合回归模型
- 获取所有特征的P值
- 找到最大的P值及其对应的特征。
# 拟合选中的因素
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
# 得出被选中的因素里P值最大的
p_values = model.pvalues
# 排除常数项（截距项）的P值
p_values = p_values.drop('const')
# 找到最大的P值及其对应的特征
max_p_value = p_values.max()
print("最大的P值为:",max_p_value)
4. 判断函数退出的标志
直到被选中的因素拟合后得到的最大 P P P 值小于 α ( 0.05 ) \alpha(0.05) α(0.05)，意味着移除任何变量都不能显著改善模型，函数退出。
def stepwise_selection(X, y,
initial_list = [],
threshold = 0.05,
mark = True):
# 初始化未被选中的因素
excluded = list(initial_list)
# 初始化被选中的因素
included = list(set(X.columns) - set(excluded))
while mark:
# 拟合选中的因素
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
# 得出被选中的因素里P值最大的
p_values = model.pvalues
# 排除常数项（截距项）的P值
p_values = p_values.drop('const')
# 找到最大的P值及其对应的特征
max_p_value = p_values.max()
print("最大的P值为:",max_p_value)
if max_p_value < threshold:
mark = False
print("最终模型:")
print(model.summary())
# 获取回归系数
b = model.params
else:
max_p_feature = p_values.idxmax()
print("最大的P值对应的特征为:",max_p_feature)
# 从被选中的因素中去除
included.remove(max_p_feature)
print("更新后的因素为:",included)
return included, b
5. 逐步回归的完整代码
import pandas as pd
import numpy as np
import statsmodels.api as sm
# 数据
x1 = np.array([7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10])
x2 = np.array([26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68])
x3 = np.array([6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8])
x4 = np.array([60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12])
y = np.array([78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 115.9, 83.8, 113.3, 109.4])
# 自变量矩阵
X = pd.DataFrame({'x1': x1, 'x2': x2, 'x3': x3, 'x4': x4})
def stepwise_selection(X, y,
initial_list=[],
threshold=0.05,
mark = True):
# 初始化未被选中的因素
excluded = list(initial_list)
# 初始化被选中的因素
included = list(set(X.columns) - set(excluded))
while mark:
# 拟合选中的因素
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
# 得出被选中的因素里P值最大的
p_values = model.pvalues
# 排除常数项（截距项）的P值
p_values = p_values.drop('const')
# 找到最大的P值及其对应的特征
max_p_value = p_values.max()
print("最大的P值为:",max_p_value)
if max_p_value < threshold:
mark = False
print("最终模型:")
print(model.summary())
# 获取回归系数
b = model.params
else:
max_p_feature = p_values.idxmax()
print("最大的P值对应的特征为:",max_p_feature)
# 从被选中的因素中去除
included.remove(max_p_feature)
print("更新后的因素为:",included)
return included, b
result, b = stepwise_selection(X, y)
print(result)
print(b)
运行结果：
评估变量显著性的过程：
最终的模型：
回归系数：
6. 生成三维图的代码
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
# 计算Z的值
Z = b['const'] + b['x1'] * x1 + b['x2'] * x2
# 创建一个新的3D绘图对象
fig = plt.figure()
'''
111 的具体含义是：
第一位 1:整个图形只有 1 行。
第二位 1:整个图形只有 1 列。
第三位 1:子图在这个 1x1 网格中的第 1 个位置。
'''
ax = fig.add_subplot(111, projection='3d')
# 绘制散点图
ax.scatter(x1, x2, Z, c='r', marker='o', label='Data Points')
# 创建网格以绘制曲面
x1_range = np.linspace(min(x1), max(x1), 100)
x2_range = np.linspace(min(x2), max(x2), 100)
x1_grid, x2_grid = np.meshgrid(x1_range, x2_range)
# 计算曲面上的Z值
Z_surface = b['const'] + b['x1'] * x1_grid + b['x2'] * x2_grid
# 绘制曲面图
'''
alpha
意义:alpha 参数用于设置曲面的透明度。
取值范围:alpha 的取值范围是 0 到 1 之间。
0 表示完全透明，即不可见。
1 表示完全不透明。
作用：通过调整 alpha 参数，你可以在同一视图中更好地叠加多个图形，使得它们不会完全遮挡对方。
cmap
意义:cmap 参数用于设置曲面的颜色映射(colormap)。
常见的颜色映射:viridis, plasma, inferno, magma, cividis, jet, rainbow, coolwarm, hot 等。
作用：颜色映射用于根据 Z 值来着色曲面，帮助更清晰地展示高度或强度的变化。
'''
ax.plot_surface(x1_grid, x2_grid, Z_surface, alpha=0.5, cmap='viridis')
# 设置标签
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_zlabel('Z')
# 添加图例
ax.legend()
# 显示图形
plt.show()
运行结果：
