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

---

前言

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

一、多元线性回归

多元线性回归模型

一般称由 $y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_k{x}_k$ 确定的模型：
$$\{\begin{array}{c}Y=X\beta +ϵ\\ E(ϵ)=0,COV(ϵ,ϵ)={\sigma }^2{I}_n\end{array}$$
为 $k$ 元线性回归模型，并简记为 $(Y,X\beta ,{\sigma }^2{I}_n)$。
$$\begin{array}{r}Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right],X=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& ...& x_{1k}\\ 1& x_{21}& x_{22}& ...& x_{2k}\\ ...& ...& ...& ...& ...\\ 1& x_{n1}& x_{n2}& ...& x_{nk}\end{array}\right],\beta =\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ...\\ {\beta }_k\end{array}\right],ϵ=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ...\\ {ϵ}_n\end{array}\right]\end{array}$$
$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_k{x}_k$称为回归平面方程。

线性模型 $(Y,X\beta ,{\sigma }^2{I}_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 }_1{x}_{i1}-\cdots -{\beta }_k{x}_{ik}{)}^2$$
• 选择 ${\beta }_0,\cdots ,{\beta }_k$ 使 $Q$ 达到最小
• 解得估计值$\hat{\beta }=(X^{T}X{)}^{-}1(X^{T}Y)$
• 得到的${\hat{\beta }}_i$代入回归平面方程得：
  $$y={\hat{\beta }}_0+{\hat{\beta }}_1{x}_1+\cdots +{\hat{\beta }}_k{x}_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,\cdot \cdot \cdot ,x_k$ 之间显著地有线性关系；否则就接受$H_0$，认为 $y$ 与 $x_1,\cdots ,x_k$ 之间线性关系不显著。
• 其中$$\begin{gathered} U=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2(\text{回归平方和}) \\ {Q}_e=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2(\text{残差平方和}) \end{gathered}$$

$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=\left[\begin{array}{c}{\hat{\beta }}_0\\ {\hat{\beta }}_1\\ ...\\ {\hat{\beta }}_p\end{array}\right],Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right],X=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1p}\\ 1& x_{21}& x_{22}& \cdots & x_{2p}\\ ...& ...& ...& ...& ...\\ 1& x_{n1}& x_{n2}& ...& x_{np}\end{array}\right]$$

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');

运行结果：

残差分析：

三维图：