回归

用于在难以选择建模模型时，分析自变量和因变量的关系

一、多元线性回归

（一）原理理论

1. 模型

多元线性回归分析的模型为

其中$\beta_0,\beta_1,\cdots,\beta_m,\sigma^2$都是与$x_1,x_2,\cdots,x_m$无关的位置参数，其中$\beta_0,\beta_1,\cdots,\beta_m,$称为回归系数

现得到$n$个独立观测数据$[bi,a{\small i1},\cdots,a{i2}]$，其中$bi$为$y$的观察值，$a{\small i1},\cdots,a{im}$分别是$x_1,x_2,\cdots,x_m$的观测值，$i=1,\cdots,n\quad,n>m$，即

记

则上式可表示为

其中$E_n$为$n$阶单位矩阵

2. 参数估计

模型式中的参数$\beta_0,\beta_1,\cdots,\beta_m$用最小二乘法进行估计，即应选取估计值$\hat{\beta}{\small j}$，使当$\beta_j=\hat{\beta}{\small j}\quad j=0,1,\cdots,m$时，误差平方和

达到最小，为此，令

得

经整理化为以下正规方程组

正规方程组的矩阵形式为：$X^{\rm T}X\beta=X^{\rm T}Y$

当矩阵$X$列满秩时，$X^{\rm T}X$为可逆方阵，$X^{\rm T}X\beta=X^{\rm T}Y$的解为：

将$\hat{\beta}$代回原模型得到y的估计值

而这组数据的拟合值为

记$\hat{Y}=X\hat{\beta}=[\hat{b}_1,\cdots,\hat{b}{\small n}]^{\rm T}$，拟合误差$e=Y-\hat{Y}$称为残差，可作为随机误差$\varepsilon$的估计，

而

为残差平方和（或剩余平方和）

3. 统计分析

不加证明的给出以下结果

（1）$\hat{\beta}$是$\beta$的线性无偏最小方差估计，$\hat{\beta}$的期望=$\beta$，在$\beta$的线性无偏估计中,$\hat{\beta}$的方差最小

（2）$\hat{\beta}$服从正态分布： $\hat{\beta}\sim N(\beta,\sigma^2(X^{\rm T}X)^{-1})$ ,记$(X^{\rm T}X)^{-1} = (c{\small ij})_{n\times n}$

（3）对残差平方和$Q$，$EQ=(n-m-1)\sigma^2$，且 $\dfrac{Q}{\sigma^2}\sim \chi^2(n-m-1)$，由此得到$\sigma^2$的无偏估计

其中$s^2$是剩余方差（残差的方差）,$s$称为剩余标准差

（4）对总平方和$SST$进行分解，有

其中，

为残差平方和，反应随机误差对$y$的影响；

称为回归平方和，反应自变量对$y$的影响。上述的分解利用了正规方程组。

4. 回归模型的假设检验

因变量$y$和自变量$x_1,\cdots,x_m$之间是否存在模型所示的线性关系是需要检验的

显然，如果所有的$|\hat{\beta}{\small j}|\ (j=1,\cdots,m)$都很小，则$y$与$x_1,\cdots,x_m$的线性关系就不明显了，所以可以令原假设为：

当$H_0$成立时由分解式$SST=Q+U$定义的$Q、U$满足

在显著性水平$\alpha$下，对于上$\alpha$分位数$F{\small \alpha}(m,n-m-1)$，若$F<F{\small \alpha}(m,n-m-1)$，接受$H_0$，否则拒绝

接受$H_0$只能说明$y$与$x_1,\cdots,x_m$的线性关系不明显，可能存在非线性关系，如：平方关系

其他用于衡量$y$与$x_1,\cdots,x_m$相关程度的指标：回归平方和在总平方和中的比值定义复判定系数

其中$R=\sqrt{R^2}$称为复相关系数，$R$越大，$y$与$x_1,\cdots,x_m$相关关系越密切。通常$R>0.8$或$R>0.9$才认为相关关系成立

5. 利用回归模型进行预测

当回归模型和系数通过检验后，可由给定$[x1,\cdots,x_m]$的取值$[a{\small 01},\cdots,a{0m}]$预测$y$的取值$b_0$，$b_0$是随机的，显然其预测值（点估计）为

给定$\alpha$可以算出$b_0$的预测区间(区间估计)，结果较复杂，但当$n$较大且$a{\small 0i}$接近平均值$\barx_i$时，$b_0$的预测区间可简化为

其中$z{\small {\frac{\alpha}{2}}}$是标准正态分布的上$\dfrac{\alpha}{2}$分位数

