线性回归模型的建立与求解——以波士顿房价问题为例

线性回归模型

模型引出

• 在统计学中，回归分析（regression analysis)指的是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。回归分析按照涉及的变量的多少，分为一元回归和多元回归分析；按照因变量的多少，可分为简单回归分析和多重回归分析；按照自变量和因变量之间的关系类型，可分为线性回归分析和非线性回归分析。
• 在大数据分析中，回归分析是一种预测性的建模技术，它研究的是因变量（目标）和自变量（预测器）之间的关系。这种技术通常用于预测分析以及发现变量之间的因果关系。例如，司机的鲁莽驾驶与道路交通事故数量之间的关系，最好的研究方法就是回归。

模型原理

模型表示：线性回归模型表示为：
$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\dots +{\beta }_n{x}_n+ϵ$$
其中，$y$表示因变量（待预测的变量），$x_1,x_2,...,x_n$是自变量，${\beta }_0$表示截距，${\beta }_1,{\beta }_2...{\beta }_n$表示自变量参数，$ϵ$表示误差项。

参数估计：对于线性回归，关键在于求解参数，常用高斯提出的最小二乘法，具体步骤为：

1. 定义残差平法和(RSS)

   残差平方和（RSS）是预测值与实际值之间差值的平方和，表示模型预测的不准确程度。定义为：
   $$\mathrm{RSS}=\sum _{i=1}^m{\left(y_i-{\hat{y}}_i\right)}^2$$
   其中：

   • $y_i$ 是第 $i$ 个观测值的实际值。

   • ${\hat{y}}_i$是第$i$个观测值的预测值，即：
     $${\hat{y}}_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_n{x}_{in}$$
     因此，RSS 可以进一步展开为：
     $$\mathrm{RSS}=\sum _{i=1}^m{\left(y_i-({\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_n{x}_{in})\right)}^2$$

2. 最小化RS

   为了找到最优的回归系数${\beta }_0,{\beta }_1...{\beta }_n$，我们需要最小化RSS。最小化RSS的过程就是通过导数找到使RSS最小的回归系数。

   将RSS对每个回归系数${\beta }_j$求导，并使导数为零，以找到最小点：
   $$\frac{\partial \mathrm{RSS}}{\partial {\beta }_j}=-2\sum _{i=1}^m\left(y_i-{\hat{y}}_i\right)x_{ij}=0$$
   其中$j=0,1,2,...,n$对应每个回归系数。

   求导后，得到一组关于${\beta }_j$的方程，称为正则方程：
   $$\sum _{i=1}^m{y}_i\cdot x_{ij}=\sum _{i=1}^m\left({\beta }_0{x}_{i0}+{\beta }_1{x}_{i1}+\cdots +{\beta }_n{x}_{in}\right)x_{ij}$$
   其中$x_{i0}$对应于截距项${\beta }_0$。

3. 参数估计公式

   将正则方程组写成矩阵形式，可以更方便地求解回归系数。首先定义以下矩阵：

   • 观测值向量$Y$：
     $$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$$

   • 设计矩阵$X$:
     $$\mathbf{X}=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \dots & x_{1n}\\ 1& x_{21}& x_{22}& \dots & x_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{m1}& x_{m2}& \dots & x_{mn}\end{array}\right]$$

   • 回归系数向量$\beta$:
     $$\mathbf{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_n\end{array}\right]$$

   • 误差向量$ϵ$:
     $$\mathbf{ϵ}=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_m\end{array}\right]$$

于是，线性回归模型可以表示为：
$$\mathbf{y}=\mathbf{X}\beta +\mathbf{ϵ}$$
使用最小二乘法估计 $\beta$ 的解为：
$$\hat{\beta }=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$$
这个公式是通过最小化残差平方和RSS 得到的，提供了回归系数的最佳估计值。利用这个公式，我们可以计算出每个自变量的系数，并使用这些系数对新的数据进行预测。

使用线性回归问题解决波士顿房价问题

问题描述

预测波士顿不同地区的房价（MEDV）是一个经典的回归分析问题。我们将使用波士顿房价数据集中的特征变量，通过线性回归模型来预测房价。线性回归模型可以通过最小二乘法（OLS）来估计回归系数。

数据集包括14个变量，其中13个自变量（特征）和1个因变量（房价中位数，MEDV）。自变量包括城镇人均犯罪率（CRIM）、住宅用地所占比例（ZN）、非住宅用地所占比例（INDUS）等。目标是建立一个多元线性回归模型，以便使用这些自变量来预测房价。具体信息如下表所示

