Matlab实现多元回归实例

（一）一般多元回归

一般在生产实践和科学研究中，人们得到了参数xx1,,xn和因变量y的数据，需要求出关系式yfx，这时就可以用到回归分析的方法。如果只考虑

f是线性函数的情形，当自变量只有一个时，即，xx1,,xn中n1时，称

为一元线性回归，当自变量有多个时，即，xx1,,xn中n2时，称为多元线性回归。

进行线性回归时，有4个基本假定： ① 因变量与自变量之间存在线性关系； ② 残差是独立的； ③ 残差满足方差奇性； ④ 残差满足正态分布。

在Matlab软件包中有一个做一般多元回归分析的命令regeress，调用格式如下：

[b, bint, r, rint, stats] = regress(y,X,alpha) 或者

[b, bint, r, rint, stats] = regress(y,X) 此时，默认alpha ＝ 0.05. 这里，y是一个n1的列向量，X是一个nm1的矩阵，其中第一列是全1向量（这一点对于回归来说很重要，这一个全1列向量对应回归方程的常数项），一般情况下，需要人工造一个全1列向量。回归方程具有如下形式：

y01x1mxm

其中，是残差。

在返回项[b,bint,r,rint,stats]中， ①b01m是回归方程的系数；

②bint是一个m2矩阵，它的第i行表示i的(1-alpha)置信区间； ③r是n1的残差列向量；

④rint是n2矩阵，它的第i行表示第i个残差ri的(1-alpha)置信区间； 注释：残差与残差区间杠杆图，最好在0点线附近比较均匀的分布，而不呈现一定的规律性，如果是这样，就说明回归分析做得比较理想。

⑤ 一般的，stast返回4个值：R2值、F_检验值、阈值f，与显著性概率相关的p值（如果这个p值不存在，则，只输出前3项）。注释：

1

（1）一般说来，R2值越大越好。

（2）人们一般用以下统计量对回归方程做显著性检验：F_检验、t_检验、以及相关系数检验法。Matlab软件包输出F_检验值和阈值f。一般说来，F_检验值越大越好，特别的，应该有F_检验值f。

（3）与显著性概率相关的p值应该满足palpha。如果palpha，则说明回归方程中有多余的自变量，可以将这些多余的自变量从回归方程中剔除（见下面逐步回归的内容）。

这几个技术指标说明拟合程度的好坏。这几个指标都好，就说明回归方程是有意义的。

例1（Hamilton，1987）数据如下： 序号 Y X1 X2 1 12.37 2.23 9.66 2 12.66 2.57 8.94 3 12.00 3.87 4.40 4 11.93 3.10 6.64 5 11.06 3.39 4.91 6 13.03 2.83 8.52 7 13.13 3.02 8.04 8 11.44 2.14 9.05 9 12.86 3.04 7.71 10 10.84 3.26 5.11 11 11.20 3.39 5.05 12 11.56 2.35 8.51 13 10.83 2.76 6.59 14 12.63 3.90 4.90 15 12.46 3.16 6.96 第一步 分析数据 在Matlab软件包中分析是否具有线性关系，并作图观察，M—文件opt_hanmilton_1987：

x1=[2.23,2.57,3.87,3.10,3.39,2.83,3.02,2.14,3.04,3.26,3.39,2.35,2.76,3.90,3.16]; x2=[9.66,8.94,4.40,6.64,4.91,8.52,8.04,9.05,7.71,5.11,5.05,8.51,6.59,4.90,6.96]; y=[12.37,12.66,12.00,11.93,11.06,13.03,13.13,11.44,12.86,10.84,11.20,11.56,10.83,12.63,12.46]; corrcoef(x1,y); corrcoef(x2,y); plot3(x1,x2,y,'*');

得到结果： ans =

1.0000 0.0025

2

0.0025 1.0000 ans =

1.0000 0.4341 0.4341 1.0000

即，corrcoef(x1,y)＝0.0025，corrcoef(x2,y)＝0.4341，说明没有非常明显的单变量线性关系。图形如下：