对$b_0$的区间估计方法可用于给出已知数据残差$e_i=b_i-\hat{b}_i\ (i=1,\cdots,n)$的置信区间，$e_i$服从均值为0的正态分布（$e_i\sim N(0,\sigma^2)$），所以若某个$e_i$的置信区间不包含零点，则认为这个数据是异常的，可予以剔除

（二）MATLAB中的多元线性回归分析

MATLAB统计工具箱中有一个实现多元线性回归的函数regress

Code : b = regress(Y, X)

其中，Y和X为按(22)式排列的数据，b为回归系数估计值

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

其中，Y和X为按(22)式排列的数据，alpha为显著性水平（缺省值为0.05），b和bint为回归系数估计值和它们的置信区间，r和rint为残差(向量)及其置信区间

stats是用于检验回归模型的统计量，有四个数值：

第一个是$R^2$；第二个是$F$；第三个是与$F$对应的概率$p\ \ ,p<\alpha$拒绝$H_0$，回归模型成立；第四个是残差的方差$s^2$

残差及其置信区间可以用rcoplot(r, rint)画图

二、多元二项式回归

MATLAB统计工具箱提供了一个作多元二项式回归的命令rstool，它能够产生一个交互式画面，并输出有关信息

Code : rstool(X, Y, medol, alpha)

其中alpha为显著性水平 $\alpha$（缺省值为0.05），model可以选择如下四种模型，缺省值为'linear'线性模型

image-20240216012623970

$[y,x1,\cdots,x_m]$的$n$个独立观测数据仍然记为。。$[b_i,a{\small i1},\cdots,a{im}]\quad i=1,\cdots,n$

Y,XX分别为n维列向量和$n\times m$的矩阵

注意！

（1）此处的XX数据矩阵与线性回归分析中的数据矩阵X不同，X是构建的表达式系数矩阵，其第一列为全1列，XX是列向量观测值横向拼接成的矩阵

（2）在完全二次多项式回归中，二次项系数的排列次序是先交叉项的系数，最后是纯二次项的系数

在交互式图像界面可以通过键选或输入相应值来得出其预测值和置信区间

并且可以通过导出按钮导出显影的拟合系数到变量beta，RMSE到变量rmse，残差保存到变量residuals中

RMSE(Root Mean Squared Error)是均方根误差，是一种衡量预测值与实际值之间误差的方法，通常用于评估模型的预测性能，计算公式为：

外部实现代码为

1
2
3
4
5
6
y_actual = [actual values];
y_predicted = [predicted values];
% 计算残差
residuals = y_actual - y_predicted;
% 计算 RMSE
rmse = sqrt(mean(residuals.^2));

但是rstool方法内置了数据标准化、添加截距项等步骤，所以外部实现的rmse不能直接同rstool输出的rmse进行比对的

三、非线性回归

非线性回归是指因变量$y$对应系数$\beta_1,\cdots,\beta_m$（而不是自变量）是非线性的

MATLAB统计工具箱中提供了四个命令用于非线性回归问题

Command Description
nlinfit 计算回归系数
nlparci 计算回归系数的置信区间
nlpredci 计算预测值及其置信区间
nlintool 得到交互式画面

命令方法解析

  • 非线性回归：[beta, r, j] = nlinfit(xr, yr, <xfx>, beta0)
    • 其中beta、r和j分别承接回归系数beta，残差r和jacobian矩阵
    • xr、yr为x和y的原始数据，y为列向量，x为列向量横向拼接成的矩阵，
    • <xfx>是匿名函数构建的目标拟合公式，beta0为beta的初始值，若未给可以自由选择

  • 回归系数的置信区间：<ValName> = nlparci(beta, r, 'jacobian', j)

    • 这里使用最常见的jacobian矩阵方法进行拟合构建，所以Model选择’jacobian’
    • 其他的条件由nlinfit命令得到

  • 预测值及其置信区间半径：[Val_pred, delta] = nlpredci(<xfx>, Xt, beta, r, 'jacobian', j)

    • val_pred和delta承接预测值及其置信区间半径
    • xt为预测值所需的自变量，构建方法同xr
    • 这里同样使用最常见的jacobian矩阵方法进行拟合构建，所以Model选择’jacobian’

  • tool 交互式非线性拟合：nlintool(xr, yr, <xfx>, beta0)

    • 交互式方法可以直接列用交互式图像完成上面的部分，只需要导入初始数据，默认使用’jacobian’模式

    • 利用‘导出’功能，可以直接导出下面的参数：

      | 参数 | 参数CI | 预测 | 预测CI | RMSE | 残差 |
      | —— | ——— | ——- | ———- | —— | ————- |
      | beta | betaci | ypred | ypredci | rmse | residuals |

      注意，这里的ypred和ypredci是交互式界面键入所得的值，而不能直接将其与实际结果的预测输出

