第十二章回归分析

-131-第十二章回归分析前面我们讲过曲线拟合问题。曲线拟合问题的特点是，根据得到的若干有关变量的一组数据，寻找因变量与（一个或几个）自变量之间的一个函数，使这个函数对那组数据拟合得最好。通常，函数的形式可以由经验、先验知识或对数据的直观观察决定，要作的工作是由数据用最小二乘法计算函数中的待定系数。从计算的角度看，问题似乎已经完全解决了，还有进一步研究的必要吗?从数理统计的观点看，这里涉及的都是随机变量，我们根据一个样本计算出的那些系数，只是它们的一个（点）估计，应该对它们作区间估计或假设检验，如果置信区间太大，甚至包含了零点，那么系数的估计值是没有多大意义的。另外也可以用方差分析方法对模型的误差进行分析，对拟合的优劣给出评价。简单地说，回归分析就是对拟合问题作的统计分析。具体地说，回归分析在一组数据的基础上研究这样几个问题：（i）建立因变量y与自变量mxxx,,,21之间的回归模型（经验公式）；（ii）对回归模型的可信度进行检验；（iii）判断每个自变量),,2,1(mixi对y的影响是否显著；（iv）诊断回归模型是否适合这组数据；（v）利用回归模型对y进行预报或控制。§1多元线性回归回归分析中最简单的形式是xy10，yx,均为标量，10,为回归系数，称一元线性回归。它的一个自然推广是x为多元变量，形如mmxxy110(1)2m，或者更一般地)()(110xfxfymm（2）其中),,(1mxxx，),,1(mjfj是已知函数。这里y对回归系数),,,(10m是线性的，称为多元线性回归。不难看出，对自变量x作变量代换，就可将（2）化为（1）的形式，所以下面以（1）为多元线性回归的标准型。1.1模型在回归分析中自变量),,,(21mxxxx是影响因变量y的主要因素，是人们能控制或能观察的，而y还受到随机因素的干扰，可以合理地假设这种干扰服从零均值的正态分布，于是模型记作),0(~2110Nxxymm（3）其中未知。现得到n个独立观测数据),,,(1imiixxy，mnni,,,1，由（3）得niNxxyiiimmii,,1),,0(~2110（4）记-132-nmnmxxxxX111111，nyyY1（5）Tn][1，Tm][10（4）表为),0(~2NXY（6）1.2参数估计用最小二乘法估计模型（3）中的参数。由（4）式这组数据的误差平方和为niTiXYXYQ12)()()(（7）求使)(Q最小，得到的最小二乘估计，记作ˆ，可以推出YXXXTT1)(ˆ（8）将ˆ代回原模型得到y的估计值mmxxyˆˆˆˆ110（9）而这组数据的拟合值为ˆˆXY，拟合误差YYeˆ称为残差，可作为随机误差的估计，而niniiiiyyeQ1122)ˆ(（10）为残差平方和（或剩余平方和），即)ˆ(Q。1.3统计分析不加证明地给出以下结果：（i）ˆ是的线性无偏最小方差估计。指的是ˆ是Y的线性函数；ˆ的期望等于；在的线性无偏估计中，ˆ的方差最小。（ii）ˆ服从正态分布))(,(~ˆ12XXNT（11）（iii）对残差平方和Q，2)1(mnEQ，且)1(~22mnQ（12）由此得到2的无偏估计22ˆ1mnQs（13）2s是剩余方差（残差的方差），s称为剩余标准差。-133-（iv）对Y的样本方差niiyyS12)(进行分解，有UQS，niiyyU12)ˆ(（14）其中Q是由（10）定义的残差平方和，反映随机误差对y的影响，U称为回归平方和，反映自变量对y的影响。1.4回归模型的假设检验因变量y与自变量mxx,,1之间是否存在如模型（1）所示的线性关系是需要检验的，显然，如果所有的|ˆ|j),,1(mj都很小，y与mxx,,1的线性关系就不明显，所以可令原假设为),,1(0:0mjHj当0H成立时由分解式（14）定义的QU,满足)1,(~)1/(/mnmFmnQmUF(15)在显著性水平下有1分位数)1,(1mnmF，若)1,(1mnmFF，接受0H；否则，拒绝。注意拒绝0H只说明y与mxx,,1的线性关系不明显，可能存在非线性关系，如平方关系。还有一些衡量y与mxx,,1相关程度的指标，如用回归平方和在样本方差中的比值定义SUR2（16）]1,0[R称为相关系数，R越大，y与mxx,,1相关关系越密切，通常，R大于0.8（或0.9）才认为相关关系成立。1.5回归系数的假设检验和区间估计当上面的0H被拒绝时，j不全为零，但是不排除其中若干个等于零。所以应进一步作如下m个检验),,1(mj：0:)(0jjH由（11）式，),(~ˆ2jjjjcN，jjc是1)(XXT对角线上的元素，用2s代替2，由（11）~（13）式，当)(0jH成立时)1(~)1/(/ˆmntmnQctjjjj（17）对给定的，若)1(||21mnttj，接受)(0jH；否则，拒绝。(17)式也可用于对j作区间估计（mj,,1,0），在置信水平1下，j的置信区间为-134-])1(ˆ,)1(ˆ[2121jjjjjjcsmntcsmnt（18）其中1mnQs。1.6利用回归模型进行预测当回归模型和系数通过检验后，可由给定的),,(0010mxxx预测0y，0y是随机的，显然其预测值（点估计）为mmxxy001100ˆˆˆˆ（19）给定可以算出0y的预测区间（区间估计），结果较复杂，但当n较大且ix0接近平均值ix时，0y的预测区间可简化为]ˆ,ˆ[210210suysuy（20）其中21u是标准正态分布的21分位数。对0y的区间估计方法可用于给出已知数据残差iiiyyeˆ),,1(ni的置信区间，ie服从均值为零的正态分布，所以若某个ie的置信区间不包含零点，则认为这个数据是异常的，可予以剔除。1.7Matlab实现Matlab统计工具箱用命令regress实现多元线性回归，用的方法是最小二乘法，用法是：b=regress(Y,X)其中Y,X为按（5）式排列的数据，b为回归系数估计值mˆ,,ˆ,ˆ10。[b,bint,r,rint,stats]=regress(Y,X,alpha)这里Y,X同上，alpha为显著性水平（缺省时设定为0.05），b,bint为回归系数估计值和它们的置信区间，r,rint为残差（向量）及其置信区间，stats是用于检验回归模型的统计量，有三个数值，第一个是2R（见（16）式），第二个是F（见（15）式），第3个是与F对应的概率p，p拒绝0H，回归模型成立。残差及其置信区间可以用rcoplot(r,rint)画图。例1合金的强度y与其中的碳含量x有比较密切的关系，今从生产中收集了一批数据如下表：x0.100.110.120.130.140.150.160.170.18y42.041.545.045.545.047.549.055.050.0试先拟合一个函数)(xy，再用回归分析对它进行检验。解先画出散点图：x=0.1:0.01:0.18;y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];plot(x,y,'+')可知y与x大致上为线性关系。设回归模型为xy10（21）-135-用regress和rcoplot编程如下：clc,clearx1=[0.1:0.01:0.18]';y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]';x=[ones(9,1),x1];[b,bint,r,rint,stats]=regress(y,x);b,bint,stats,rcoplot(r,rint)得到b=27.4722137.5000bint=18.685136.259475.7755199.2245stats=0.798527.74690.0012即4722.27ˆ0，6194.140ˆ1，0ˆ的置信区间是[18.6851,36.2594]，1ˆ的置信区间是[75.7755,199.2245]；7985.02R，7469.27F，0012.0p。可知模型（21）成立。观察命令rcoplot(r,rint)所画的残差分布，除第8个数据外其余残差的置信区间均包含零点，第8个点应视为异常点，将其剔除后重新计算，可得b=30.7820109.3985bint=26.280535.283476.9014141.8955stats=0.918867.85340.0002应该用修改后的这个结果。例2某厂生产的一种电器的销售量y与竞争对手的价格1x和本厂的价格2x有关。下表是该商品在10个城市的销售记录。1x元1201401901301551751251451801502x元10011090150210150250270300250Y个10210012077469326696585试根据这些数据建立y与1x和2x的关系式，对得到的模型和系数进行检验。若某市本厂产品售价160（元），竞争对手售价170（元），预测商品在该市的销售量。解分别画出y关于1x和y关于2x的散点图，可以看出y与2x有较明显的线性关系，而y与1x之间的关系则难以确定，我们将作几种尝试，用统计分析决定优劣。设回归模型为22110xxy（22）编写如下程序：x1=[120140190130155175125145180150]';x2=[10011090150210150250270300250]';y=[10210012077469326696585]';x=[ones(10,1),x1,x2];[b,bint,r,rint,stats]=regress(y,x);b,bint,stats得到b=66.51760.4139-0.2698bint=-32.5060165.5411-0.20181.0296-0.4611-0.0785-136-stats=0.65276.57860.0247可以看出结果不是太好：0247.0p，取05.0时回归模型（22）可用，但取01.0则模型不能用；6527.02R较小；10ˆ,ˆ的置信区间包含了零点。下面将试图用21,xx的二次函数改进它。1.8多项式回归如果从数据的散点图上发现y与x呈较明显的二次（或高次）函数关系，或者用线性模型（1）的效果不太好，就可以选用多项式回归。1.8.1一元多项式回归一元多项式回归可用命令polyfit实现。例3将17至29岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄对这种运动能力的影响。现得到一组数据如下表：年龄17192123252729第一人20.4825.1326.1530.026.120.319.35第二人24.3528.1126.331.426.9225.721.3试建立二者之间的关系。解数据的散点图明显地呈现两端低中间高的形状，所以应拟合一条二次曲线。选用二次模型0122axaxay（23）编写如下程序：x0=17:2:29;x0=[x0,x0];y0=[20.4825.1326.1530.026.120.319.35...24.3528.1126.331.426.9225.721.3];[p,s]=polyfit(x0,y0,2);p得到p=-0.20038.9782-72.2150即2003.02a，9782.81a，2150.720a。上面的s