实验五 基于matlab 的数据拟合与回归分析一、 实验类型验证性实验二、 实验学时2学时三、 实验目的1、掌握利用MATLAB 中实现单因变量的多元线性回归分析(经典多元线性回归分析)的方法；2、掌握利用MATLAB 中实现多因变量的多元线性回归分析(多对多线性回归分析)的方法。
3、掌握MATLAB 有关逐步回归的命令。
四、 实验内容及要求 实验内容：1 学习MATLAB 中有关经典多元线性回归分析的命令；(1) [b,bint,r,rint,stats] = regress(y,X,alpha)建立回归分析模型01122,i i i ip ip i y b b x b x b x ε=+++++其中()()2E 0,D i i εεσ==,1,2,,i n =.注：在该命令中，设计矩阵()(1)ij n p X x ⨯+=(X 的第1列全为1 ) 或响应值向量()1i n y y ⨯=中的无效值Nan 将被免忽略。
输入参数至少有2个，alpha 是检验的显著性水平，默认值为0.05。
输出参数至少有1个，b 是回归系数的估计值；bint 是各回归系数的置信度为1-alpha 的置信区间(第1列是它们的下界，第2列是它们的上界)；r 是残差，rint 是残差的置信区间；stats 给出一个1*4的向量，依次是判定系数2R ，F 统计量的观测值及检验的p 值，以及误差的方差的估计值2ˆσ. 这里//(1)U p MMSF Q n p MSE==-- 是检验回归模型是否显著的检验统计量，当012:0p H b b b ====成立时， ~(,1)F F p n p --.可以通过分析残差的置信区间来观察原始数据是否存在异常点，若残差的置信区间不包括0，则可认为对应的观测值是异常值，将异常值全部去掉，重新建立的回归模型将更加精确。
可以调用命令rcoplot(r,rint)按观测顺序逐个画出各组观测值对应的残差和残差置信区间，从而观测异常值的位置。
(2) stats = regstats(y,X,model ,whichstats ) % 根据指定模型model 建立回归方程输入参数：响应变量()1i n y y ⨯=，设计变量()ijn pX x ⨯=，model 参数表5.1 model 参数及说明包括系数的估计值、因变量的预测值、残差、判定系数等23个选项供用户选择，通过这一界面，matlab将用户所选择的选项对应的回归分析结果显示在matlab工作空间中。
图5.1 regstats弹出的图形用户界面当带有输出选项stats时，直接返回一个结构体stats,它有24字段，包含了图形界面中的所有选项，此时，不弹出用户界面。
若输入中带有第4个参数whichstats,则stats仅包含whichstats 指定的统计量。
表5.2 stats字段名称及说明(3) 乘方法。
输入参数至少两个，X(n p ⨯)是设计矩阵，y(1n ⨯)响应变量的观测值。
wfun 给出权函数的名称，tune 给出权函数相应的调节系数，const 控制模型中是否包含常数项。
若const=’on’或1,则模型中包含常数项，此时自动在X 的第1列加入1列1，若const=’off’或0，则不包含常数项。
表5.3 权函数与调节系数其中r =，resid 是上一步迭代的残差，/0.6745s MAD =，MAD 是残差绝对值的中位数，若X 有p 列，计算MAD 时应残差绝对值的前p 个值最小值去掉。
h 是最小二乘拟合得到的中心化杠杆值。
输出：至少1个参数，b 是回归模型中的系数向量；stats 返回用于进行模型诊断的结构体变量。
其字段名称和意义如下：表5.4 stats 的字段名称及说明2 [beta,Sigma,E,CovB,logL]= mvregress(X,Y ,Name,Value) 建立多对多的线性回归模型 输入：至少两个参数。
d 维响应变量的观测值矩阵()ij n d Y y ⨯=，设计矩阵X 可以是矩阵或元胞数组，当1d =时，()ijn kX x ⨯=，当模型包含常数项时，X 的第1列应全为1；当d>1，并且Y 的每一个维度有相同的设计矩阵时，()ijn kX x ⨯=或是一个只包含单个n k ⨯的设计矩阵的元胞数组；当 d>1，并且Y 的每一个维度有不同的设计矩阵时，X 将是一个长度为d 的元胞数组，每个元胞都是一个n k ⨯的设计矩阵。
Name 和Value 参数可指定算法与相应的一些选项。
其值如表5.5所示。
表5.5 mvregress 函数的name 选项及其值(部分)的协方差矩阵的估计；E 是残差矩阵；CovB 是beta 的协方差矩阵；,logL 是最后一次迭代对应的对数似然函数的值。
3 学习MATLAB 有关逐步回归的命令。
表5.6 逐步回归命令函数4 (1) 设某种水泥在凝固时所释放的热量Y(卡/克)与水泥中下列4种化学成分有关1x 为233CaO Al O ⋅的成分(%)；2x 为23CaO SiO ⋅的成分(%)； 3x 为23234CaO Al O Fe O ⋅⋅的成分(%)；4x 为22CaO SiO ⋅的成分(%)。
共观测了13组数据(见表5.7)。
试求出Y 与1234,,,x x x x 的回归方程，并对该回归方程和各个回归系数进行检验。
表5.7 水泥数据(2) 设发电量1，工业总产值2与钢材产量1，水泥产量2，机械工业总产值3，棉纱产量4x ，机制纸产量5x 之间有线性相关关系。
现收集了1949到1978年共30年的数据（见表4.2）。
试建立12,Y Y 与12345,,,,x x x x x 的线性模型。
表5.8 发电量与经济发展数据（1）实验过程要求用Matlab 软件完成；（2）实验结果逐个导入Word 文档，并按问题作出解释； （3）实验报告按照既定格式书写。
五、 实验的重点和难点实验重点： MA TLAB 中关于多元线性回归模型的实现方法。
实验难点： 无六、 实验准备1、实验室电脑要求安装有Matlab软件；2、上实验课程的学生要对涉及到的统计概念有所了解；3、要求学生事先查阅并熟悉Matlab的相关命令。
七、实验程序及结果分析1 学习MATLAB中有关经典多元线性回归分析的命令（1）程序如下：y=[7613.51 7850.91 8381.86 9142.81 10813.6 8631.43 8124.94 9429.79 10230.81 10163.61 9737.56 8561.06 7781.82 7110.97]';x1=[7666 7704 8148 8571 8679 7704 6471 5870 5289 3815 3335 2927 2758 2591]'; x2=[16.22 16.85 17.93 17.28 17.23 17 19 18.22 16.3 13.37 11.62 10.36 9.83 9.25]'; X=[ones(size(y)) x1.^2 x2.^2 x1 x2 x1.*x2];[b,bint,r,rint,stats] = regress(y,X)结果：结果说明：b是回归系数的估计值；bint是各回归系数的置信度为1-alpha的置信区间(第1列是它们的下界，第2列是它们的上界)；r是残差，rint是残差的置信区间；stats给出一个1*4的向量，依R，F统计量的观测值及检验的p值，以及误差的方差的估计值2ˆ .次是判定系数2(2)程序如下：x=[ones(10,1) (1:10)'];y=x*[10;1]+normrnd(0,0.1,10,1);[b,bint]=regress(y,x,0.05)结果：（3）程序如下：y=[3.73 6.27 5.93 5.77 5.72 5.80 5.87 5.78 5.96 ...6.03 6.507.00 6.80 6.68 7.03 7.67 7.59 6.96 ...7.17 6.99 7.17 6.99 6.64 6.71 7.01 7.40 7.49 ...7.75 8.17 8.09 8.11 8.48 8.99];x=1:33;scatter(x,y);%画图，离散点hold on;p=regress(y',[ones(33,1) x']);%线性拟合r=robustfit(x,y);%稳健拟合plot(x,p(1)+p(2)*x,':',x,r(1)+r(2)*x,'r');%画图legend('原始数据','线性拟合','稳健拟合')%设置图例结果：2. 学习MATLAB有关多对多的多元线性回归分析的命令程序如下：x=[143 144 145 147 148 150 153 154 155 156 157 158 159 160 161 162]'; X=[ones(16,1),x];Y=[87 85 88 91 92 90 93 95 98 98 97 95 97 99 100 102]';[beta,Sigma,E,CovB,logL]= mvregress(X,Y,'covtype','full')%建立多对多的线性回归模型结果：结果说明：beta是回归方程的系数估计矩阵；Sigma是对响应向量的协方差矩阵的估计；E是残差矩阵；CovB是beta的协方差矩阵；,logL是最后一次迭代对应的对数似然函数的值。
3.学习MATLAB有关逐步回归的命令。
程序如下：y=[11.2 13.4 40.7 5.3 24.8 12.7 20.9 35.7 8.7 ...9.6 14.5 26.9 15.7 36.2 18.1 28.9 14.9 25.8 21.7 25.7];%犯罪率y，低收入家庭百分比x1，失业率x2，人口总数x3x1=[16.5 20.5 26.3 16.5 19.2 16.5 20.2 21.3 17.2 14.3 ...18.1 23.1 19.1 24.7 18.6 24.9 17.9 22.4 20.2 16.9];x2=[6.2 6.4 9.3 5.3 7.3 5.9 6.4 7.6 4.9 6.4 6.0 7.4 ...5.8 8.66.5 8.3 6.7 8.6 8.4 6.7];x3=[587 643 635 692 1248 643 1964 1531 713 749 7895 ...762 2793 741 625 854 716 921 595 3353];figure(1),plot(x1,y,'*');figure(2),plot(x2,y,'*');figure(3),plot(x3,y,'*');X1=[x1',x2',x3'];stepwise(X1,y)结果：4.(1)程序如下：a=[7 26 6 60;1 29 15 52;11 56 8 20;11 31 8 47;7 52 6 33;11 55 9 22;1 31 22 44;1 31 22 44;2 54 18 22;21 47 4 26;1 40 23 34;11 66 9 12;10 68 8 12];x1=a(:,1);x2=a(:,2);x3=a(:,3);x4=a(:,4);y=[78.5;74.3;104.3;87.6;95.9;109.2;72.5;72.5;93.1;115.9;83.8;113.3;109.4]; X=[ones(size(x1)),x1,x2,x3,x4]；[b,bint,r,rint,stats] = regress(y,X)结果如下：b =69.76611.41180.4621-0.0399-0.2226bint =-97.7881 237.3202-0.3563 3.1798-1.2765 2.2008-1.8203 1.7406-1.9175 1.4722r =0.43281.8950-2.1031-1.2393-0.19323.7441-2.3311-2.33111.17110.71442.62350.5346-2.9178rint =-3.7577 4.6233-3.1648 6.9548-5.8987 1.6924-6.4122 3.9337-5.3272 4.9408-1.2481 8.7364-7.3498 2.6876-7.3498 2.6876-3.9131 6.2554-2.6165 4.0453-2.0795 7.3264-4.7525 5.8216-7.2782 1.4427stats =0.9833 117.4117 0.0000 6.5281说明：b 是线性方程的系数估计值，并且第一值表示常数系数。