第八章 统计回归模型

第八章 统计回归模型时间：创作：欧阳法回归分析是研究一个变量 Y 与其它若干变量 X 之 间相关关系的一种数学工具 .它是在一组试验或观测 数据的基础上，寻找被随机性掩盖了的变量之间的依 存关系.粗略的讲，可以理解为用一种确定的函数关 系去近似代替比较复杂的相关关系 .这个函数称为回 归函数.回归分析所研究的主要问题是如何利用变量 X 、 Y 的观察值 （样本），对回归函数进行统计推断，包括 对它进行估计及检验与它有关的假设等.回归分析包含的内容广泛 . 此处将讨论多项式回 归、多元线性回归、非线性回归以及逐步回归.一、多项式回归（1） 一元多项式回归 一元多项式回归模型的一般形式为y P + P x +. + P xm + 01m如果从数据的散点图上发现y与x呈现较明显的 二次（或高次）函数关系，则可以选用一元多项式回归.1. 用函数 polyfit 估计模型参数，其具体调用格式 如下：p二polyfit(x,y,m) p返回多项式系数的估计值；m 设定多项式的最高次数；x, y为对应数据点值.p,S=polyfit(x,y,m) S 是一个矩阵，用来估计预测 误差.2. 输出预估值与残差的计算用函数 polyval 实 现，其具体调用格式如下：Y=polyval(p,X) 求 polyfit 所得的回归多项式在 X 处的预测值 Y.Y,DELTA=polyval(p,X,S) p，S 为 polyfit 的输 出， DELTA 为误差估计 .在线性回归模型中， Y DELTA以50%的概率包含函数在X处的真值.3模型预测的置信区间用polyconf实现，其具体 调用格式如下：Y,DELTA=polyconf(p,X,S,alpha) 求 polyfit 所得 的回归多项式在 X 处的预测值 Y 及预测值的显著性 为 1-alpha 的置信区间 YDELTA，alpha 缺省时为 0.05.4.交互式画图工具polytool，其具体调用格式如下：polytool(x,y,m)； polytool(x,y,m,alpha)； 用 m 次多项式拟合 x，y 的值，默认值为 1， alpha 为显著性水平，默认值为 0.05.例 1 观测物体降落的距离 s 与时间 t 的关系，得 到数据如下表，求 s.t (s)1/302/303/304/305/306/307/30s (cm)11.8615.6720.6026.6933.7141.9351.13t (s)8/309/3010/3011/3012/3013/3014/30s (cm)61.4972.9085.4499.08113.77129.54146.48解 根据数据的散点图，应拟合为一条二次曲线 .选用二次模型，具体代码如下：%输入数据t=1/30:1/30:14/30;s=11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.4972.90 85.44 99.08 113.77 129.54 146.48;%多项式系数拟合p,S=polyfit(t,s,2);则得回归模型为：s = 489.2946t 2 + 65.8896t + 9.1329 -%y 的拟合值及预测值 y 的置信半径 deltay,dalta=polyconf(p,t,S);得结果如下：y=Columns 1 through 1111.8729 15.7002 20.6148 26.6168 33.7060 41.882651.1465 61.4978 72.9363 85.4622 99.0754Columns 12 through 14113.7759 129.5637 146.4389 dalta=Columns 1 through 110.0937 0.0865 0.0829 0.0816 0.0817 0.0823 0.08270.0827 0.0823 0.0817 0.0816Columns 12 through 140.0829 0.0865 0.0937 %交互式画图 polytool（t,s,2）;polytool 所得的交互式图形如图 8-1 所示.图 8-1（2） 多元二项式回归 多元二项式回归模型的一般形式为 y P +B x +. + B x + Y B x x + 0 1 1 m m jk j k1 j ,k F (m,n-m-1)时拒绝H , F越大，1-a0说明回归方程越显著；与F对应的概率P a时拒绝H，回归模型成立；alpha表示显著性水平(缺省时为00.05).残差及其置信区间可以用命令rcoplot(r,rint)画出.例3已知某湖泊八年来湖水中COD浓度实测值(y)与影响因素，如湖区工业产值(x)、总人口数(x)、捕12鱼量牡)、降水量(x?的资料，建立y的水质分析模型.湖水浓度与影响因素数据表x11.3761.3751.3871.4011.4121.4281.4451.477X20.4500.4750.4850.5000.5350.5450.5500.575X32.1702.5542.6762.7132.8233.0883.1223.262X40.89221.16100.53460.95891.02391.04991.10651.1387y5.195.305.605.826.006.066.456.95解作出因变量y与各自变量的样本散点图作散点图的目的主要是观察因变量 y 与各自变量 间是否有比较好的线性关系，以便选择恰当的数学模 型形式.图 8-3、图 8-4、图 8-5、图 8-6 分别为 y 与 x 、x 、x 、x 的散点图.从图中可以看出这些点大致 1234 分布在一条直线旁边，因此有较好的线性关系，可以 采用线性回归.图8-3 y 与 x 的散点图图18-4 y 与 x 的散点图2图8-5 y 与 x 的散点图图38-6 y 与 x 的散点图4在Matlab中实现回归的具体代码如下：%输入数据x1=1.376 1.375 1.387 1.401 1.412 1.428 1.445 1.477;x2=0.450 0.475 0.485 0.500 0.535 0.545 0.550 0.575;x3=2.170 2.554 2.676 2.713 2.823 3.088 3.122 3.262; x4=0.8922 1.1610 0.5346 0.9589 1.0239 1.0499 1.1065 1.1387;x=ones(8,1) x1 x2 x3 x4;y=5.19 5.30 5.60 5.82 6.00 6.06 6.45 6.95;%多元线性回归b,bint,r,rint,stats=regress(y,x);得如下结果：b =-13.984913.19202.42280.0754-0.1897bint =-26.0019 -1.96791.4130 24.9711-14.2808 19.1264-1.4859 1.6366-0.9638 0.5844r =-0.06180.02280.01230.08900.0431-0.14730.01450.0274rint =-0.1130 -0.0107-0.1641 0.2098-0.1051 0.1297-0.2542 0.4321-0.0292 0.1153-0.2860 -0.0085-0.3478 0.3769-0.1938 0.2486stats =0.9846 47.9654 0.0047 0.0123故回归模型为y = 13.9849 +13.1920x + 2.4228x + 0.0754x - 0.1897x ,1234此外，由 stats 的值可知 r2 = 0.9846 , F 二 47.9654,P = 0.0047 。%残差分析，作残差图rcoplot(r,rint)图 8-7从残差图可以看出，除第一和第六个数据外，其 余数据的残差离零点均较近，且残差的置信区间均包 含零点.第一和第六个数据可视为异常点，将其去掉 后重新拟合可得新的回归模型.三、非线性回归非线性回归模型的一般形式为y = f (x, p ) +8 , x二(x , x ,x ) , p= (p , p，,p )，其中f对回归系数p1 2 m 1 2 k是非线性的.非线性回归包括如下几个用于计算回归参数、预 测输出、置信区间以及输出交互图像的函数.1. 非线性最小二乘参数估计对于非线性方程的的系数估计通常采用最小二乘 估计，又叫做非线性最小二乘回归 .在 Matlab 中采用 nlinfit 实现，其调用格式如下：beta=nlinfit(x,y,fun,beta0)返回非线性回归方程系数的最小二乘估计值 .非 线 性 方程 由 fun 给 定 ， fun 为 用 户提 供形 如y = f (beta, x)的函数，beta为待估系数B ； betaO为回归 系数b的初值.beta,r,J=nlinfit(x,y,fun,betaO)返回回归系数beta、残差r、Jacobi矩阵J.2. 最小二乘估计参数的置信区间 求非线性最小二乘估计的系数的置信区间用nlparci 计算。其输入为 nlinfit 函数的输出 beta，r， J。函数具体调用格式如下：ci=nlparci(beta,r,J)返回系数 beta 的 95%置信区间. ci=nlparci(beta,r,J,alpha)返回系数beta的100(1-alpha)%置信区间.3. 最小二乘估计模型的预测输出及其置信区间 非线性最小二乘估计模型的预测输出及其置信区间用 nlpredci 计算.函数具体调用格式如下：ypred,delta=nlpredci(fun,inputs,beta,r,J)返回模型在对应 inputs 处的输出预测值 ypred，给出 95% 的置信区间ypred-delta, ypred+delta； ypred,delta=nlpredci(fun,inputs,beta,r,J,alpha) 返回模型在对应 inputs 处的输出预测值 ypred， 给出 100(1-alpha)% 的置信区间ypred-delta, ypred+delta.4. 非线性拟合和预测的交互图形工具 nlintool 是非线性拟合和预测的交互图形工具 .函 数具体调用格式如下：nlintool(x,y,fun,beta0)返回 x，y 的非线性最小二乘法的曲线拟合图， 并画出 95%的置信区间；nlintool(x,y,fun,beta0,alpha) 给出曲线拟合图及1 00 (1 -alpha) %的置信区间； nlintool(x,y,fun,beta0,alpha,xname,yname) 给出曲线拟合图及1OO(1-alpha)%的置信区间，标 出x，y变量名称.例4在化工生产中获得的氯气的级分y随生产时 间x下降，假定在x 8时，y与x之间有如下形式的非 线性模型： 现收集43组数据，如下所示要求利用该数据求a，b 的值以及它们 95%的置信区间.并画出拟合曲线.此处 设a , b的初值分别为0.30, 0.02.x:8.00,8.00,10.00,10.00,10.00,10.00,12.00,12.00,12.00,14.00 ,14. 00,14. 00,16. 00,16. 00,16. 00,18. 00,18. 00, 20.00,20.00,20.00,20.00,22.00,22.00,24.00,24.00,24.00,26.00 ,26. 00,26. 00,28. 00,28. 00,30. 00,30. 00,30. 00, 32.00,32.00,34.00,36.00,36.00,38.00,38.00,40.00,42.00 y:0.49,0.49,0.48,0.47,0.48,0.47,0.46,0.46,0.45,0.43,0.45,0.43 ,0.43,0.44,0.43,0.43,0.46,0.42,0.42,0. 43,0.41, 0.41,0.40,0.42,0.40,0.40,0.41,0.40,0.41,0.41,0.40,0.40,0.40,0. 38,0.4 1,0.40 ,0.40 ,0.4 1,0.3 8 ,0.40,0.4 0,0.3 9,0.3 9解 首先编写函数fun,表示出需要拟合模型的函 数形式.%拟合模型的函数形式function y=fun(beta0,x)a=beta0(1);b=beta0(2);y=a+(0.49-a)*exp(-b*(x-8); 编写程序对数据进行最小二乘拟合和处理 %输入数据x=8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.0014.0014.0014.0016.0016.0016.0018.0018.0020.0020.0020.0020.0022.0022.0024.0024.0024.0026.0026.0026.0028.0028.0030.0030.0030.0032.0032.0034.00 36.00 36.00 38.00 38.00 40.00 42.00;y=0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43 0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41 0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39;beta0=0.30 0.02; % 系数 初值beta,r,J=nlinfit(x,y,fun,beta0); % %系数最小二乘估计值为 betabetaci=nlparci(beta,r,J); % 求 95%置信区间nlintool(x,y,fun,beta0,0.05,生产时间 x,氯气级分 y); %显示拟合结果和置信区间得如下结果：beta=0.3896 0.1011ci=0.3805 0.39880.0771 0.1251可知a与b的最小二乘估计值分别为0.3896和0.1011，它们 95%的置信区间分别为 0.3805,0.3988和 0.0771,0.1251.所以用最小二乘法对该模型系数进行 估计，最终的模型表达式为最终得到的拟合图如下所示.图 8-8在图 8-8 中，中间的实线为拟合曲线，两侧的虚 线标出其 95%置信区间.四、逐步回归实际问题中影响因变量的因素可能很多，希望从 中挑选出影响显著的自变量来建立回归模型 .这就涉 及到了变量选择的问题 .逐步回归就是一种从众多变 量中有效地选择重要变量的方法.逐步回归的命令是stepwise，它提供了 一个交互式 画面，通过此工具可以自由地选择变量，进行统计分 析，其调用格式为 stepwise（x,y,inmodel,alpha）.其中x表示自变量数据，nxm阶矩阵；y表示因变 量数据，nx1阶矩阵；inmodel表示矩阵的列数指标 （缺省时设定为全部自变量 ）；alpha 表示显著性水平 （缺省时为 0.5）.运行 stepwise 命令时产生三个图形窗口： Stepwise Plot，Stepwise Table，Stepwise History.在 Stepwise Plot 窗口，虚线表示该变量的拟合系 数与 0 无显著差异，实线表示有显著差异，红色线表 示从模型中移去的变量，蓝色线表示存在模型中的变 量，单击一条会改变其状态.Stepwise Table 窗口中列出了一个统计表，包括回 归系数及其置信区间，以及模型的统计量剩余标准差 (RMSE)、相关系数(R-square)、f值、与F对应的概率 P.例5 水泥凝固时放出的热量 y 与水泥中 4 种化学 成分 x1、x2、x3、x4 有关，今测得一组数据如下， 试用逐步回归法确定一个线性模型.序号12345678910111213x17111117113122111110x226295631525571315447406668x3615886917221842398x46052204733226442226341212y78.574.3104.387.695.9109.2102.772.593.1115.983.8113.3109.4解在Matlab中实现逐步回归的具体代码如下:%数据输入：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;x=x1 x2 x3 x4;y=78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1115.9 83.8 113.3 109.4;%逐步回归%先在初始模型中取全部自变量： stepwise(x,y)得图8-9 所示的结果.图8-9从图8-9中可以看出，x3、x4的P值较大，其对 模型影响不显著.若要移去x3、x4两个变量，在Matlab中输入：stepwise(x,y,1,2)得图 8-10 所示的结果.图 8-10移去变量 x3 和 x4 后模型具有显著性，虽然剩余 标准差(RMSE)没有太大的变化，但是统计量f的值明 显增大，因此新的回归模型更好.对变量 y 和 x1、x2 作线性回归，在 Matlab 中输 入：X=ones(13,1) x1 x2;b=regress(y,X)得结果：b二52.57731.46830.6623故最终模型为：y=52.5773+1.4683x1+0.6623x2.讨论题1. 根据经验，在人的身高相等的情况下，血压的 收缩压y与体重xl,年龄x2有关，现在收集了 13个 男子的有关数据，如下表所示，试建立 y 与 x1、x2 的回归模型.X176.091.585.582.579.080.574.579.085.076.582.095.092.5x250202030305060504055404020y1201411241261171251231251321231321551472. 下表为 1980-1991 年间以1987年不变价计算的 美国个人消费支出 y 与美国国内生产总值 x 的数据 (单位： 10亿美元).1) 在直角坐标系下，作 x 与 y 的散点图，并判断y与x是否存在线性相关关系？2) 试求y与x的一元线性回归方程.3）对所得的回归方程作显著性检验o.o5）.4）若 x =4500，试求对应的 y 的点预测和包含00概率为 95%的区间预测.年份yx年份yx19802447.13776.319862969.14404.519812476.93843.119873052.24539.919822503.73760.319883162.44718.619832619.43906.619893223.34838.019842746.14148.519903260.44877.519852865.84279.819913240.84821.03. 某地人事部门研究中学教师的薪金与他们的资历、性别、教育程度、及培训情况等因素之间的关 系，要建立一个数学模型，分析人事策略的合理性， 特别是考察女教师是否受到不公正的待遇，以及她们 的婚姻状况是否会影响收入 .为此，从当地教师中随 机选中 3414 位进行观察，然后从中保留了 90 个观察 对象，得到关键数据如下表.Zx1x2x3x4x5x6x79987000000101514110000102818110100125019110000102819010100102819000000101827000001107230000000129030110000120430010000135231012010120431000100110438000000111841110000112742000000125942110100112742110000112742000100109547000001111352000001146252012010118254110000140454000100118254000000159455112110145966000100123767110100Zx1x2x3x4x5x6x712376701010014967501000014247811010014247901000013479111010013439200000113109400010018141030021101534103000000143010311000014391111101001946114113110221611411411018341141141111416117000001205213911010020871400021112264154002111220115811401129921591151111695162010000179216711010016901730000011827174000001260417511211017201990100001720209000000215920901410018522100100002104213110100185222000000118522220000002210222110000226622301000020272231100001852227000100185223200000119952350000012616245113110232425311010018522570100012054260000000Zx1x2x3x4x5x6x726172841131101948287110000172029001000126043081121101852391101011942319000100202732511000019423261101001720329110100204833700000023343461121111720355000001194235711000021173801100012742387112111274040311211119424061101002266437010000243645301000020674580100002000464112110表中符号介绍如下：Z-月薪（元）；x1-工作时间（月）；x2=1-男性，x2=0-女性；x3=1-男性 或单身女性，x3=0-已婚女 性；x4-学历（取值0-6，值越大表示学历越高）；x5=1- 受雇于重点中学，x5=0-其它；x6=1-受过培训的毕业 生，x6=0-未受过培训的毕业生或受过培训的肄业 生；x7=1 -已两年以上未从事教学工作，x7=0-其它。 注意组合（x2,x3）=（1,1），（0,1），（0,0）的含义.进行变量选择，建立变量 x1-x7 与 Z 的回归模型 （不一定包括每个自变量），说明教师的薪金与哪些变 量的关系密切，是否存在性别和婚姻状况上的差异 . 为了数据处理上的方便，建议对薪金取对数后作为因 变量.4. 热敏电阻器的电阻y与温度x之间的关系模型为试根据下表中的数据给出参数0 ， 0 ， 0的估计.1 2 3xyxyxyxy50347807016370908266110442755286107513720957030115382060236508011540100600512033076519630859744105514712528725. 财政收入预测问题：财政收入与国民收入、工 业总产值、农业总产值、总人口、就业人口、固定资 产投资等因素有关.下表列出了 1952-1981 年的原始数 据，试构造预测模型.年份国民收工业总农业总总人口就业人固定资产财政收入（亿产值（亿产值（万人）口（万投资（亿入（亿元）元）（亿 元）人）元）元）1952598349461574822072944184195358645547558796213648921619547075204916026621832972481955737558529614652232898254195682571555662828230181502681957837798575646532371113928619581028123559865994266002563571959111416815096720726173338444196010791870444662072588038050619617571156434658592559013827119626779644616729525110662301963779104651469172266408526619649431250584704992773612932319651152158163272538286701753931966132219116877454229805212466196712491647697763683081415635219681187156568078534319151273031969137221016888067133225207447197016382747767829923443231256419711780315679085229356203556381972183333657898717735854354658197319783684855892113665237469119741993369689190859373693936551975212142549329242138168462692197620524309955937173883444365719772189492597194974393774547231978247555901058962593985655092219792702606511509754240581564890198027916592119498705418965688261981292768621273100072732804968106. 水泥凝固时放出的热量 y 与水泥中 4 种化学成时间：创作：欧阳法分 x1、x2、x3、x4 有关，今测得一组数据如下，试 用逐步回归法确定一个线性模型.序号123456789101112X71111171131221111X22629563152557131544740662X3615886917221842393X4605220473322644222634124y78.574.3104.387.695.9109.2102.772.593.1115.983.8113.3时间：2021.03.09创作：欧阳法