【机器学习】多元线性回归算法和正规方程解求解

一、摘要

本文围绕多元线性回归的正规方程解展开，为初学者系统介绍了相关基本概念、求解方法、实际应用以及算法封装要点。

首先，深入阐释了正规方程解这一多元线性回归的重要求解方法，同时明确了截距和系数的概念及其差异。通过详细步骤，指导读者如何运用正规方程解决多元线性回归问题，并剖析了将截距和系数分开处理的原因，帮助读者理解其背后的数学逻辑和实际意义。

接着，全面分析了正规方程解的优缺点，使读者清晰了解该方法在不同场景下的适用性。在此基础上，介绍了正规方程解在实际应用中的具体操作方法，增强了读者将理论知识转化为实践的能力。

在算法封装层面，着重解释了在封装多元线性回归算法时，采用截距和系数分开处理方式的原因，为读者在实际编程实现中提供了思路。同时，还提及了 我们自己新建的工程包中新加入的广义线性回归模型，该模型具备强大的通用性，既能解决多元线性回归问题，也能应对简单线性回归问题，拓宽了读者的技术视野。

最后，以 fit 函数的实现过程为切入点，详细讲解了使用正规化方程求解的具体步骤，包括如何创建矩阵 $X_b$，让读者深入了解算法的底层实现细节，为进一步学习和开发奠定坚实基础。

二、多元线性回归介绍

1. 多元线性回归用于处理样本具有多个特征值的问题，常见于真实世界的数据集。

2. 与简单线性回归相比，多元线性回归的样本可能有成千上万个特征。
   

3. 多元线性回归算法的目标：损失函数定义为预测值与真实值之间的差的平方和，旨在最小化该函数。通过找到合适的 θ0 到 θn 的值，使得损失函数最小。
   

4. 多元线性回归的数学表示 ： 为 y = θ0 + θ1x1 + θ2x2 + … + θnxn，其中 θ0 是截距项，θ1 到 θn 是各特征的系数。θ0 到 θn 共 n+1 个参数需要估计。
   
   通过向量化的表示方法将 θ0 到 θn 整理成一个向量 θ，同时将 X 矩阵的每一行添加一个虚拟的第零个特征，值为 1。这样可以将预测值表示为 Y^ = Xθ，其中 Y^ 是预测值向量。
   
   
   
   最终得到的公式如上图所示，这个公式叫做多元线性回归的正规方程解，总结如下图所示：
   

三、正规方程解的求解及代码实现

