【学习笔记】Matlab和python双语言的学习(一元线性回归）

---

前言

通过模型算法，熟练对Matlab和python的应用。
学习视频链接：
https://www.bilibili.com/video/BV1EK41187QF?p=42&vd_source=67471d3a1b4f517b7a7964093e62f7e6

一、一元线性回归

• 一元线性回归是一种统计方法，用于理解和建模两个变量之间的关系。具体来说，它是用来估计一个因变量 $y$ 和一个自变量 $x$ 之间线性关系的工具。其目标是找到一条直线，能够最好地预测 $y$ 的值，基于给定的 $x$ 值。这条直线通常被称为 “回归线” 。

回归分析的一般步骤

一元线性回归的基本形式

$$y={\beta }_0+{\beta }_{1}x+ϵ$$

• $y$ 是因变量（被解释变量）。
• $x$ 是自变量（解释变量）。
• ${\beta }_0$ 是截距（当 $x=0$ 时， $y$ 的预测值）。
• ${\beta }_1$ 是斜率（表示 $x$ 每增加一个单位，$y$ 变化的量）。
• $ϵ$ 是误差项，表示因变量的实际值和预测值之间的差异。

回归方程

• 描述因变量y的期望值如何依赖于自变量x的方程称为回归方程。根据对一元线性回归模型的假设，可以得到它的回归方程为:
  $$E(y)={\beta }_0+{\beta }_{1}x$$
• 如果回归方程中的参数已知，对于一个给定的x值，利用回归方程就能计算出y的期望值
• 用样本统计量代替回归方程中的未知参数，就得到估计的回归方程，简称回归直线

参数的最小二乘法估计

• 对于回归直线，关键在于求解参数，常用高斯提出的最小二乘法，它是使因变量的观察值 $y$ 与估计值之间的离差平方和达到最小来求解的。
  $$Q=\sum (y-\hat{y}{)}^2=\sum (y-{\hat{\beta }}_0-{\hat{\beta }}_{1}x{)}^2$$
  展开可得：
  $$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$$
  求偏导可得：
  $$\{\begin{array}{l}\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{array}$$

$$\{\begin{array}{l}{\hat{\beta }}_1=\frac{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{array}$$

即可得到${\hat{\beta }}_0$及${\hat{\beta }}_1$

对回归方程的各种检验

估计标准误差的计算

• 为了度量回归方程的可靠性，通常计算估计标准误差，它度量观察值回绕着回归直线的变化程度或分散程度
• 估计平均误差：
  $$S_e=\sqrt{\frac{\Sigma (y-\hat{y}{)}^2}{n-2}}$$
• 估计标准误差越大，则数据点围绕回归直线的分散程度就越大，回归方程代表性就越小
• 估计标准误差越小，则数据点围绕回归直线的分散程度越小，回归方程的代表性越大，可靠性越高

回归直线的拟合优度

• 回归直线与各观测点的接近程度称为回归直线对数据的拟合优度
  • 总平方和 (TSS)：反映因变量的$n$个观察值与其均值的总离差 $$TSS=\sum y_i^2=\sum (y_i-{\bar{y}}_i{)}^2$$
  • 回归平方和 (ESS)：反映了 $y$ 的总变差中，由于 $x$ 与 $y$ 之间的线性关系引起的 $y$ 的变化部分$$ESS=\sum {\hat{y}}_i^2=\sum ({\hat{y}}_i-{\bar{y}}_i{)}^2$$
  • 残差平方和(RSS)：反映了除了$x$对$y$的线性影响之外的其他因素对$y$变差的作用，是不能由回归直线来解释的y的变差部分$$RSS=\sum e_i^2=\sum (y_i-{\hat{y}}_i{)}^2$$
  • 总平方和可以分解为回归平方和、残差平方和两部分：$TSS=ESS+RSS$

判定系数

• 回归平方和占总平方和的比例，用 $R^2$ 表示，其值在0到1之间
  • $R^2=0$：说明 $y$ 的变化与 $x$ 无关，$x$ 完全无助于解释 $y$ 的变差
  • $R^2=1$：说明残差平方和为0，拟合是完全的，$y$ 的变化只与 $x$ 有关

显著性检验

• 显著性检验的主要目的是根据所建立的估计方程用自变量 $x$ 来估计或预测因变量 $y$ 的取值，当建立了估计方程后，还不能马上进行估计或预测，因为该估计方程是根据样本数据得到的，它是否真实的反映了变量 $x$ 和 $y$ 之间的关系，则需要通过检验后才能证实。
• 根据样本数据拟合回归方程时，实际上就已经假定变量 $x$ 和 $\nu$ 之间存在着线性关系，并假定误差项是一个服从正态分布的随机变量，且具有相同的方差，但这些假设是否成立需要检验
• 显著性检验包括两方面
  • 线性关系检验
  • 回归系数检验

线性关系检验

1. 提出假设：$H_0:{\beta }_1=0$ 两个变量之间的线性关系不显著
2. 计算检验统计量F：
   $$\mathrm{F}=\frac{\mathrm{ESS}/1}{\mathrm{RSS}/(\mathrm{n}-2)}=\frac{\mathrm{MSR}}{\mathrm{MSE}}\sim \mathrm{F}(1,\mathrm{n}-2)$$
3. 确定显著性水平$\alpha$
4. 作出决策：
   • 用 $F$ 分布：查找临界值 $F_{\alpha }(1,n-2)$ 在F 分布表中的值
     • $F>F_{\alpha }$，拒绝 $H_0$ ，表明两个变量之间的线性关系是显著。
     • $F<{\mathcal{F}}_a$，不拒绝 $H_0$ ，没有证据表明两个变量之间的线性关系显著。
   • 用 $P$ 值（查表）：
     • 若$P<\alpha$，拒绝 $H_0$ ，表明两个变量之间的线性关系显著
     • 若$P>\alpha$ ，不拒绝 $H_0$，没有证据表明两个变量之间的线性关系显著。

回归系数检验
$$S_{{\hat{\beta }}_1}=\frac{S_e}{\sqrt{\sum x_i^2-\frac{1}{n}(\sum x_i{)}^2}}{S}_e=\sqrt{\frac{\sum (y_i-{\hat{y}}_i{)}^2}{n-k-1}}=\sqrt{MSE}$$

1. 提出假设：

   • 两个变量之间的线性关系不显著 $H_0：{\beta }_1=0$
   • 两个变量之间的线性关系显著 $H_1：{\beta }_1\ne 0$

2. 计算检验统计量 $t$ ：
   $$t=\frac{{\hat{\beta }}_1}{s_{{\hat{\beta }}_1}}\sim t(n-2)$$

3. 确定显著性水平$a$

4. 作出决策：

   • 用 $F$ 分布：查找临界值 $t_{\alpha /2}(n-2)$ 在 $F$ 分布表中的
     • ${\mathrm{t}>\mathrm{t}}_{\alpha /2}$，拒绝 $H_0$ ，回归系数等于 0 的可能性小于 α，表明两个变量之间的线性关系是显著的。
     • $t<{\mathrm{t}}_{\alpha /2}$，不拒绝 $H_0$ ，没有证据表明两个变量之间的线性关系显著。
   • 用 $P$ 值：
     • 若 $P<\alpha$ ，拒绝 $H_0$，表明两个变量之间的线性关系是显著的
     • 若 $P>\alpha$，不拒绝 $H_0$，二者不存在显著的线性关系。

二、示例

某团队测了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=\beta X+ϵ$$
其中：

• $Y$ 是响应变量（因变量）向量。
• $X$ 是设计矩阵（自变量）矩阵。
• $\beta$ 是回归系数向量。
• $ϵ$ 是误差项，假设其独立同分布，且符合正态分布$ϵ\sim N(0,{\sigma }^2)$。

回归系数的置信区间可以通过以下步骤计算：

1. 估计回归系数 $\hat{\beta }:$
   $$\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$$
2. 计算残差向量：
   $$r=Y-X\hat{\beta }$$
3. 计算均方误差(Mean Squared Error, MSE):
   $$\mathrm{MSE}=\frac{r^{T}r}{n-p}$$
   其中$n$是观测值的数量，$p$是回归系数的数量(包括截距项).
4. 计算回归系数估计的标准误差
   $$\mathrm{SE}(\hat{\beta })=\sqrt{\mathrm{MSE}\cdot \mathrm{diag}((X^{T}X{)}^{-1})}$$
   这里 $\mathrm{diag}((X^{T}X{)}^{-1})$ 表示矩阵 $(X^{T}X{)}^{-1}$ 的对角元素。
5. 确定 $t$ 分布的临界值 $t_{\alpha /2,n-p}$ :
   根据所需的置信水平(例如 95%)，从t分布表或使用统计函数(如 stats.t.ppf) 找到对应的临界值。

残差的置信区间公式

回归模型的残差向量 $r$ 表示为：
$$r=Y-\hat{Y}=Y-X\hat{\beta }$$
残差的标准误差

残差的标准误差(Residual Standard Error, RSE) 可以通过以下公式计算：
$$\mathrm{RSE}=\sqrt{\frac{\sum r_i^2}{n-p}}$$

其中：

• $n$ 是观测值的数量。
• $p$ 是回归系数的数量(包括截距项)。

残差的置信区间

残差的置信区间基于 $t$ 分布，这意味看我们需要找到一个 $t$ 统计量，该统计量依赖于给定的置信水平(例如 95% 置信区间) 和自由度(通常为$n-p$) 。

残差 $r_i$ 的置信区间可以表示为$$r_i\pm t_{\alpha /2,n-p}\cdot \mathrm{RSE}$$
其中：

$t_{\alpha /2,n-p}$ 是t分布在自由度 $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进行了算法的实现。