变量	描述
CRIM	城镇人均犯罪率
ZN	住宅用地所占比例
INDUS	城镇中非住宅用地所占比例
CHAS	虚拟变量,用于回归分析
NOX	环保指数
RM	每栋住宅的房间数
AGE	1940 年以前建成的自住单位的比例
DIS	距离 5 个波士顿的就业中心的加权距离
RAD	距离高速公路的便利指数
TAX	每一万美元的不动产税率
PTRATIO	城镇中的教师学生比例
B	城镇中的黑人比例
LSTAT	地区中有多少房东属于低收入人群
MEDV	自住房屋房价中位数（也就是均价）

数据集链接：https://pan.baidu.com/s/1JV3-tbWLLBRIrO0APNMeWw 提取码：1111

建立数学模型

线性回归模型可以表示为：
$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_n{x}_n+ϵ$$
其中：

• $y$是因变量（MEDV，即房价中位数）。
• $x_1,x_2,...,x_n$是自变量（13个特征变量）。
• ${\beta }_0$是截距。
• ${\beta }_1,\beta 2,...,\beta n$是回归系数。
• $ϵ$是误差项，表示模型与实际数据之间的差异.

最小二乘法求解

最小二乘法的目的是通过最小化残差平方和（RSS）来找到最优的回归系数${\beta }_j$。