1. 正规方程解是通过矩阵运算直接求解 θ 的方法，时间复杂度为 O(n^3)，时间复杂度较高，对于大规模数据集可能效率较低。实际应用中，更常用的方法是梯度下降法。
2. 求解过程涉及矩阵求导和求逆运算，具体推导过程较为复杂。
3. 以波士顿房价数据为例来编写程序和测试
   • 在PyChram中新建metrics.py文件，用于实现MSE、RMSE、MAE求解代码

     import numpy as np
     
     
     def mean_squared_error(y_true, y_predict):
         """计算y_true和y_predict之间的MSE"""
         assert len(y_true) == len(y_predict), \
             "the size of y_true must be equal to the size of y_predict"
         return np.sum((y_true - y_predict) ** 2) / len(y_true)
     
     
     def root_mean_squared_error(y_true, y_predict):
         """计算y_true和y_predict之间的RMSE"""
         return np.sqrt(mean_squared_error(y_true, y_predict))
     
     
     def mean_absolute_error(y_true, y_predict):
         """计算y_true和y_predict之间的MAE"""
         assert len(y_true) == len(y_predict), \
             "the size of y_true must be equal to the size of y_predict"
         return np.sum(np.abs(y_true - y_predict)) / len(y_true)
     
     
     def r2_score(y_true, y_predict):
         """计算y_true和y_predict之间的R方值"""
     
         return 1 - mean_squared_error(y_true, y_predict) / np.var(y_true)
     
     

   • 在PyChram中编写正规方程解代码

     import numpy as np
     from metrics import r2_score
     
     
     class LinearRegression:
         # 构造函数
         def __init__(self):
             """初始化Linear Regression模型"""
             self.interception_ = None  # 定义正规方程中的截距
             self.coef_ = None  # 定义正规方程中的系数
             self._theta = None  # 定义正规方程中的内部参数
     
         def fit_normal(self, X_train, y_train):
             """根据训练数据集X_train, y_train训练Linear Regression模型"""
             assert X_train.shape[0] == y_train.shape[0], \
                 "the size of X_train must be equal to the size of y_train"
     
             X_train = np.array(X_train)
             X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
     
             # 检查 X_b 的数据类型
             print("X_b的数据类型:", X_b.dtype)
     
             # 检查 X_b 中是否有非数值元素
             if np.issubdtype(X_b.dtype, np.number):
                 print("X_b 全部为数值类型")
             else:
                 print("X_b 存在非数值类型元素")
     
             # 尝试将 X_b 转换为数值类型
             X_b = X_b.astype(np.float64)
     
             self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
             self.interception_ = self._theta[0]
             self.coef_ = self._theta[1:]
     
             return self
     
         def predict(self, X_predict):
             """给定待预测数据集X_predict, 返回表示X_predict的结果向量"""
             assert self.interception_ is not None and self.coef_ is not None, \
                 "must fit before predict!"
             assert X_predict.shape[1] == len(self.coef_), \
                 "the feature number of X_predict must be equal to X_train"
     
             X_train = np.array(X_predict)
             X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
     
             # 检查 X_b 的数据类型
             print("X_b的数据类型:", X_b.dtype)
     
             # 检查 X_b 中是否有非数值元素
             if np.issubdtype(X_b.dtype, np.number):
                 print("X_b 全部为数值类型")
             else:
                 print("X_b 存在非数值类型元素")
     
             # 尝试将 X_b 转换为数值类型
             X_b = X_b.astype(np.float64)
     
             return X_b.dot(self._theta)
     
         def score(self, X_test, y_test):
             """根据测试数据集 X_test 和 y_test 确定当前模型的准确度"""
             y_predict = self.predict(X_test)
             return r2_score(y_test, y_predict)
     
         def __repr__(self):
             return "LinearRegression()"
     

   • 运行测试，先导入波士顿房价测试数据

     import openml
     import numpy as np
     
     # 从 openml 获取波士顿房价数据集
     dataset = openml.datasets.get_dataset(531)
     X, y, categorical_indicator, attribute_names = dataset.get_data(
         target=dataset.default_target_attribute, dataset_format='dataframe'
     )
     
     # 分布在50那里的一些点，可能不是真实的点，比如问卷调查中通过会设置一些上限点，而往往这些不是真实存在的额点，因此可以去除
     y_normal = y[y < 50.0]
     x_normal = X[y < 50.0]
     
     import sys
     # 替换为你的 PyCharm 工程实际路径
     project_path = 'D:/PycharmProjects/pythonProject/'
     if project_path not in sys.path:
         sys.path.append(project_path)
     
     # 拆分训练集和测试集
     from model_selection import train_test_split
     X_train1, y_train1, X_test1, y_test1 = train_test_split(np.array(x_normal), np.array(y_normal), seed=666)
     

     在Jupyter中先运行，如下图所示：
     
   • 在Jupyter运行测试我们编写的模型

     from LinearRegressionModel import LinearRegression
     reg3 = LinearRegression()
     reg3.fit_normal(X_train1,y_train1)
     

     执行结果：
     
     打印出截距、系数：
     
     利用R方指标来衡量我们的模型最值的效果：
     
     R方值为：0.812979405621282，相比简单线性方程中的预测结果要高出许多。在13维的空间中，我们的模型预测的结果要好得多。