文章目录
前言
通过模型算法,熟练对Matlab和python的应用。
学习视频链接:
https://www.bilibili.com/video/BV1EK41187QF?p=42&vd_source=67471d3a1b4f517b7a7964093e62f7e6
一、一元线性回归
- 一元线性回归是一种统计方法,用于理解和建模两个变量之间的关系。具体来说,它是用来估计一个因变量 y y y 和一个自变量 x x x 之间线性关系的工具。其目标是找到一条直线,能够最好地预测 y y y 的值,基于给定的 x x x 值。这条直线通常被称为 “回归线” 。
回归分析的一般步骤
一元线性回归的基本形式
y = β 0 + β 1 x + ϵ y=\beta_0+\beta_1x+\epsilon y=β0+β1x+ϵ
- y y y 是因变量(被解释变量)。
- x x x 是自变量(解释变量)。
- β 0 \beta_0 β0 是截距(当 x = 0 x = 0 x=0 时, y y y 的预测值)。
- β 1 \beta_1 β1 是斜率(表示 x x x 每增加一个单位, y y y 变化的量)。
- ϵ \epsilon ϵ 是误差项,表示因变量的实际值和预测值之间的差异。
回归方程
- 描述因变量y的期望值如何依赖于自变量x的方程称为回归方程。根据对一元线性回归模型的假设,可以得到它的回归方程为:
E ( y ) = β 0 + β 1 x E(y)=\beta_0+\beta_1x E(y)=β0+β1x - 如果回归方程中的参数已知,对于一个给定的x值,利用回归方程就能计算出y的期望值
- 用样本统计量代替回归方程中的未知参数,就得到估计的回归方程,简称回归直线
参数的最小二乘法估计
- 对于回归直线,关键在于求解参数,常用高斯提出的最小二乘法,它是使因变量的观察值
y
y
y 与估计值之间的离差平方和达到最小来求解的。
Q = ∑ ( y − y ^ ) 2 = ∑ ( y − β ^ 0 − β ^ 1 x ) 2 Q=\sum(y-\hat{y} )^2=\sum(y-\hat{\beta}_0-\hat{\beta}_1x)^2 Q=∑(y−y^)2=∑(y−β^0−β^1x)2
展开可得:
Q = ∑ ( y − y ^ ) 2 = ∑ y 2 + n β ^ 0 2 + β ^ 1 2 ∑ x 2 + 2 β ^ 0 β ^ 1 ∑ x − 2 β ^ 0 ∑ y − 2 β ^ 1 ∑ x y Q=\sum(y-\hat{y})^{2}=\sum y^{2}+n\hat{\beta}_{0}^{2}+\hat{\beta}_{1}^{2}\sum x^{2}+2\hat{\beta}_{0}\:\hat{\beta}_{1}\sum x-2\hat{\beta}_{0}\:\sum y-2\hat{\beta}_{1}\sum xy Q=∑(y−y^)2=∑y2+nβ^02+β^12∑x2+2β^0β^1∑x−2β^0∑y−2β^1∑xy
求偏导可得:
{ ∑ y = n β ^ 0 + β ^ 1 Σ x ∑ x y = 2 β ^ 0 Σ x + β ^ 1 Σ x 2 \begin{cases}\sum y=n\hat{\beta}_0+\hat{\beta}_1\Sigma x\\\sum xy=2\hat{\beta}_0\Sigma x+\hat{\beta}_1\Sigma x^2\end{cases} {∑y=nβ^0+β^1Σx∑xy=2β^0Σx+β^1Σx2
{ β ^ 1 = n ∑ x y − ∑ x ∑ y n ∑ x 2 − ( ∑ x ) 2 β ^ 0 = y ˉ − β ^ 1 x ˉ \begin{cases}\hat{\beta}_1=\dfrac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x)^2}\\\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}\end{cases} ⎩ ⎨ ⎧β^1=n∑x2−(∑x)2n∑xy−∑x∑yβ^0=yˉ−β^1xˉ
即可得到 β ^ 0 \hat{\beta}_0 β^0及 β ^ 1 \hat{\beta}_1 β^1
对回归方程的各种检验
估计标准误差的计算
- 为了度量回归方程的可靠性,通常计算估计标准误差,它度量观察值回绕着回归直线的变化程度或分散程度
- 估计平均误差:
S e = Σ ( y − y ^ ) 2 n − 2 S_e=\sqrt{\frac{\Sigma(y-\hat{y})^2}{n-2}} Se=n−2Σ(y−y^)2 - 估计标准误差越大,则数据点围绕回归直线的分散程度就越大,回归方程代表性就越小
- 估计标准误差越小,则数据点围绕回归直线的分散程度越小,回归方程的代表性越大,可靠性越高
回归直线的拟合优度
- 回归直线与各观测点的接近程度称为回归直线对数据的拟合优度
- 总平方和 (TSS):反映因变量的 n n n个观察值与其均值的总离差 T S S = ∑ y i 2 = ∑ ( y i − y ˉ i ) 2 TSS=\sum y_i^2=\sum(y_i-\bar{y}_i)^2 TSS=∑yi2=∑(yi−yˉi)2
- 回归平方和 (ESS):反映了 y y y 的总变差中,由于 x x x 与 y y y 之间的线性关系引起的 y y y 的变化部分 E S S = ∑ y ^ i 2 = ∑ ( y ^ i − y ˉ i ) 2 ESS=\sum\hat{y}_i^2=\sum(\hat{y}_i-\bar{y}_i)^2 ESS=∑y^i2=∑(y^i−yˉi)2
- 残差平方和(RSS):反映了除了 x x x对 y y y的线性影响之外的其他因素对 y y y变差的作用,是不能由回归直线来解释的y的变差部分 R S S = ∑ e i 2 = ∑ ( y i − y ^ i ) 2 RSS=\sum e_i^2=\sum(y_i-\hat{y}_i)^2 RSS=∑ei2=∑(yi−y^i)2
- 总平方和可以分解为回归平方和、残差平方和两部分: T S S = E S S + R S S TSS=ESS+RSS TSS=ESS+RSS
判定系数
- 回归平方和占总平方和的比例,用
R
2
R^2
R2 表示,其值在0到1之间
- R 2 = 0 R^{2}=0 R2=0:说明 y y y 的变化与 x x x 无关, x x x 完全无助于解释 y y y 的变差
- R 2 = 1 R^{2}=1 R2=1:说明残差平方和为0,拟合是完全的, y y y 的变化只与 x x x 有关
显著性检验
- 显著性检验的主要目的是根据所建立的估计方程用自变量 x x x 来估计或预测因变量 y y y 的取值,当建立了估计方程后,还不能马上进行估计或预测,因为该估计方程是根据样本数据得到的,它是否真实的反映了变量 x x x 和 y y y 之间的关系,则需要通过检验后才能证实。
- 根据样本数据拟合回归方程时,实际上就已经假定变量 x x x 和 ν \nu ν 之间存在着线性关系,并假定误差项是一个服从正态分布的随机变量,且具有相同的方差,但这些假设是否成立需要检验
- 显著性检验包括两方面
- 线性关系检验
- 回归系数检验
线性关系检验
- 提出假设: H 0 : β 1 = 0 {H}_0:\mathcal{\beta}_1=0 H0:β1=0 两个变量之间的线性关系不显著
- 计算检验统计量F:
F = E S S / 1 R S S / ( n − 2 ) = M S R M S E ∼ F ( 1 , n − 2 ) \mathrm{F}\:=\:\frac{\mathrm{ESS}/1}{\mathrm{RSS}/(\mathrm{n}-2)}\:=\:\frac{\mathrm{MSR}}{\mathrm{MSE}}\:\sim\:\mathrm{F}\:(1,\mathrm{n}-2) F=RSS/(n−2)ESS/1=MSEMSR∼F(1,n−2) - 确定显著性水平 α \alpha α
- 作出决策:
- 用
F
F
F 分布:查找临界值
F
α
(
1
,
n
−
2
)
{F} _\alpha ( 1,n-2)
Fα(1,n−2) 在F 分布表中的值
- F > F α {F}>{F}_\alpha F>Fα,拒绝 H 0 {H}_0 H0 ,表明两个变量之间的线性关系是显著。
- F < F a {F}<\mathcal{F}_a F<Fa,不拒绝 H 0 {H}_0 H0 ,没有证据表明两个变量之间的线性关系显著。
- 用
P
P
P 值(查表):
- 若 P < α {P}<\alpha P<α,拒绝 H 0 {H}_0 H0 ,表明两个变量之间的线性关系显著
- 若 P > α {P}>\alpha P>α ,不拒绝 H 0 {H}_0 H0,没有证据表明两个变量之间的线性关系显著。
- 用
F
F
F 分布:查找临界值
F
α
(
1
,
n
−
2
)
{F} _\alpha ( 1,n-2)
Fα(1,n−2) 在F 分布表中的值
回归系数检验
S
β
^
1
=
S
e
∑
x
i
2
−
1
n
(
∑
x
i
)
2
S
e
=
∑
(
y
i
−
y
^
i
)
2
n
−
k
−
1
=
M
S
E
S_{\widehat{\beta}_{1}}=\frac{S_{e}}{\sqrt{\sum x_{i}^{2}-\frac{1}{n}(\sum x_{i})^{2}}}\quad S_{e}=\sqrt{\frac{\sum(y_{i}-\hat{y}_{i})^{2}}{n-k-1}}=\sqrt{MSE}
Sβ
1=∑xi2−n1(∑xi)2SeSe=n−k−1∑(yi−y^i)2=MSE
-
提出假设:
- 两个变量之间的线性关系不显著 H 0 : β 1 = 0 {H}_0:\beta_1=0 H0:β1=0
- 两个变量之间的线性关系显著 H 1 : β 1 ≠ 0 {H}_1:\beta_1\neq0 H1:β1=0
-
计算检验统计量 t t t :
t = β ^ 1 s β ^ 1 ∼ t ( n − 2 ) t=\frac{\widehat{\beta}_{1}}{s_{\widehat{\beta}_{1}}}\sim t(n-2) t=sβ 1β 1∼t(n−2) -
确定显著性水平 a a a
-
作出决策:
- 用
F
F
F 分布:查找临界值
t
α
/
2
(
n
−
2
)
t_{\alpha / 2}( n- 2)
tα/2(n−2) 在
F
{F}
F 分布表中的
- t > t α / 2 \mathrm{~t>t}_{\alpha/2} t>tα/2,拒绝 H 0 {H}_0 H0 ,回归系数等于 0 的可能性小于 α,表明两个变量之间的线性关系是显著的。
- t < t α / 2 t< \mathrm{t}_{\alpha/2} t<tα/2,不拒绝 H 0 {H}_0 H0 ,没有证据表明两个变量之间的线性关系显著。
- 用
P
P
P 值:
- 若 P < α {P}<\alpha P<α ,拒绝 H 0 {H}_0 H0,表明两个变量之间的线性关系是显著的
- 若 P > α {P}>\alpha P>α,不拒绝 H 0 {H}_0 H0,二者不存在显著的线性关系。
- 用
F
F
F 分布:查找临界值
t
α
/
2
(
n
−
2
)
t_{\alpha / 2}( n- 2)
tα/2(n−2) 在
F
{F}
F 分布表中的
二、示例
某团队测了16名成年女子的身高与腿长所得数据如下
分析身高和腿长的关系
三、代码实现----Matlab
% 1、输入数据
% 输入X的样本值
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
% 插入\beta0对应列
X=[ones(16,1) x];
% 输入Y的样本值
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
% 2、回归分析及检验:
[b,bint,r,rint,stats]=regress(Y,X);
% 输出我们需要的数据
b,bint,stats,rint,r
% 结果:
% b =
% -16.0730
% 0.7194
% bint =
% -33.7071 1.5612
% 0.6047 0.8340
% stats =
% 0.9282 180.9531 0.0000 1.7437
% 即β_0=-16.073 ,β_1=0.7194;
% β_0 的置信区间为[-33.7017 ,1.5612],β_1 的置信区间为[0.6047 ,0.834];
% r_2=0.9282,F=180.9531,p=0.0000
% 由p<0.05,可知回归模型 y=-16.073+0.7194x 成立
% 3、残差分析,做残差图
rcoplot(r,rint)
% 从残差图可以看出,除第二个数据外,其余数据的残差离零点均较近,且残差的置信区间均包含零点,这说明回归模型
% y=-16.073+0.7194x能较好的符合原始数据,而第二个数据可视为异常点
% 4、预测及作图
z=b(1)+b(2)*x
plot(x,Y,'k+',x,z,'r')
运行结果:
残差图:
数据拟合图:
从运行结果可以发现,第二个数据异常。数据拟合程度较好。
四、代码实现----python
回归系数的置信区间公式
假设我们有一个简单的线性回归模型:
Y
=
β
X
+
ϵ
Y=\beta X + \epsilon
Y=βX+ϵ
其中:
- Y Y Y 是响应变量(因变量)向量。
- X X X 是设计矩阵(自变量)矩阵。
- β \beta β 是回归系数向量。
- ϵ \epsilon ϵ 是误差项,假设其独立同分布,且符合正态分布 ϵ ∼ N ( 0 , σ 2 ) \epsilon\sim N(0,\sigma^{2}) ϵ∼N(0,σ2)。
回归系数的置信区间可以通过以下步骤计算:
- 估计回归系数
β
^
:
\hat{\beta}:
β^:
β ^ = ( X T X ) − 1 X T Y \hat{\beta}=(X^TX)^{-1}X^TY β^=(XTX)−1XTY - 计算残差向量:
r = Y − X β ^ r=Y-X\hat{\beta} r=Y−Xβ^ - 计算均方误差(Mean Squared Error, MSE):
M S E = r T r n − p \mathrm{MSE}=\frac{r^Tr}{n-p} MSE=n−prTr
其中 n n n是观测值的数量, p p p是回归系数的数量(包括截距项). - 计算回归系数估计的标准误差
S E ( β ^ ) = M S E ⋅ d i a g ( ( X T X ) − 1 ) \mathrm{SE}(\hat{\beta})=\sqrt{\mathrm{MSE}\cdot\mathrm{diag}((X^TX)^{-1})} SE(β^)=MSE⋅diag((XTX)−1)
这里 diag ( ( X T X ) − 1 ) \operatorname{diag}((X^TX)^{-1}) diag((XTX)−1) 表示矩阵 ( X T X ) − 1 (X^TX)^{-1} (XTX)−1 的对角元素。 - 确定
t
t
t 分布的临界值
t
α
/
2
,
n
−
p
t_{\alpha/2,n-p}
tα/2,n−p :
根据所需的置信水平(例如 95%),从t分布表或使用统计函数(如 stats.t.ppf) 找到对应的临界值。
残差的置信区间公式
回归模型的残差向量
r
r
r 表示为:
r
=
Y
−
Y
^
=
Y
−
X
β
^
r=Y-\hat{Y}=Y-X\hat{\beta}
r=Y−Y^=Y−Xβ^
残差的标准误差
残差的标准误差(Residual Standard Error, RSE) 可以通过以下公式计算:
R
S
E
=
∑
r
i
2
n
−
p
\mathrm{RSE}=\sqrt{\frac{\sum r_i^2}{n-p}}
RSE=n−p∑ri2
其中:
- n n n 是观测值的数量。
- p p p 是回归系数的数量(包括截距项)。
残差的置信区间
残差的置信区间基于 t t t 分布,这意味看我们需要找到一个 t t t 统计量,该统计量依赖于给定的置信水平(例如 95% 置信区间) 和自由度(通常为 n − p n-p n−p) 。
残差
r
i
r_i
ri 的置信区间可以表示为
r
i
±
t
α
/
2
,
n
−
p
⋅
R
S
E
r_i\pm t_{\alpha/2,n-p}\cdot\mathrm{RSE}
ri±tα/2,n−p⋅RSE
其中:
t α / 2 , n − p t_{\alpha/2,n-p} tα/2,n−p 是t分布在自由度 n − p n-p n−p 和显著水平 α \alpha α 下的临界值。
python代码
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] #设置字体
mpl.rcParams['axes.unicode_minus'] = False # - 号设置
# 1、输入数据
# 输入X的样本值
x = np.array([143, 145, 146, 147, 149, 150, 153, 154, 155, 156, 157, 158, 159, 160, 162, 164])
# 插入\beta0对应列
X = sm.add_constant(x)
# 输入Y的样本值
Y = np.array([88, 85, 88, 91, 92, 93, 93, 95, 96, 98, 97, 96, 98, 99, 100, 102])
# 2、回归分析及检验
model = sm.OLS(Y, X)
results = model.fit()
# 输出我们需要的数据
print("回归系数 (b):")
print(results.params)
print("\n回归系数的区间估计 (bint):")
print(results.conf_int())
print("\n模型统计信息 (stats):")
print("R-squared:", results.rsquared)
print("F-statistic:", results.fvalue)
print("p-value:", results.f_pvalue)
print("残差方差 (RMSE):", results.mse_resid)
# 5. 预测
y_pred = results.predict(X)
# 6. 计算残差
residuals = Y - y_pred
# 获取回归系数
b = results.params
bint = results.conf_int()
# 计算残差的置信区间
alpha = 0.05 # 95% 置信区间
n = len(residuals)
df = n - len(b) # 自由度
t_value = stats.t.ppf(1 - alpha / 2, df) # 找到 t 分布的分位点
# 计算标准误差
s_err = np.sqrt(np.sum(residuals**2) / df)
# 置信区间
rint = np.zeros((n, 2))
for i in range(n):
rint[i, 0] = residuals[i] - t_value * s_err
rint[i, 1] = residuals[i] + t_value * s_err
# 7. 绘制残差图
plt.figure(figsize=(10, 6))
plt.scatter(range(n), residuals, color='blue', marker='o', label='Residuals')
plt.hlines(0, xmin=0, xmax=n, color='red', linewidth=2)
plt.plot(range(n), rint[:, 0], 'r_', label='Lower bound')
plt.plot(range(n), rint[:, 1], 'r_', label='Upper bound')
plt.xlabel('Observation')
plt.ylabel('Residuals')
plt.title('Residuals with Confidence Interval')
plt.legend()
plt.show()
# 4、预测及作图
z = results.predict(X)
plt.scatter(x, Y, color='black', marker='+', label='实际数据')
plt.plot(x, z, color='red', label='拟合线')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
运行结果:
总结
本文介绍了一元线性回归,详细介绍了 回归分析的一般步骤。分别使用matlab和python进行了算法的实现。