也看不出有线性关系，但是，旋转图形，可以看出所有点几乎在一个平面上。

这说明，y,x1,x2在一个平面上，满足线性关系：

a1x1a2x2bya

或者，换成一个常见的形式 ya0a1x1a2x2

其中，是残差。于是，在Matlab软件包中做线性多元回归，写一个M—文件opt_regress_hamilton：

x1=[2.23,2.57,3.87,3.10,3.39,2.83,3.02,2.14,3.04,3.26,3.39,2.35,2.76,3.90,3.16]'; x2=[9.66,8.94,4.40,6.64,4.91,8.52,8.04,9.05,7.71,5.11,5.05,8.51,6.59,4.90,6.96]'; y=[12.37,12.66,12.00,11.93,11.06,13.03,13.13,11.44,12.86,10.84,11.20,11.56,10.83,12.63,12.46]'; e=ones(15,1); x=[e,x1,x2];

[b,bint,r,rint,stats]=regress(y,x,0.05) rcoplot(r,rint)

其中，rcoplot（Residual case order plot）表示画出残差与残差区间的杠杆图。执

3

行后得到： b =

-4.5154 3.0970 1.0319 bint =

-4.6486 -4.3822 3.0703 3.1238 1.0238 1.0399 r =

0.0113 -0.0087 -0.0102 -0.0069 0.0101 -0.0106 -0.0037 -0.0105 0.0049 -0.0136 0.0057 0.0163 -0.0023 0.0110 0.0071 rint =

-0.0087 -0.0303 -0.0301 -0.0299 -0.0106 -0.0313 -0.0252 -0.0299 -0.0174 -0.0331 -0.0161 -0.0027 -0.0236 -0.0079 -0.0156 stats =

1.0e+004 *

0.0001

0.0314 0.0128 0.0098 0.0162 0.0308 0.0102 0.0178 0.0089 0.0272 0.0058 0.0275 0.0354 0.0190 0.0299 0.0298 3.9222 0 0.0000

4

即，y4.5153.097x11.0319x2。

置信度95％，且R21.0,F_检验值392220，与显著性概率0.05相关的

p0.00000.05，这说明，回归方程中的每个自变量的选取，都是有意义的。

残差杠杆图：

从杠杆图看出，所有的残差都在0点附近均匀分布，区间几乎都位于0.03,0.03之间，即，没有发现高杠杆点，也就是说，数据中没有强影响点、异常观测点。 综合起来看，以上回归结果（回归函数、拟合曲线或曲面）近乎完美。

（二）逐步回归

假设已有数据X 和Y，在Matlab软件包中，使用stepwise命令进行逐步回归，得到回归方程Ya1X1a2X2anXn，其中是随机误差。stepwise命令的使用格式如下：stepwise(X,Y)

注意：应用stepwise命令做逐步回归，数据矩阵X的第一列不需要人工加一个全1向量，程序会自动求出回归方程的常数项（intercept）。

在应用stepwise命令进行运算时，程序不断提醒将某个变量加入（Move in）回归方程，或者提醒将某个变量从回归方程中剔除（Move out）。

注释：①使用stepwise命令进行逐步回归，既有剔除变量的运算，也有引入变量的运算，它是目前应用较为广泛的一种多元回归方法。②在运行stepwise(X,Y)命令时，默认显著性概率0.05。

例2（Hald,1960）Hald数据是关于水泥生产的数据。某种水泥在凝固时放出的热量Y（单位：卡／克）与水泥中4种化学成分所占的百分比有关：

x1:3CaoAl2o3x2:3CaoSio2x3:4CaoAl2o3Fe2o3x4:2CaoSio2

5

在生产中测得13组数据： 序号 X1 1 7 2 1 3 11 4 11 5 7 6 11 7 3 8 1 9 2 10 21 11 1 12 11 13 10 求出关系式YfX。

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 解：（1）本问题涉及的数据是5维的，不能画图观察。先做异常值分析。 X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12;10,68,8,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]'; A=[X,Y]; mahal(A,A)

程序执行后得到结果： ans =

