【学习笔记】多元线性回归模型 —— Matlab

【学习笔记】多元线性回归模型 —— Matlab

引用

CSDN

1.

https://blog.csdn.net/m0_65032457/article/details/141332893

本文将介绍多元线性回归模型的理论基础、参数估计方法、模型检验方法以及在Matlab中的实现。通过一个水泥凝固时放热的示例，详细展示了如何使用逐步回归法确定最佳回归模型，并给出了相应的Matlab代码和结果分析。

前言

通过模型算法，熟练对Matlab的应用。

学习视频链接：
https://www.bilibili.com/video/BV1EK41187QF?p=44&vd_source=67471d3a1b4f517b7a7964093e62f7e6

一、多元线性回归

多元线性回归模型

一般称由 $y = \beta_0 + \beta_1x_1 + \cdots + \beta_kx_k$ 确定的模型：

$$
\left{
\begin{matrix}
Y=X\beta+\epsilon \
E(\epsilon)=0, COV(\epsilon, \epsilon)=\sigma^2I_n
\end{matrix}
\right.
$$

为 $k$ 元线性回归模型，并简记为 $(Y, X\beta, \sigma^2I_n)$。

$$
Y=\begin{bmatrix}
y_1 \
y_2 \
... \
y_n
\end{bmatrix}, X=\begin{bmatrix}
1 & x_{11} & x_{12} & ... & x_{1k} \
1 & x_{21} & x_{22} & ... & x_{2k} \
... & ... & ... & ... & ... \
1 & x_{n1} & x_{n2} & ... & x_{nk}
\end{bmatrix}, \beta=\begin{bmatrix}
\beta_0 \
\beta_1 \
... \
\beta_k
\end{bmatrix}, \epsilon=\begin{bmatrix}
\epsilon_1 \
\epsilon_2 \
... \
\epsilon_n
\end{bmatrix}
$$

$y = \beta_0 + \beta_1x_1 + \cdots + \beta_kx_k$ 称为回归平面方程。

线性模型 $(Y, X\beta, \sigma^2I_n)$ 考虑的主要问题

• 对参数 $\beta$ 和 $\sigma^2$ 作点估计，建立 $y$ 与 $x_1, x_2, \cdots, x_k$ 之间的数量关系
• 对模型参数、模型结果等做检验
• 对 $\gamma$ 的值作预测，即对 $\gamma$ 作点(区间)估计

多元线性回归模型的参数估计

• 用最小二乘法求 $\beta_0, \cdots, \beta_k$ 的估计量：作离差平方和
  $$
  Q=\sum_{i=1}^n(y_i-\beta_0-\beta_1x_{i1}-\cdots-\beta_kx_{ik})^2
  $$
• 选择 $\beta_0, \cdots, \beta_k$ 使 $Q$ 达到最小
• 解得估计值 $\hat{\beta}=(X^TX)^{-1}(X^TY)$
• 得到的 $\hat{\beta}_i$ 代入回归平面方程得：
  $$
  y=\hat{\beta}_0+\hat{\beta}_1x_1+\cdots+\hat{\beta}_kx_k
  $$
  称为经验回归平面方程，$\hat{\beta}_i$ 称为经验回归系数

多元线性回归模型和回归系数的检验

1) F 检验法

• 当 $H_0$ 成立时，$F=\frac{U/k}{Q_e/(n-k-1)} \sim F(k, n-k-1)$
• 如果 $F > F_{1-\alpha}(k, n-k-1)$，则拒绝 $H_0$，认为 $y$ 与 $x_1, \cdots, x_k$ 之间显著地有线性关系；否则就接受 $H_0$，认为 $y$ 与 $x_1, \cdots, x_k$ 之间线性关系不显著。
• 其中 $U=\sum_{i=1}^{n}(\hat{y}i-\bar{y})^2$ (回归平方和)
  $Q_e=\sum{i=1}^{n}(y_i-\hat{y}_i)^2$ (残差平方和)

2) r 检验法

