课件-5-Matlab软件包与多元回归(共17页)

精选优质文档-倾情为你奉上Matlab软件包与多元回归 回归分析的方法以及“回归（Regress）”这个名称的起源，统计学史上一般归功于英国生物学家、统计学家Francis Galton（著名的Darwin的表哥，cousin of ）。 Francis Galton（18221911） K. Pearson（18571936） 问题：考察人体的某项指标，比如说：身高。假设父辈的身高为，孩子的身高为，身高具有遗传性，但是，父辈的身高又不能完全决定孩子的身高。局部说来，父辈是高个子，子女也是高个子，父辈是矮个子，子女也是矮个子。总体观察，这个结论不对（否则，长期进化之后，人类的身高应该两极分化：很高的人与很矮的人）。F Galton的学生K Pearson观察了1078对夫妇，以每对夫妇的平均身高为，取他们的一个成年儿子的身高为，全体父母身高的平均值吋。K Pearson观察得到：父母的身高中等，即，（吋）时，（吋），看起来遗传起作用；父母的身高较高时，（吋）时，（吋），子女的身高变小了；父母的身高较矮时，（吋）时，（吋），子女的身高变大了。观察说明，父母的身高到了高和矮的边界时，总有某种力量将子女的身高拉向中心平均值，这就叫回归（Regress）。正是由于这种现象，Galton用“Regress”一词来描述与的关系，尽管，这种“向中心回归”的现象只在特殊的领域观察得到，不具有普遍性，但是，统计学界使用习惯了，所以，沿用至今。回归（Regress）的这种性质，可以得到理论上的说明。设，的相关系数，则可以证明，线性回归方程中的系数满足，当时，这一点解释了Galton指出的“Regress”现象：父辈身高的方差，子女身高的方差，在一代之间变化不大，故可以假设。于是，令，从线性回归方程得到：，以及，这样，就可以得到对于具体的，有，从可以看出：当时，尽管，也增加，但是，增加一个单位，并不相应增加一个单位，相当于打了折，这就解释了“向中心回归”的现象。注：在许多实际问题中有，此时，“向中心回归”的现象就不存在了。（一）一般多元回归一般在生产实践和科学研究中，人们得到了参数和因变量的数据，需要求出关系式，这时就可以用到回归分析的方法。如果只考虑是线性函数的情形，当自变量只有一个时，即，中时，称为一元线性回归，当自变量有多个时，即，中时，称为多元线性回归。进行线性回归时，有4个基本假定： 待定参数（系数）是线性关系； 残差是独立的； 残差满足正态分布。 残差满足方差奇性（所谓方差齐性指的就是我们要比较的几组数据是独立的、且服从同方差的正态分布）；在Matlab软件包中有一个做一般多元回归分析的命令regress，调用格式如下：b, bint, r, rint, stats = regress(y,X,alpha) 或者b, bint, r, rint, stats = regress(y,X) 此时，默认置信度alpha 0.05。这里，y是一个的列向量，X是一个的矩阵，其中第一列是全1向量（这一点对于回归来说很重要，这一个全1列向量对应回归方程的常数项），一般情况下，需要人工造一个全1列向量。回归方程具有如下形式（待定参数具有线性关系）：其中，是残差。在返回项b,bint,r,rint,stats中，是回归方程的系数；是一个矩阵，它的第行表示的(1-alpha)可信区间；是的残差列向量；是矩阵，它的第行表示第个残差的(1-alpha)可信区间；注释：残差与残差区间杠杆图，最好在0点线附近比较均匀的分布，而不呈现一定的规律性，如果是这样，就说明回归分析做得比较理想。 一般的，返回4个值：值、F_检验值、阈值，与显著性概率相关的值（如果这个值不存在，则，只输出前3项）。注释：（1）一般说来，值越大越好。（2）人们一般用以下统计量对回归方程做显著性检验：F_检验、t_检验、以及相关系数检验法。Matlab软件包输出F_检验值和阈值。一般说来，F_检验值越大越好，特别的，应该有F_检验值。 我国著名统计学家许宝禄（19101970）教授证明：F检验有多方面的优良性。（3）与显著性概率相关的值应该满足。如果，则说明回归方程中有多余的自变量，可以将这些多余的自变量从回归方程中剔除（见下面逐步回归的内容）。这几个技术指标说明拟合程度的好坏。这几个指标都好，就说明回归方程是有意义的。例1（Hamilton，1987）数据如下：序号YX1X2112.372.239.66212.662.578.94312.003.874.40411.933.106.64511.063.394.91613.032.838.52713.133.028.04811.442.149.05912.863.047.711010.843.265.111111.203.395.051211.562.358.511310.832.766.591412.633.904.901512.463.166.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.0000ans = 1.0000 0.4341 0.4341 1.0000即，corrcoef(x1,y)0.0025，corrcoef(x2,y)0.4341，说明没有非常明显的单变量线性关系。图形如下：也看不出有线性关系，但是，旋转图形，可以看出所有点几乎在一个平面上。这说明，在一个平面上，满足线性关系：或者，换成一个常见的形式 其中，是残差。于是，在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.0319bint = -4.6486 -4.3822 3.0703 3.1238 1.0238 1.0399r = 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.0071rint = -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.0298stats = 1.0e+004 * 0.0001 3.9222 0 0.0000即，。可信度95，且，与显著性概率相关的，这说明，回归方程中的每个自变量的选取，都是有意义的。残差杠杆图：从杠杆图看出，所有的残差都在0点附近均匀分布，区间几乎都位于之间，即，没有发现高杠杆点，也就是说，数据中没有强影响点、异常观测点。 综合起来看，以上回归结果（回归函数、拟合曲线或曲面）近乎完美。（二）回归诊断的必要性例2（Anscombe，1973）数据如下：序号108.049.147.468.06.5828.06.958.146.7785.76313.07.588.7412.7487.71498.818.777.1188.845118.339.267.8188.476149.968.18.8487.04767.246.136.0885.25844.263.105.391912.591210.849.138.1585.561074.827.266.4487.911155.684.745.7386.89解：A 在Matlab7.1中，Import Data：A=上表中的数据。（1）画图观察，并做与的线性回归画图观察x1=A(:,2);y1=A(:,3);plot(x1,y1,-*)旋转图形，看出数据点都在一个平面上，但不一定在同一条直线上。做线性回归：x1=A(:,2);y1=A(:,3);plot(x1,y1,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y1,x)rcoplot(r,rint)var(r)程序执行后，得到结果：b = 5.2818 0.2743即，系数的95%可信区间为：bint = 2.6131 7.9504 -0.0206 0.5691回归误差的方差为：Var= 2.7663拟合指标为：stats = 0.3297 4.4277 0.0647 3.0737即，与显著性概率相关的，这说明，拟合效果不好。残差杠杆图为：残差杠杆图显示，第一个数据有异常。（2）画图观察，并做与的线性回归画图观察x1=A(:,2);y2=A(:,4);plot(x1,y1,-*) 有一个数据点不在直线上。做线性回归x1=A(:,2);y2=A(:,4);plot(x1,y2,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y2,x)rcoplot(r,rint)程序执行后得到：b = 5.7742 0.2134即，回归误差的方差为：var = 3.3036拟合指标不太好：stats = 0.1996 2.2447 0.1683 3.6707残差杠杆图为：残差杠杆图显示，有2个数据异常。（3）画图观察，并做与的线性回归画图观察 做线性回归x1=A(:,2);y3=A(:,5);plot(x1,y1,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y3,x)rcoplot(r,rint)var(r)结果如下：b = 5.0264 0.3059即，回归误差的方差为：Var=2.4250拟合指标不太好：stats = 0.4112 6.2846 0.0335 2.6944（4）画图观察，并做与的线性回归画图观察 做线性回归b = 3.0017 0.4999即，回归误差的方差为：Var=1.3742拟合指标不太好：stats = 0.6667 18.0033 0.0022 1.5269残差杠杆图为：B 在Mathematica5.0中，做最小二乘拟合：A=0.0000,8.0400,9.1400,7.4600,8.0000,6.5800,8.0000,6.9500,8.1400,6.7700,8.0000,5.7600,13.0000,7.5800,8.7400,12.7400,8.0000,7.7100,9.0000,8.8100,8.7700,7.1100,8.0000,8.8400,11.0000,8.3300,9.2600,7.8100,8.0000,8.4700,14.0000,9.9600,8.1000,8.8400,8.0000,7.0400,6.0000,7.2400,6.1300,6.0800,8.0000,5.2500, 4.0000,4.2600,3.1000,5.3900,19.0000,12.5000,12.0000,10.8400,9.1300,8.1500,8.0000,5.5600,7.0000,4.8200,7.2600,6.4400,8.0000,7.9100,5.0000,5.6800,4.7400,5.7300,8.0000,6.8900;n=LengthA;x1=TableAj,1,j,n;y1=TableAj,2,j,n;d1=Threadx1,y1;f1=Fitd1,1,x,x;Printf1 = ,f1y2=TableAj,3,j,n;d2=Threadx1,y2;f2=Fitd2,1,x,x;Printf2 = ,f2y3=TableAj,4,j,n;d3=Threadx1,y3;f3=Fitd3,1,x,x;Printf3 = ,f3x2=TableAj,5,j,n;y4=TableAj,6,j,n;d4=Threadx2,y4;f4=Fitd4,1,x,x;Printf4 = ,f4程序执行后，得到结果：f1 = 5.28176 +0. xf2 = 5.77417 +0. xf3 = 5.02645 +0. xf4 = 3.00173 +0. x注：这个例子说明，应该对数据拟合、回归结果作出评价，即，要做回归诊断。（三）逐步回归假设已有数据X 和Y，在Matlab软件包中，使用stepwise命令进行逐步回归，得到回归方程，其中是随机误差。stepwise命令的使用格式如下：stepwise(X,Y)注意：应用stepwise命令做逐步回归，数据矩阵X的第一列不需要人工加一个全1向量，程序会自动求出回归方程的常数项（intercept）。在应用stepwise命令进行运算时，程序不断提醒将某个变量加入（Move in）回归方程，或者提醒将某个变量从回归方程中剔除（Move out）。注释：使用stepwise命令进行逐步回归，既有剔除变量的运算，也有引入变量的运算，它是目前应用较为广泛的一种多元回归方法。在运行stepwise(X,Y)命令时，默认显著性概率。 例2（Hald,1960）Hald数据是关于水泥生产的数据。某种水泥在凝固时放出的热量（单位：卡克）与水泥中4种化学成分所占的百分比有关：在生产中测得13组数据：序号X1X2X3X4Y172666078.52129155274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.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;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.1441bint = -99.1786 223.9893 -0.1663 3.2685 -1.1589 2.1792 -1.6385 1.8423 -1.7791 1.4910r = 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.2943rint = -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 -2.2372 6.2191 -4.1338 6.0797 -6.9115 2.3228stat = 0.9824 111.4792 0.0000 5.9830以及残差杠杆图：于是，我们得到：并且，残差杠杆图显示，残差均匀分布在0点线附近，在stat返回的4个值中，0.9824，说明模型拟合的很好。F_检验值111.47920.000，符合要求。但是，与显著性概率相关的值5.98300.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)程序执行后得到下列逐步回归的画面： Move x4 in程序提示：将变量x4加进回归方程（Move x4 in），点击Next Step按钮，即，进行下一步运算，将第4列数据对应的变量加入回归方程。点击Next Step按健后，又得到提示：将变量x1加进回归方程（Move x1 in），点击 Next Step按钮，即，进行下一步运算，将第1列数据对应的变量加入回归方程。点击Next Step按健后，又得到提示：Move No terms，即，没有需要加入（也没有需要剔除）的变量了。 注意：在Matlab7.0软件包中，可以直接点击“All Steps”按钮，直接求出结果（省略中间过程）。intercept 最后得到回归方程（蓝色行是被保留的有效行，红色行表示被剔除的变量）：回归方程中录用了原始变量和。模型评估参数分别为：，修正的值，F_检验值176.627，与显著性概率相关的值，残差均方RMSE2.73427（这个值越小越好）。以上指标值都很好，说明回归效果比较理想。 另外，截距intercept=103.097（Intercept = the estimated value of the constant term），这就是回归方程的常数项。我们将的数据放在一起画图观察：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)程序执行后得到图形：不断旋转画面观察，图形大约是一个平面。专心-专注-专业