5.6803 3.6484 6.7002 3.3676 3.3839 4.4300 4.0080 6.5067 3.0849 7.5016 5.1768 2.4701

可以认为数据都是正常的。 （2）一般多元回归。

在Matlab软件包中写一个M—文件opt_cement_1:

X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47; 7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44; 2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12; 10,68,8,12];

6

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]'; a1=ones(13,1); A=[a1,X];

[b,bint,r,rint,stat]=regress(Y,A) rcoplot(r,rint)

程序执行后得到： b =

62.4054 1.5511 0.5102 0.1019 -0.1441 bint =

-99.1786 223.9893 -0.1663 3.2685 -1.1589 2.1792 -1.6385 1.8423 -1.7791 1.4910 r =

0.0048 1.5112 -1.6709 -1.7271 0.2508 3.9254 -1.4487 -3.1750 1.3783 0.2815 1.9910 0.9730 -2.2943 rint =

-4.0390 4.0485 -3.2331 6.2555 -5.3126 1.9707 -6.5603 3.1061 -4.5773 5.0788 -0.5623 8.4132 -6.0767 3.1794 -6.8963 0.5463 -3.5426 6.2993 -3.0098 3.5729

7

-2.2372 6.2191 -4.1338 6.0797 -6.9115 2.3228 stat =

0.9824 111.4792 0.0000 5.9830 以及残差杠杆图：

于是，我们得到：

Y62.40541.5511x10.5102x20.1019x30.1441x4

并且，残差杠杆图显示，残差均匀分布在0点线附近，在stat返回的4个值中，

R2＝0.9824，说明模型拟合的很好。F_检验值＝111.4792>0.000，符合要求。但是，与显著性概率相关的p值＝5.9830>0.05，这说明，回归方程中有些变量可以剔除。

（3）逐步回归

在Matlab软件包中写一个M—文件opt_cement_2: X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47; 7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44; 2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12; 10,68,8,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]'; stepwise(X,Y)

程序执行后得到下列逐步回归的画面：

8

Coefficients with Error BarsX1 Coeff. t-stat p-val 1.86875 3.5500 0.0046X2 0.789125 4.6862 0.0007X3 -1.25578 -2.0984 0.0598Move x4 in X4-3-2-10123 -0.738162 -4.7748 0.000617Model HistoryRMSE1615141

程序提示：将变量x4加进回归方程（Move x4 in），点击Next Step按钮，即，进行下一步运算，将第4列数据对应的变量x4加入回归方程。点击Next Step按健后，又得到提示：将变量x1加进回归方程（Move x1 in），点击 Next Step按钮，即，进行下一步运算，将第1列数据对应的变量x1加入回归方程。点击Next Step按健后，又得到提示：Move No terms，即，没有需要加入（也没有需要剔除）的变量了。

注意：在Matlab7.0软件包中，可以直接点击“All Steps”按钮，直接求出结果（省略中间过程）。

Coefficients with Error BarsX1 Coeff. t-stat p-val 1.43996 10.4031 0.0000X2 0.41611 2.2418 0.0517X3intercept -1-0.500.511.5 -0.410043 -2.0581 0.0697X4 -0.613954 -12.6212 0.0000220Model HistoryRMSE10123 最后得到回归方程（蓝色行是被保留的有效行，红色行表示被剔除的变量）：

0 9

Y103.0971.43996x10.613954x4

回归方程中录用了原始变量x1和x4。模型评估参数分别为：

20.964212，F_检验值＝176.627，与显著性R20.972471，修正的R2值R概率相关的p值＝1.581061080.05，残差均方RMSE＝2.73427（这个值越小越好）。以上指标值都很好，说明回归效果比较理想。

另外，截距intercept=103.097（Intercept = the estimated value of the constant term），这就是回归方程的常数项。

我们将x1,x4,y的数据放在一起画图观察：

X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12;10,68,8,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]'-103.097;

plot3(X(:,1),X(:,4),Y)

程序执行后得到图形：

100-10-20-3060504030201020151050

不断旋转画面观察，图形大约是一个平面。

10