• 定义
  $$
  R=\sqrt{\frac{U}{L_{yy}}}=\sqrt{\frac{U}{U+Q_e}}
  $$
  称为 $y$ 与 $x_1, \cdots, x_k$ 的多元相关系数或复相关系数。由于
  $$
  F=\frac{n-k-1}{k}\frac{R^2}{1-R^2}
  $$
  故用 $F$ 和用 $R$ 检验是等效的。

二、示例

• 水泥凝固时放出的热量 $y$ 与水泥中4种化学成分 $x_1, x_2, x_3, x_4$ 有关，今测得一组数据如下，试用逐步回归法确定一个线性模型，并找出影响水泥凝固时放出热量的必要因素

多元线性回归的逐步回归

• “最优” 的回归方程就是包含所有对 $Y$ 有影响的变量，而不包含对 $Y$ 影响不显著的变量回归方
• 逐步回归分析法的思想：
• 从一个自变量开始，根据自变量对 $Y$ 作用的显著程度，从大到小地依次逐个引入回归方程
• 当引入的自变量由于后引入变量而变得不显著时，要将其剔除掉。
• 引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步。
• 对于每一步都要进行 $Y$ 值检验，以确保每次引入新的显著性变量前回归方程中只包含对 $Y$ 作用显著的变量。
• 这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程时为止。

三、代码实现——Matlab

1. 多元回归的实现

Matlab 的 regress 函数

regress 函数的基本语法如下：

[b,bint,r,rint,stats]=regress(Y,X,alpha)

$b = \begin{bmatrix}
\hat{\beta}0 \
\hat{\beta}1 \
... \
\hat{\beta}p
\end{bmatrix}, Y=\begin{bmatrix}
y_1 \
y_2 \
... \
y_n
\end{bmatrix}, X=\begin{bmatrix}
1 & x{11} & x{12} & \cdots & x{1p} \
1 & x_{21} & x_{22} & \cdots & x_{2p} \
... & ... & ... & ... & ... \
1 & x_{n1} & x_{n2} & ... & x_{np}
\end{bmatrix}$

2. 逐步回归的实现

Matlab 的 stepwise 函数

stepwise 函数的基本语法如下：

stepwise(x, y, inmodel, alpha)

• $x$ — 自变量数据，$n \times m$ 阶矩阵
• $y$ — 因变量数据，$n \times 1$ 阶矩阵
• inmodel — 矩阵的列数的指标，给出初始模型中包括的子集(缺省时设定为全部自变量)
• alpha — 默认为 0.05

3. Matlab 代码

clc;clear
x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';
x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';
x3=[6 15 8 8 6 9 17 22 18 4 23 9 8]';
x4=[60 52 20 47 33 22 6 44 22 26 34 12 12]';
y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4]';
x=[x1 x2 x3 x4];  % 自变量矩阵，包括四个自变量
stepwise(x,y)  % 使用逐步回归分析方法，选取最佳模型

X1、X2、X3、X4 全部选中时，F = 111.479

只选中 X1、X2 时，F = 229.504

只选中 X3、X4 时，F = 72.2674

可见 X1、X2 是影响 Y 的必要因素

只对 X1、X2 进行线性回归：

% 水泥凝固时放热分析
% （1）数据输入
clc;clear
x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';
x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';
y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4]';
x=[ones(size(x1)),x1, x2];  % 自变量矩阵，包括四个自变量
% （2）求结果
[b,bint,r,rint,stats]=regress(y,x)  % 进行多元线性回归分析
% （3）画残差图
rcoplot(r,rint)  % 绘制残差图，用于评估回归模型的拟合情况
% 3、作图
% 定义自变量的范围
x1 = -10:0.1:10;  % x1的取值范围
x2 = -10:0.1:10;  % x2的取值范围
% 计算对应的因变量值
[X1, X2] = meshgrid(x1, x2);  % 创建网格点坐标矩阵
Z = b(1) + b(2)*X1 + b(3)*X2;  % 计算对应的因变量值
% 绘制三维图形
figure;
surf(X1, X2, Z);
xlabel('x1');
ylabel('x2');
zlabel('z');
title('Z = 52.5773 + 1.4683*X1 + 0.6623*X2');

运行结果：

残差分析：

三维图：

本文原文来自CSDN