所以nlintool交互式画面更容易求解基础的问题和RMSE等没有确定函数的问题

nlinfit \ nlparci \ nlpredci更适合建立一个复杂问题的模型，也更方面对大量的数据进行预测

四、实战例题

Eg.s1 合金强度拟合的回归分析

合金强度y与其中的碳含量x有比较密切的关系，今从生产中收集了一批数据如下表，试先拟合一个函数y(x)，再用回归分析对它进行检验

Data:

x 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18
y 42.0 41.5 45.0 45.5 45.0 47.5 49.0 55.0 50.0 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
clear
clc
%% Data
x1 = [0.10:0.01:0.18]'; % 仅有一个变量，尝试一元线性模型
y1 = [42.0 41.5 45.0 45.5 45.0 47.5 49.0 55.0 50.0]';
plot(x1,y1,'+') % 大致预览数据
%% 拟合
a = [ones(size(x1)),x1];
[b,bint,r,rint,stats] = regress(y1,a);
b,bint % 查看b值是否在置信区间内
stats %
subplot(2,2,1)
rcoplot(r, rint) % 第8个数据出现问题，删除该数据
%% ReData
clear
x1 = [0.10:0.01:0.18]'; % 仅有一个变量，尝试一元线性模型
y1 = [42.0 41.5 45.0 45.5 45.0 47.5 49.0 55.0 50.0]';
x1(8)=[];
y1(8)=[];
%% Re拟合
a = [ones(size(x1)),x1];
[b,bint,r,rint,stats] = regress(y1,a);
b,bint % 查看b值是否在置信区间内
stats %
subplot(2,2,2)
rcoplot(r, rint) % 第8个数据出现问题(异常点)，删除该数据
%% Others
clc
% CalTest最小二乘计算相应值
b1 = a\y1
% 拟合图像
x = min(x1):0.001:max(x1);
y = b(1)+b(2)*x;
subplot(2,2,4)
plot(x,y,x1,y1,'r+')
% MSE
y1_fit = b(1)+b(2)*x1;
disp(['MSE: ' num2str(mse(y1-y1_fit))])

output

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
b =
   27.4722
  137.5000
bint =
   18.6851   36.2594
   75.7755  199.2245
stats =

    0.7985   27.7469    0.0012    4.0883
    
b =
   30.7820
  109.3985
bint =
   26.2805   35.2834
   76.9014  141.8955
stats =
    0.9188   67.8534    0.0002    0.8797

b1 =
   30.7820
  109.3985
MSE: 0.65977

Pic

Eg.s2 电器的销售量预测

某厂生产的一种电器的销售量y与竞争对手的价格x1和本厂的价格x2有关。下表是该商品在10个城市的销售记录。试根据这些数据建立y与x1和x2的关系式，对得到的模型和系数进行检验。若某市本厂产品售价160 (元)，竞争对手售价170(元)，预测商品在该市的销售量。

Data:

$x_1$元 120 140 190 130 155 175 125 145 180 150
$x_2$元 100 110 90 150 210 150 250 270 300 250
$y$个 102 100 120 77 46 93 26 69 65 85

分别用多元线性回归 和 多元二项式回归进行预测

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
clear,clc
%% Data:
x1 = [120 140 190 130 155 175 125 145 180 150]';
x2 = [100 110 90 150 210 150 250 270 300 250]';
y = [102 100 120 77 46 93 26 69 65 85]';
%% ShowRawData
subplot(1,3,1)
plot(x1,y,'r+')
subplot(1,3,2)
plot(x2,y,'r+')
%% FitRegressAnalysis 拟合回归分析
a = [ones(size(x1)), x1, x2];
[b,bint,r,rint,stats] = regress(y,a);
b,bint,stats
subplot(1,3,3)
rcoplot(r,rint)
% 如果计算RMSE，其值将不与rstool得出的结果相同
% 原因是没有进行rstool内置的数据标准化、添加截距项等
% 所以regress得出的方法不能跟rstool得出的RMSE直接比对
%% x1:170 x2:160 y?
y_pred_1716 = b(1) + b(2)*170 + b(3)*160
%% Rstool 多元二项式回归
rstool([x1,x2],y,'purequadratic')
rstool([x1,x2],y,'linear')
% beta
% rmse

output

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
b =
   66.5176
    0.4139
   -0.2698
bint =
  -32.5060  165.5411
   -0.2018    1.0296
   -0.4611   -0.0785
stats =
    0.6527    6.5786    0.0247  351.0445

y_pred =
   93.7182
   
--------
beta =

 -312.5871
    7.2701
   -1.7337
   -0.0228
    0.0037
rmse =
   16.6436

Pic

Eg.s3 运动员问题