• 残差平方和（RSS）定义为：
  $$\mathrm{RSS}=\sum _{i=1}^m{\left(y_i-{\hat{y}}_i\right)}^2=\sum _{i=1}^m{\left(y_i-\left({\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_n{x}_{in}\right)\right)}^2$$

• 通过对每个 ${\beta }_j$求导，并设导数为零，可以得到一组关于 ${\beta }_j$ 的方程，称为正则方程。将正则方程组写成矩阵形式，可以用以下公式计算回归系数：
  $$\hat{\beta }=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$$
  其中:

• $X$是$m\times (n+1)$的设计矩阵，包含截距项和各个自变量。

• $Y$是$m\times 1$的因变量向量。

使用Python语言进行模型的求解

数据导入

#导入Python常用数据分析的库

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()    #设置画图空间为 Seaborn 默认风格

names=['CRIM','ZN','INDUS','CHAS','NOX','RM','GE','DIS','RAD','TAX','PRTATIO','B','LSTAT','PRICE']

boston=pd.read_csv("boston_house_prices (1).csv",names=names)
boston.head(10)

• pd.read_csv 是 pandas 提供的一个函数，用于从 CSV 文件中读取数据。"boston_house_prices (1).csv" 是数据文件的路径。
• names=names 参数指定了读取数据时使用的列名列表。如果 CSV 文件中没有包含列名，使用 names 参数会将这些列名作为 DataFrame 的列名。

查看数据信息

boston.info()

info()是 pandas DataFrame 的一个方法，用于获取关于 DataFrame 的简要信息。运行这行代码将显示以下内容：

• DataFrame的行数和列数：包括数据框的总行数（数据条目数）和列数。
• 每一列的名称：列的标题，通常是通过 CSV 文件中定义的或在创建 DataFrame 时指定的.
• 每一列的数据类型：如 int64, float64, object（通常表示字符串），等。
• 每列的非空值数量：显示每列中的非空（非缺失）值的数量，可以用来识别缺失数据的情况。
• 缺失值数量：可以通过非空值数量和总行数的差值间接了解每列的缺失值数量。
• 内存消耗：DataFrame 占用的内存大小（以字节为单位），帮助评估数据集的内存需求。

具体结果如下图所示

由上图可知，数据中没有缺失值，不需要进行缺失值的处理。

如果数据中存在缺失值，常用的缺失值处理方法

• 删除缺失值

  df.dropna()  # 删除所有含有缺失值的行
  df.dropna(axis=1)  # 删除所有含有缺失值的列
  

• 填充缺失值

  df.fillna(value=0)  # 用0填充所有缺失值
  df.fillna({'column_name': value})  # 用特定值填充某一列的缺失值
  df.fillna(method='ffill')  # 用前一个值填充缺失值
  df.fillna(method='bfill')  # 用后一个值填充缺失值
  df.fillna(df.mean())  # 用均值填充缺失值
  df.fillna(df.median())  # 用中位数填充缺失值
  df.fillna(df.mode().iloc[0])  # 用众数填充缺失值
  

描述性数据分析

boston.describe()

• describe方法可以得出下表所示数据

统计量	说明
count	每列的非缺失值数量
mean	每列数据的均值
std	每列数据的标准差，衡量数据的离散程度
min	每列数据的最小值
25%	25 百分位数，即数据中 25% 的值小于或等于此值
50%	50 百分位数，即中位数
75%	75 百分位数，即数据中 75% 的值小于或等于此值
max	每列数据的最大值

查看各字段的相关性

corrboston = boston.corr()
corrboston

plt.figure(figsize=(10,10))    #设置画布
sns.heatmap(corrboston,annot=True,cmap='RdGy')
plt.savefig('热力图.svg', format='svg')
plt.show()

• boston.corr()：计算 boston 中所有数值列之间的相关系数矩阵。相关系数矩阵显示了每对列之间的相关性。

• sns.heatmap(corrboston, annot=True, cmap='RdGy')：使用 Seaborn 库绘制相关系数矩阵的热力图。corrboston：要绘制的相关系数矩阵数据。annot=True：在热力图的每个单元格上显示相关系数的数值。cmap='RdGy'：指定热力图的颜色映射方案。'RdGy' 是一个红色到灰色的渐变色图。

• plt.savefig('热力图.svg', format='svg')：将生成的图形保存为 SVG 格式的文件，文件名为 '热力图.svg'。SVG（可缩放矢量图形）格式适合高质量打印和矢量图处理。

查看是否穿过查尔斯河对房价的影响

bostonCHAS = boston[['CHAS','PRICE']]    #先将CHAS和PRICE两列数据取出

bostonCHAS1=bostonCHAS.pivot_table(values='PRICE',    #计算的值
                               index='CHAS',       #透视的行，分组的依据
                               aggfunc='mean')            #聚合函数

# 对透视表进行降序排列
bostonCHAS1 = bostonCHAS1.sort_values(by='PRICE',     # 排序依据
                        ascending=False                 # 是否升序排列
                       )

bostonCHAS1

根据结果可以得出，被查尔斯河穿过的豪宅比没被穿过的豪宅价格更高。

各字段与价格的散点图

x_data = boston[['CRIM','ZN','INDUS','CHAS','NOX','RM','GE','DIS','RAD','TAX','PRTATIO','B','LSTAT']] # 导入所有特征变量
y_data = boston[['PRICE']] # 导入目标值（房价）

plt.figure(figsize=(18,10))

for i in range(13):
    plt.subplot(4,4,i+1)
    plt.scatter(x_data.values[:,i],y_data,s = 5)    #.values将DataFrame对象X_df转成ndarray数组
    plt.xlabel(names[i])
    plt.ylabel('Price')
    plt.title(str(i+1)+'. '+names[i]+' - Price')  
    
plt.tight_layout()
plt.show()

通过上图可以看到不是所有的字段与价格都有较强的相关关系，但本例中不涉及多元线性回归的向后删除，仅做最简单的多元性性回归的分析处理。

预测性数据分析

选取线性回归字段

import statsmodels.api as sm
# 导入 statsmodels 库，并使用别名 sm。statsmodels 提供了用于统计建模和分析的功能。

# 复制数据集，以避免对原始数据集进行修改
lr_house_price = boston.copy()

# 提取目标变量 'PRICE'，并将其赋值给 y
y = lr_house_price['PRICE']

# 删除 'PRICE' 列，得到包含所有自变量的特征数据 X
X = lr_house_price.drop('PRICE', axis=1)

# 计算自变量之间的相关系数矩阵，并取得其绝对值。然后判断是否有高于 0.8 的相关性。
# 这个步骤用于检查多重共线性问题，如果某些自变量之间的相关性过高，可能需要进一步处理。
high_corr = X.corr().abs() > 0.8

# 在自变量数据 X 中添加一列常数项（全为1的列），以便估计回归模型中的截距项。
# statsmodels 的线性回归模型默认不包括截距项，因此需要手动添加这一列。
X = sm.add_constant(X)

调用OLS函数，利用最小二乘法来得到线性回归模型的参数值

model=sm.OLS(y,X).fit()
model.summary()

• model = sm.OLS(y, X).fit()：使用 statsmodels 库中的 OLS 类来创建一个普通最小二乘法（OLS）线性回归模型。这里的 y 是因变量（房价中位数），X 是包含自变量的设计矩阵（包括常数项）。.fit() 方法用于拟合这个模型，返回一个包含拟合结果的 RegressionResults 对象。

• model.summary()：调用 RegressionResults 对象的 summary() 方法，生成一个详细的回归结果汇总。这个方法会输出模型的各种统计信息和评估指标，用于帮助解释和评估回归模型的表现。可以得到下表所示的数据：

  Dep. Variable:	PRICE	R-squared:	0.741
  Model:	OLS	Adj. R-squared:	0.734
  Method:	Least Squares	F-statistic:	108.1
  Date:	Fri, 26 Jul 2024	Prob (F-statistic):	6.72e-135
  Time:	09:08:22	Log-Likelihood:	-1498.8
  No. Observations:	506	AIC:	3026.
  Df Residuals:	492	BIC:	3085.
  Df Model:	13
  Covariance Type:	nonrobust

  	coef	std err	t	P>|t|	[0.025	0.975]
  const	36.4595	5.103	7.144	0.000	26.432	46.487
  CRIM	-0.1080	0.033	-3.287	0.001	-0.173	-0.043
  ZN	0.0464	0.014	3.382	0.001	0.019	0.073
  INDUS	0.0206	0.061	0.334	0.738	-0.100	0.141
  CHAS	2.6867	0.862	3.118	0.002	0.994	4.380
  NOX	-17.7666	3.820	-4.651	0.000	-25.272	-10.262
  RM	3.8099	0.418	9.116	0.000	2.989	4.631
  GE	0.0007	0.013	0.052	0.958	-0.025	0.027
  DIS	-1.4756	0.199	-7.398	0.000	-1.867	-1.084
  RAD	0.3060	0.066	4.613	0.000	0.176	0.436
  TAX	-0.0123	0.004	-3.280	0.001	-0.020	-0.005
  PRTATIO	-0.9527	0.131	-7.283	0.000	-1.210	-0.696
  B	0.0093	0.003	3.467	0.001	0.004	0.015
  LSTAT	-0.5248	0.051	-10.347	0.000	-0.624	-0.425

  Omnibus:	178.041	Durbin-Watson:	1.078
  Prob(Omnibus):	0.000	Jarque-Bera (JB):	783.126
  Skew:	1.521	Prob(JB):	8.84e-171
  Kurtosis:	8.281	Cond. No.	1.51e+04

各项数据的含义

名称	作用
决定系数（R-squared）	衡量模型对因变量的解释能力。值越接近1，说明模型对数据的解释能力越强。
调整后的决定系数（Adj. R-squared）	考虑了模型中的自变量个数，更准确地反映模型的解释能力。
回归系数（coef）	表示每个自变量对因变量的影响程度。正值或负值反映了自变量与因变量的关系。
标准误（std err）	回归系数估计的标准误差，反映系数估计的不确定性。标准误越小，说明系数估计越精确。
t 值（t）	回归系数与其标准误的比值，用于判断系数的显著性。
p值（p>|t|）	检验回归系数是否显著。p 值小于 0.05 通常表示该变量对因变量有显著影响。
置信区间（95% [0.025, 0.975]）	给出回归系数的估计范围。区间越窄，说明对系数的估计越精确。
F 统计量（F-statistic）	测试模型整体的显著性。衡量所有自变量对因变量的共同影响。
F 统计量的 p 值（Prob (F-statistic)）	用于检验模型整体是否显著。如果 p 值小于 0.05，说明模型中的至少一个自变量对因变量有显著影响。
Durbin-Watson 统计量	检验残差的自相关性。值在 1.5 到 2.5 之间通常表明残差自相关性较低。
AIC（Akaike Information Criterion）	用于模型选择，数值越小，模型越好。
BIC（Bayesian Information Criterion）	用于模型选择，数值越小，模型越好。
Omnibus	测试残差的正态性。Omnibus 统计量及其 p 值用于检验残差是否符合正态分布。
Jarque-Bera (JB)	测试残差的正态性。Jarque-Bera 统计量及其 p 值用于检验残差是否符合正态分布。
偏度（Skew）	衡量残差分布的偏斜程度。值为0表示正态分布。
峰度（Kurtosis）	衡量残差分布的尖峭程度。值为0表示正态分布。
条件数（Condition Number）	检验多重共线性。条件数非常大可能存在多重共线性问题。

重点看一下 P>|t| 这列数据，我们发现有两项数据INDUS和GE的值比较大，意味着该自变量对因变量的影响不显著，可能出现多重共线性等问题，这里我们把这两项数据去掉，再进行一次拟合。

X=X.drop(['INDUS','GE'],axis=1)
model=sm.OLS(y,X).fit()
model.summary()

根据结果可知，拟合效果有了一定的提高。

预测

# 预测
predictions = model.predict(X)
#绘制真实值与预测值比较的图像
from matplotlib.font_manager import FontProperties

# 设置中文字体
font = FontProperties(fname='C:/Windows/Fonts/SimHei.ttf')  # 修改为你安装的字体路径

# 比较预测值与真实值，并绘制折线图
plt.figure(figsize=(12, 6))
plt.plot(y.values, label='真实值', marker='o')
plt.plot(predictions, label='预测值', marker='x')
plt.title('真实值与预测值比较', fontproperties=font)
plt.xlabel('样本索引', fontproperties=font)
plt.ylabel('房价', fontproperties=font)
plt.legend(prop=font)

# 保存为矢量图（SVG 格式）
plt.savefig('price_comparison.svg', format='svg')

plt.show()

根据图形可知，拟合效果良好。