线性回归模型
模型引出
- 在统计学中,回归分析(regression analysis)指的是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。回归分析按照涉及的变量的多少,分为一元回归和多元回归分析;按照因变量的多少,可分为简单回归分析和多重回归分析;按照自变量和因变量之间的关系类型,可分为线性回归分析和非线性回归分析。
- 在大数据分析中,回归分析是一种预测性的建模技术,它研究的是因变量(目标)和自变量(预测器)之间的关系。这种技术通常用于预测分析以及发现变量之间的因果关系。例如,司机的鲁莽驾驶与道路交通事故数量之间的关系,最好的研究方法就是回归。
模型原理
模型表示:线性回归模型表示为:
y
=
β
0
+
β
1
x
1
+
β
2
x
2
+
…
+
β
n
x
n
+
ϵ
y=\beta_0+\beta_1x_1+\beta_2x_2+\ldots+\beta_nx_n+\epsilon
y=β0+β1x1+β2x2+…+βnxn+ϵ
其中,
y
y
y表示因变量(待预测的变量),
x
1
,
x
2
,
.
.
.
,
x
n
x_{1},x_{2},...,x_{n}
x1,x2,...,xn是自变量,
β
0
\beta_0
β0表示截距,
β
1
,
β
2
.
.
.
β
n
\beta_{1},\beta_{2}...\beta_n
β1,β2...βn表示自变量参数,
ϵ
\epsilon
ϵ表示误差项。
参数估计:对于线性回归,关键在于求解参数,常用高斯提出的最小二乘法,具体步骤为:
-
定义残差平法和(RSS)
残差平方和(RSS)是预测值与实际值之间差值的平方和,表示模型预测的不准确程度。定义为:
R S S = ∑ i = 1 m ( y i − y ^ i ) 2 \mathrm{RSS}=\sum_{i=1}^m\left(y_i-\hat{y}_i\right)^2 RSS=i=1∑m(yi−y^i)2
其中:-
y i y_i yi 是第 i i i 个观测值的实际值。
-
y ^ i \hat{y}_i y^i是第 i i i个观测值的预测值,即:
y ^ i = β 0 + β 1 x i 1 + β 2 x i 2 + ⋯ + β n x i n \hat{y}_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_nx_{in} y^i=β0+β1xi1+β2xi2+⋯+βnxin
因此,RSS 可以进一步展开为:
R S S = ∑ i = 1 m ( y i − ( β 0 + β 1 x i 1 + β 2 x i 2 + ⋯ + β n x i n ) ) 2 \mathrm{RSS}=\sum_{i=1}^m\left(y_i-(\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_nx_{in})\right)^2 RSS=i=1∑m(yi−(β0+β1xi1+β2xi2+⋯+βnxin))2
-
-
最小化RS
为了找到最优的回归系数 β 0 , β 1 . . . β n \beta_{0},\beta_{1}...\beta_n β0,β1...βn,我们需要最小化RSS。最小化RSS的过程就是通过导数找到使RSS最小的回归系数。
将RSS对每个回归系数 β j \beta_j βj求导,并使导数为零,以找到最小点:
∂ R S S ∂ β j = − 2 ∑ i = 1 m ( y i − y ^ i ) x i j = 0 \frac{\partial\mathrm{RSS}}{\partial\beta_j}=-2\sum_{i=1}^m\left(y_i-\hat{y}_i\right)x_{ij}=0 ∂βj∂RSS=−2i=1∑m(yi−y^i)xij=0
其中 j = 0 , 1 , 2 , . . . , n j=0,1,2,...,n j=0,1,2,...,n对应每个回归系数。求导后,得到一组关于 β j \beta_{j} βj的方程,称为正则方程:
∑ i = 1 m y i ⋅ x i j = ∑ i = 1 m ( β 0 x i 0 + β 1 x i 1 + ⋯ + β n x i n ) x i j \sum_{i=1}^my_i\cdot x_{ij}=\sum_{i=1}^m\left(\beta_0x_{i0}+\beta_1x_{i1}+\cdots+\beta_nx_{in}\right)x_{ij} i=1∑myi⋅xij=i=1∑m(β0xi0+β1xi1+⋯+βnxin)xij
其中 x i 0 x_{i0} xi0对应于截距项 β 0 \beta_{0} β0。 -
参数估计公式
将正则方程组写成矩阵形式,可以更方便地求解回归系数。首先定义以下矩阵:
-
观测值向量 Y Y Y:
y = [ y 1 y 2 ⋮ y m ] \mathbf{y}=\begin{bmatrix}y_1\\y_2\\\vdots\\y_m\end{bmatrix} y= y1y2⋮ym -
设计矩阵 X X X:
X = [ 1 x 11 x 12 … x 1 n 1 x 21 x 22 … x 2 n ⋮ ⋮ ⋮ ⋱ ⋮ 1 x m 1 x m 2 … x m n ] \mathbf{X}=\begin{bmatrix}1&x_{11}&x_{12}&\dots&x_{1n}\\1&x_{21}&x_{22}&\dots&x_{2n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_{m1}&x_{m2}&\dots&x_{mn}\end{bmatrix} X= 11⋮1x11x21⋮xm1x12x22⋮xm2……⋱…x1nx2n⋮xmn -
回归系数向量 β \beta β:
β = [ β 0 β 1 ⋮ β n ] \boldsymbol\beta=\begin{bmatrix}\beta_0\\\beta_1\\\vdots\\\beta_n\end{bmatrix} β= β0β1⋮βn -
误差向量 ϵ \epsilon ϵ:
ϵ = [ ϵ 1 ϵ 2 ⋮ ϵ m ] \boldsymbol{\epsilon}=\begin{bmatrix}\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_m\end{bmatrix} ϵ= ϵ1ϵ2⋮ϵm
-
于是,线性回归模型可以表示为:
y
=
X
β
+
ϵ
\mathbf{y}=\mathbf{X}\beta+\boldsymbol{\epsilon}
y=Xβ+ϵ
使用最小二乘法估计
β
\beta
β 的解为:
β
^
=
(
X
⊤
X
)
−
1
X
⊤
y
\mathbf{\hat{\beta}}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}
β^=(X⊤X)−1X⊤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
=
β
0
+
β
1
x
1
+
β
2
x
2
+
⋯
+
β
n
x
n
+
ϵ
y=\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_nx_n+\epsilon
y=β0+β1x1+β2x2+⋯+βnxn+ϵ
其中:
- y y y是因变量(MEDV,即房价中位数)。
- x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1,x2,...,xn是自变量(13个特征变量)。
- β 0 \beta_{0} β0是截距。
- β 1 , β 2 , . . . , β n \beta_{1},\beta{2},...,\beta{n} β1,β2,...,βn是回归系数。
- ϵ \epsilon ϵ是误差项,表示模型与实际数据之间的差异.
最小二乘法求解
最小二乘法的目的是通过最小化残差平方和(RSS)来找到最优的回归系数 β j \beta_{j} βj。
-
残差平方和(RSS)定义为:
R S S = ∑ i = 1 m ( y i − y ^ i ) 2 = ∑ i = 1 m ( y i − ( β 0 + β 1 x i 1 + β 2 x i 2 + ⋯ + β n x i n ) ) 2 \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_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_nx_{in}\right)\right)^2 RSS=i=1∑m(yi−y^i)2=i=1∑m(yi−(β0+β1xi1+β2xi2+⋯+βnxin))2 -
通过对每个 β j \beta_{j} βj求导,并设导数为零,可以得到一组关于 β j \beta_{j} βj 的方程,称为正则方程。将正则方程组写成矩阵形式,可以用以下公式计算回归系数:
β ^ = ( X ⊤ X ) − 1 X ⊤ y \mathbf{\hat{\beta}}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} β^=(X⊤X)−1X⊤y
其中: -
X X X是 m × ( n + 1 ) m\times(n+1) m×(n+1)的设计矩阵,包含截距项和各个自变量。
-
Y Y Y是 m × 1 m\times1 m×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()
根据图形可知,拟合效果良好。