将17至29岁的运动员每两岁一组分为7组,每组两人测量其旋转定向能力,以考察年龄对这种运动能力的影响。现得到一-组数据如下，试建立两者之间的关系

Data:

Age 17 19 21 23 25 27 29
Person1 20.48 25.13 26.15 30.0 26.1 20.3 19.35
Person2 24.35 28.11 26.3 31.4 26.92 25.7 21.3

注意当前问题，虽然看似是含有两个自变量的多元问题，实际上是一个Age对应两个Person的值，所以在此处尝试使用多项式拟合poly的方法进行建模

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
clc,clear
%% Data
age = [17:2:29];
ages = [age,age];
h1 = [20.48 25.13 26.15 30.0 26.1 20.3 19.35];
h2 = [24.35 28.11 26.3 31.4 26.92 25.7 21.3];
h0 = [h1,h2];
%% poly
[p, s] = polyfit(ages, h0, 2);
p
[y, delta] = polyconf(p,age,s); % 带入模型进行预测
y % 拟合值
delta; % 预测区间的半径

polytool(ages, h0, 2) %  预测模型交互式图像

output

1
2
3
4
5
p =
   -0.2003    8.9782  -72.2150

y =
   22.5243   26.0582   27.9896   28.3186   27.0450   24.1689   19.6904

Pic

Eg7.20 化学动力反应模型

在研究化学动力学反应过程中，建立了一个反应速度和反应物含量的数学模型，形式为

其中$\beta_1,\cdots,\beta_5$是位置的参数，$x_1,x_2,x_3$是三种反应物(氢、n戊烷，异构戊烷)的含量，$y$是反应速度。

今测得一组数据入下表，试由此确定参数$\beta_1,\cdots,\beta_5$，并给出其置信区间

$\beta_1,\cdots,\beta_5$的参考值为$[0.1,0.05,0.02,1,2]$

此处可以看作给出了回归系数的初始值，若没有提供则可以随机取

Data:

序号 反应速度 氢 戊烷 异构戊烷
1 8.55 470 300 10
2 3.79 285 80 10
3 4.82 470 300 120
4 0.02 470 80 120
5 2.75 470 80 10
6 14.39 100 190 10
7 2.54 100 80 65
8 4.35 470 190 65
9 13.00 100 300 54
10 8.50 100 300 120
11 0.05 100 80 120
12 11.32 285 300 10
13 3.13 285 190 120

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
clear,clc
%% Data
yr = [8.55 3.79 4.82 0.02 2.75 14.39 ...
      2.54 4.35 13.00 8.50 0.05 11.32 3.13]';
x1 = [470 285 470 470 470 100 100 ...
      470 100 100 100 285 285]';
x2 = [300 80 300 80 80 190 80 190 ...
      300 300 80 300 190]';
x3 = [10 10 120 120 10 10 65 ...
      65 54 120 120 10 120]';
xr = [x1, x2, x3];
beta0 = [0.1,0.05,0.02,1,2]';
%% Regress
% 用匿名函数定义拟合函数
x = xr;
hs = @(beta,x) ((beta(4)*x(:,2)-(x(:,3))/(beta(5))) ...
               ./(1+beta(1)*x(:,1)+beta(2)*x(:,2)+beta(3)*x(:,3)));
% 引入初始值
beta0;
% [nlinfit]非线性回归计算回归系数beta，其中r和j是后续命令所需的中间量
[beta, r, j] = nlinfit(xr,yr,hs,beta0)
% [nlparci]计算回归系数的置信区间
betaci = nlparci(beta,r,"jacobian",j)
% [nlpredci]计算预测y值及其置信区间半径
[y_pred, delta] = nlpredci(hs,xr,beta,r,'jacobian',j)

%% Disp
% [nlintool]展示拟合图像结果的交互式画面
nlintool(xr,yr,hs,beta0)
% rmse

output

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
beta =
    0.0628
    0.0400
    0.1124
    1.2526
    1.1914
betaci =
   -0.0377    0.1632
   -0.0312    0.1113
   -0.0609    0.2857
   -0.7467    3.2519
   -0.7381    3.1208
r =
    0.1321
   -0.1642
   -0.0909
    0.0310
    0.1142
    0.0498
   -0.0262
    0.3115
   -0.0292
    0.1096
    0.0716
   -0.1501
   -0.3026
y_pred =
    8.4179
    3.9542
    4.9109
   -0.0110
    2.6358
   14.3402
    2.5662
    4.0385
   13.0292
    8.3904
   -0.0216
   11.4701
    3.4326
delta =
    0.2805
    0.2474
    0.1766
    0.1875
    0.1578
    0.4236
    0.2424
    0.1638
    0.3426
    0.3281
    0.3699
    0.3237
    0.1749

Pic