Matlab实现多元的回归实例

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）一般说来，R2值越大越好。

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

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

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

例1（Hamilton，1987）数据如下：

序号 Y X1 X2 精彩文档

实用标准文案

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 第一步 分析数据

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 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 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 在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 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）表示画出残差与残差区间的杠杆图。执行后得到： b = -4.5154

精彩文档

实用标准文案

3.0970 1.0319 bint =

-4.6486 -4.3822 3.0703 3.1238 1.0238 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

精彩文档

1.0399 实用标准文案

rint =

-0.0087 0.0314 -0.0303 0.0128 -0.0301 0.0098 -0.0299 0.0162 -0.0106 0.0308 -0.0313 0.0102 -0.0252 0.0178 -0.0299 0.0089 -0.0174 0.0272 -0.0331 0.0058 -0.0161 0.0275 -0.0027 0.0354 -0.0236 0.0190 -0.0079 0.0299 -0.0156 0.0298 stats = 1.0e+004 *

0.0001 3.9222 0 0.0000 即，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

在生产中测得13组数据：

序号 1 2 3 4 5 6 7 8 9 10 11 12 13 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 求出关系式YfX。

精彩文档

实用标准文案

解：（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];

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 -3.2331 -5.3126 -6.5603 精彩文档

4.0485 6.2555 1.9707 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 -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)

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

精彩文档

实用标准文案

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 HistoryRMSE10

0123

最后得到回归方程（蓝色行是被保留的有效行，红色行表示被剔除的变量）：

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

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

精彩文档