聚类分析 主成分分析和典型相关分析 含matlab程序.pdf

-226-第十二章第十二章 回归分析回归分析 前面我们讲过曲线拟合问题。曲线拟合问题的特点是，根据得到的若干有关变量的一组数据，寻找因变量与（一个或几个）自变量之间的一个函数，使这个函数对那组数据拟合得最好。通常，函数的形式可以由经验、先验知识或对数据的直观观察决定，要作的工作是由数据用最小二乘法计算函数中的待定系数。从计算的角度看，问题似乎已经完全解决了，还有进一步研究的必要吗?从数理统计的观点看，这里涉及的都是随机变量，我们根据一个样本计算出的那些系数，只是它们的一个（点）估计，应该对它们作区间估计或假设检验，如果置信区间太大，甚至包含了零点，那么系数的估计值是没有多大意义的。另外也可以用方差分析方法对模型的误差进行分析，对拟合的优劣给出评价。简单地说，回归分析就是对拟合问题作的统计分析。具体地说，回归分析在一组数据的基础上研究这样几个问题：（i）建立因变量y与自变量mxxx,21L之间的回归模型（经验公式）；（ii）对回归模型的可信度进行检验；（iii）判断每个自变量),2,1(mixiL=对y的影响是否显著；（iv）诊断回归模型是否适合这组数据；（v）利用回归模型对y进行预报或控制。1 数据表的基础知识 1.1 样本空间 在本章中，我们所涉及的均是样本点变量样本点变量类型的数据表。如果有m个变量mxxx,21L，对它们分别进行了n次采样（或观测），得到n个样本点 ),(21imiixxxL，ni,2,1L=则所构成的数据表X可以写成一个mn维的矩阵。=TnTmnijeexXM1)(式中mTimiiiRxxxe=),(21L，ni,2,1L=，ie被称为第i个样本点。样本的均值为 ),(21mxxxxL=，=niijjxnx11，mj,2,1L=样本协方差矩阵及样本相关系数矩阵分别为 TknkkmmijxexensS)()(11)(1=jjiiijmmijsssrR)(其中 -227-=nkjkjikiijxxxxns1)(11 1.2 数据的标准化处理 （1）数据的中心化处理 数据的中心化处理是指平移变换，即 jijijxxx=*，ni,2,1L=；mj,2,1L=该变换可以使样本的均值变为 0，而这样的变换既不改变样本点间的相互位置，也不改变变量间的相关性。但变换后，却常常有许多技术上的便利。（2）数据的无量纲化处理 在实际问题中，不同变量的测量单位往往是不一样的。为了消除变量的量纲效应，使每个变量都具有同等的表现力，数据分析中常用的消量纲的方法，是对不同的变量进行所谓的压缩处理，即使每个变量的方差均变成 1，即 jijijsxx/*=其中=nijijjxxns12)(11。还可以有其它消量纲的方法，如 max/*ijiijijxxx=，min/*ijiijijxxx=jijijxxx/*=，)minmax/(*ijiijiijijxxxx=（3）标准化处理 所谓对数据的标准化处理，是指对数据同时进行中心化压缩处理，即 jjijijsxxx=*，ni,2,1L=，mj,2,1L=。2 一元线性回归 2.1 模型 一元线性回归的模型为 +=xy10，（1）式中，10,为回归系数，是随机误差项，总是假设),0(2N，则随机变量),(210 xNy+。若对y和x分别进行了n次独立观测，得到以下n对观测值 ),(iixy，ni,2,1L=（2）这n对观测值之间的关系符合模型 iixy+=10，ni,2,1L=（3）这里，ix是自变量在第i次观测时的取值，它是一个非随机变量，并且没有测量误差。对应于ix，iy是一个随机变量，它的随机性是由i造成的。),0(2Ni，对于不同的观测，当ji时，i与j是相互独立的。2.2 最小二乘估计方法 -228-2.2.1 最小二乘法 用最小二乘法估计10,的值，即取10,的一组估计值10,，使iy与xyi10+=的误差平方和达到最小。若记 =niiixyQ121010)(),(则 =niiixyQQ121010,10)(),(min),(10 显然0),(10Q，且关于10,可微，则由多元函数存在极值的必要条件得 0)(21100=niiixyQ 0)(21101=niiiixyxQ 整理后，得到下面的方程组 =+=+=niiiniiniiniiniiyxxxyxn1121101110 （4）此方程组称为正规方程组，求解可以得到 =xyxxyyxxniiniii101211)()(（5）称10,为10,的最小二乘估计，其中，yx,分别是ix与iy的样本均值，即 =niixnx11，=niiyny11 关于1的计算公式还有一个更直观的表示方法，即 =niiniiixxyyxx1211)()(-229-=niiniiniiiniiniiyyxxyyxxxxyy121211212)()()()()(xyxyrss=式中=niixxxns122)(11，=niiyyyns122)(11，xyr是x与y的样本相关系数。显然，当iiyx,都是标准化数据时，则有0=x，0=y，1=xs，1=ys。所以，有 00=，xyr=1 回归方程为 xryxy=由上可知，对标准化数据，1可以表示y与x的相关程度。2.2.2 10,的性质 作为一个随机变量，1有以下性质。11是iy的线性组合，它可以写成 =niiiyk11 （6）式中，ik是固定的常量，=niiiixxxxk12)(。证明 事实上 =niininiiiiniiniiixxxxyyxxxxyyxx12111211)()()()()(由于 0)()(1=xnxnyxxynii 所以 ininiiiyxxxx=1121)(2因为1是随机变量),2,1(niyiL=的线性组合，而iy是相互独立、且服从正态分布的，所以，1的抽样分布也服从正态分布。3点估计量1是总体参数1的无偏估计，有 -230-=niiiniiiyEkykEE111)()(=+=+=niiiniiiniixkkxEk1110101)(由于 0)(1121=niniiiniixxxxk 1)()()(1211121=niiniiiininiiiiniixxxxxxxxxxxxk 所以 11)(=E 4估计量1的方差为 =niixx1221)()(Var （7）这是因为 =niiniiniiiniiikkykyk1221221211)(VarVar)(Var 由于 =niiniiniininiiiniixxxxxxxxxxk1212212121212)(1)()(1)(因此，式（7）得证。5 对于总体模型中的参数1，在它的所有线性无偏估计量中，最小二乘估计量1具有最小的方差。记任意一个线性估计量 =niiiyc11 式中ic是任意常数，ic不全为零，ni,2,1L=。要求1是1的无偏估计量，即 111)()(=niiiyEcE 另一方面，由于iixyE10)(+=，所以又可以写成 -231-=+=+=niiiniiniiixccxcE11101101)()(为保证无偏性，ic要满足下列限制 01=niic，01=niiixc 定义iiidkc+=，其中ik是式（6）中的组合系数，id是任意常数，则 +=niiiniiniiniidkdkc11212212212)(Var 由于 =niniiniiiiniiiiniiikxxxxckckdk1121211)()(0)(1)(1)(1212121211=niiniiniiniiniiniiixxxxkxxcxxc 而 )(Var)(1122122=niiniixxk 所以 =+=niid12211)(Var)(Var =niid12的最小值为零，所以，当=niid120时，1的方差最小。但是，只有当0id时，即iikc 时，才有=niid120。所以，最小二乘估计量1在所有无偏估计量中具有最小的方差。同理，可以得出相应于点估计量0的统计性质。对于一元线性正态误差回归模型来说，最小二乘估计量0是iy的线性组合，所以，它的抽样分布也是正态的。它是总体参数0的无偏估计量，即 00)(=E 同样可以证明 )(1)(12220=+=niixxxnVar （8）-232-且0是0的线性无偏的最小方差估计量。2.2.3 其它性质 用最小二乘法拟合的回归方程还有一些值得注意的性质：1残差和为零。残差 iiiyye=，ni,2,1L=由第一个正规方程，得 0)(11101=niiniixye （9）2拟合值iy 的平均值等于观测值iy的平均值，即 yynynniinii=1111 （10）按照第一正规方程，有 0)(110=niiixy 所以 =+=niiniiniiyxy11101)(3当第i次试验的残差以相应的自变量取值为权重时，其加权残差和为零，即 01=niiiex （11）这个结论由第二个正规方程0)(110=niiiixyx即可得出。4当第i次试验的残差以相应的因变量的拟合值为权重时，其加权残差和为零，即 01=iniiey （12）这是因为 0)(1110110=+=+=niiiniiniiiexeex 5最小二乘回归线总是通过观测数据的重心),(yx的。事实上，当自变量取值为x时，由式（5）xy10=所以 yxxyxy=+=+=1110)(2.3 拟合效果分析 当根据一组观测数据得到最小二乘拟合方程后，必须考察一下，是否真的能由所得-233-的模型（iixy10+=）来较好地拟合观测值iy？用iixy10+=能否较好地反映（或者说解释）iy值的取值变化？回归方程的质量如何？误差多大？对这些，都必须予以正确的评估和分析。2.3.1 残差的样本方差 记残差 iiiyye=，ni,2,1L=残差的样本均值为 0)(11=niiiyyne 残差的样本方差为 =niiiniiniiyyneneenMSE121212)(2121)(21 由于有01=niie和01=niiiex的约束，所以，残差平方和有)2(n个自由度。可以证明，在对=niie12除以其自由度)2(n后得到的MSE，是总体回归模型中)(2iVar=的无偏估计量。记 =niieenMSES1221 （13）一个好的拟合方程，其残差总和应越小越好。残差越小，拟合值与观测值越接近，各观测点在拟合直线周围聚集的紧密程度越高，也就是说，拟合方程xy10+=解释y的能力越强。另外，当eS越小时，还说明残差值ie的变异程度越小。由于残差的样本均值为零，所以，其离散范围越小，拟合的模型就越为精确。2.3.2 判定系数（拟合优度）对应于不同的ix值，观测值iy的取值是不同的。建立一元线性回归模型的目的，就是试图以x的线性函数（x10+）来解释y的变异。那么，回归模型xy10+=究竟能以多大的精度来解释y的变异呢？又有多大部分是无法用这个回归方程来解释的呢？nyyy,21L的变异程度可采用样本方差来测度，即 =niiyyns122)(11 根据式（10），拟合值nyyy,21L的均值也是y，其变异程度可以用下式测度 =niiyyns122)(11 下面看一下2s与2 s之间的关系，有 -234-=+=niiiiniiniiiniiyyyyyyyyyy1121212)(2)()()(由于 =+=niiiiniiiiyxxyyyyy110101)()(0)()()(11011011100=+=niiiniiiiniiixyyxyxxy 因此，得到正交分解式为 =+=niiiniiniiyyyyyy121212)()()(（14）记 =niiyySST12)(，这是原始数据iy的总变异平方和，其自由度为1=ndfT；=niiyySSR12)(，这是用拟合直线iixy10+=可解释的变异平方和，其自由度为1=Rdf；=niiiyySSE12)(，这是残差平方和，其的自由度为2=ndfE。所以，有 SSESSRSST+=，ERTdfdfdf+=从上式可以看出，y的变异是由两方面的原因引起的；一是由于x的取值不同，而给y带来的系统性变异；另一个是由除x以外的其它因素的影响。注意到对于一个确定的样本（一组实现的观测值），SST是一个定值。所以，可解释变异SSR越大，则必然有残差SSE越小。这个分解式可同时从两个方面说明拟合方程的优良程度：（1）SSR越大，用回归方程来解释iy变异的部分越大，回归方程对原数据解释得越好；（2）SSE越小，观测值iy绕回归直线越紧密，回归方程对原数据的拟合效果越好。因此，可以定义一个测量标准来说明回归方程对原始数据的拟合程度，这就是所谓的判定系数，有些文献上也称之为拟合优度。判定系数是指可解释的变异占总变异的百分比，用2R表示，有 )1(2SSTSSESSTSSRR=（15）从判定系数的定义看，2R有以下简单性质：（1）102 R；（2）当12=R时，有SSTSSR=，也就是说，此时原数据的总变异完全可以由拟合值的变异来解释，并且残差为零（0=SSE），即拟合点与原数据完全吻合；（3）当02=R时，回归方程完全不能解释原数据的总变异，y的变异完全由与x-235-无关的因素引起，这时SSTSSE=。测定系数时一个很有趣的指标：一方面它可以从数据变异的角度指出可解释的变异占总变异的百分比，从而说明回归直线拟合的优良程度；另一方面，它还可以从相关性的角度，说明原因变量y与拟合变量y 的相关程度，从这个角度看，拟合变量y 与原变量y的相关度越大，拟合直线的优良度就越高。看下面的式子 ),()()()()()(212122112122yyryyyyyyyeyyyyySSTSSRRniiniiniiiiniinii=+=（16）在推导中，注意有 0)(111=niiniiiniiieyyeyye 所以，2R又等于y与拟合变量y 的相关系数平方。还可以证明，2R等于y与自变量x的相关系数，而相关系数的正、负号与回归系数1的符号相同。2.4 显著性检验 2.4.1 回归模型的线性关系检验 在拟合回归方程之前，我们曾假设数据总体是符合线性正态误差模型的，也就是说，y与x之间的关系是线性关系，即 iiixy+=10，),0(2Ni，ni,2,1L=然而，这种假设是否真实，还需进行检验。对于一个实际观测的样本，虽然可以用判定系数2R说明y与y 的相关程度，但是，样本测度指标具有一定的随机因素，还不足以肯定y与x的线性关系。假设y与x之间存在线性关系，则总体模型为 iiixy+=10，ni,2,1L=如果01，则称这个模型为全模型。用最小二乘法拟合全模型，并求出误差平方和为 =niiiyySSE12)(现给出假设0:10=H。如果0H假设成立，则 iiy+=0 这个模型被称为选模型。用最小二乘法拟合这个模型，则有 01=yxy=00 因此，对所有的ni,2,1L=，有 -236-yyi 该拟合模型的误差平方和为 SSTyynii=12)(因此，有 SSTSSE 这就是说，全模型的误差总是小于（或等于）选模型的误差的。其原因是在全模型中有较多的参数，可以更好地拟合数据。假若在某个实际问题中，全模型的误差并不比选模型的误差小很多的话，这说明0H假设成立，即1近似于零。因此，差额)(SSESST 很少时，表明0H成立。若这个差额很大，说明增加了x的线性项后，拟合方程的误差大幅度减少，则应否定0H，认为总体参数1显著不为零。假设检验使用的统计量为 MSEMSRnSSESSRF=)2/(1/式中 1/SSRdfSSRMSRR=)2/(/=nSSEdfSSEMSEE 若假设0:10=H成立，由于SSESSRSST+=，则2/SSE与2/SSR是独立的随机变量，且 )2(/22nSSE，)1(/22SSR 这时 )2,1(=nFMSEMSRF 综上所述，为了检验是否可以用x的线性方程式来解释y，可以进行下面的统计检验。记iy关于ix的总体回归系数为1，则F检验的原假设0H与备则假设1H分别是 0:10=H，0:11H 检验的统计量为 )2,1(=nFMSEMSRF （17）对于检验水平，按自由度（11=n，22=nn）查F分布表，得到拒绝域的临界值)2,1(nF。决策规则为 若)2,1(nFF，则接受0H假设，这时认为1显著为零，无法用x的线性关系式来解释y。若)2,1(nFF，则否定0H，接受1H。这时认为1显著不为零，可以用x的线性关系来解释y。习惯上说，线性回归方程的F检验通过了。需要注意的是，即使F检验通过了，也不说明 iiixy+=10 -237-就是一个恰当的回归模型，事实上，当0H假设被拒绝后，只能说明y与x之间存在显著的线性关系，但很有可能在模型中还包括更多的回归变量，而不仅仅是一个回归变量x。一般地，回归方程的假设检验包括两个方面：一个是对模型的检验，即检验自变量与因变量之间的关系能否用一个线性模型来表示，这是由F检验来完成的；另一个检验是关于回归参数的检验，即当模型检验通过后，还要具体检验每一个自变量对因变量的影响程度是否显著。这就是下面要讨论的t检验。在一元线性分析中，由于自变量的个数只有一个，这两种检验是统一的，它们的效果完全是等价的。但是，在多元线性回归分析中，由于变量的个数只有一个，这两种检验是统一的，它们的效果完全是等价的。但是，在多元线性回归分析中，这两个建议的意义是不同的。从逻辑上说，一般常在F检验通过后，再进一步进行t建议。2.4.2 回归系数的显著性建议 回归参数的建议是考察每一个自变量对因变量的影响是否显著。换句话说，就是要检验每一个总体参数是否显著不为零。首先看对01=的检验。1代表ix变化一个单位对iy的影响程度。对1的检验就是要看这种影响程度与零是否有显著差异。由于)(,(12211=niixxN=niixx1221)()(Var的点估计为 =niixxMSES1212)()(容易证明统计量 )2()(111ntS 事实上，由于 )(Var/)()(Var/)()(11111111SS=其分子)(Var/)(111服从标准正态分布，而分母项有 )2()(/)(/)(Var)(2212212112=nSSEMSExxxxMSESniinii 已知)2(/22nSSE，所以 -238-)2()(111ntS 1的抽样分布清楚后，可以进行1是否显著为零的检验。0:10=H，0:11H 检验统计量为 )(111St=检验统计量1t在01=假设为真时，服从自由度为)2(n的t分布。对于给定的检验水平，则通过t分布表可查到统计量1t的临界值)2(2nt。决策规则是：若)2(21ntt，则接受0H，认为1显著为零；若)2(21ntt，则拒绝0H，认为1显著不为零。当拒绝了0H，认为1显著不为零时，又称1通过了t检验。另一方面，由于 =ntt，则拒绝0H，认为0显著不为零。此外，根据 =,1L，由（20）得=+=niNxxyiiimmii,1),0(2110LL （21）记=nmnmxxxxXLMLMML111111，=nyyYM1 （22）Tn1L=，Tm10L=（20）表为+=),0(2nENXY （23）其中nE为n阶单位矩阵。-240-3.2 参数估计 模型（20）中的参数m,10L仍用最小二乘法估计，即应选取估计值j，使当jj=时，mj,2,1,0L=时，误差平方和=niniimmiiixxyQ1121102)(L （24）达到最小。为此，令 0=jQ，nj,2,1,0L=得 =mjxxxyQxxyQniijimmiijniimmii,2,1,0)(20)(2111011100LLL （25）经整理化为以下正规方程组 =+=+=+=niiimniimmniiimniiimniimniiiniimimniiiniiniiniiniimmniiniiyxxxxxxxyxxxxxxxyxxxn11212211110111112121211110111221110LLL （26）正规方程组的矩阵形式为 YXXXTT=（27）当矩阵X列满秩时，XXT为可逆方阵，（27）式的解为 YXXXTT1)(=（28）将代回原模型得到y的估计值 mmxxy110+=L （29）而这组数据的拟合值为XY=，拟合误差YYe=称为残差，可作为随机误差的估计，而 =niniiiiyyeQ1122)(（30）为残差平方和（或剩余平方和），即)(Q。3.3 统计分析 不加证明地给出以下结果：（i）是的线性无偏最小方差估计。指的是是Y的线性函数；的期望等于；在的线性无偏估计中，的方差最小。-241-（ii）服从正态分布 )(,(12XXNT （31）记nnijTcXX=)()(1。（iii）对残差平方和Q，2)1(=mnEQ，且 )1(22mnQ （32）由此得到2的无偏估计 221=mnQs （33）2s是剩余方差（残差的方差），s称为剩余标准差。（iv）对总平方和=niiyySST12)(进行分解，有 UQSST+=，=niiyyU12)(（34）其中Q是由（24）定义的残差平方和，反映随机误差对y的影响，U称为回归平方和，反映自变量对y的影响。上面的分解中利用了正规方程组。3.4 回归模型的假设检验 因变量y与自变量mxx,1L之间是否存在如模型（20）所示的线性关系是需要检验的，显然，如果所有的|j),1(mjL=都很小，y与mxx,1L的线性关系就不明显，所以可令原假设为 ),1(0:0mjHjL=当0H成立时由分解式（34）定义的QU,满足 )1,()1/(/=mnmFmnQmUF (35)在显著性水平下有上分位数)1,(mnmF，若)1,(mnmFF，接受0H；否则，拒绝。注意注意 接受0H只说明y与mxx,1L的线性关系不明显，可能存在非线性关系，如平方关系。还有一些衡量y与mxx,1L相关程度的指标，如用回归平方和在总平方和中的比值定义复判定系数 SUR=2 （36）2RR=称为复相关系数，R越大，y与mxx,1L相关关系越密切，通常，R大于0.8（或 0.9）才认为相关关系成立。3.5 回归系数的假设检验和区间估计 当上面的0H被拒绝时，j不全为零，但是不排除其中若干个等于零。所以应进-242-一步作如下m个检验),1,0(mjL=：0:)(0=jjH 由（31）式，),(2jjjjcN，jjc是1)(XXT中的第),(jj元素，用2s代替2，由（31）（33）式，当)(0jH成立时)1()1/(/=mntmnQctjjjj （37）对给定的，若)1(|2mnttj，接受)(0jH；否则，拒绝。（37）式也可用于对j作区间估计（mj,1,0L=），在置信水平1下，j的置信区间为)1(,)1(22jjjjjjcsmntcsmnt+（38）其中1=mnQs。3.6 利用回归模型进行预测 当回归模型和系数通过检验后，可由给定的),(0010mxxxL=预测0y，0y是随机的，显然其预测值（点估计）为 mmxxy001100+=L （39）给定可以算出0y的预测区间（区间估计），结果较复杂，但当n较大且ix0接近平均值ix时，0y的预测区间可简化为 ,2020szyszy+（40）其中2z是标准正态分布的上2分位数。对0y的区间估计方法可用于给出已知数据残差iiiyye=),1(niL=的置信区间，ie服从均值为零的正态分布，所以若某个ie的置信区间不包含零点，则认为这个数据是异常的，可予以剔除。4 Matlab 中的回归分析 4.1 多元线性回归 Matlab 统计工具箱用命令 regress 实现多元线性回归，用的方法是最小二乘法，用法是：b=regress(Y,X)其中 Y,X 为按（22）式排列的数据，b 为回归系数估计值m,10L。b,bint,r,rint,stats=regress(Y,X,alpha)这里 Y,X 同上，alpha 为显著性水平（缺省时设定为 0.05），b,bint 为回归系数估计值和它们的置信区间，r,rint 为残差（向量）及其置信区间，stats 是用于检验回归模型的统计量，有四个数值，第一个是2R（见（36）式），第二个是F（见（35）式），第三个-243-是与F对应的概率p，p拒绝0H，回归模型成立，第四个是残差的方差2s（见（33）式）。残差及其置信区间可以用 rcoplot(r,rint)画图。例 1 合金的强度y与其中的碳含量x有比较密切的关系，今从生产中收集了一批数据如下表 1。表 1 x 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 y 42.0 41.5 45.0 45.5 45.0 47.5 49.0 55.0 50.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+=（41）用 regress 和 rcoplot 编程如下：clc,clear x1=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.4722 137.5000 bint=18.6851 36.2594 75.7755 199.2245 stats=0.7985 27.7469 0.0012 4.0883 即4722.270=，5000.1371=，0的置信区间是18.6851,36.2594，1的置信区间是75.7755,199.2245；7985.02=R，7469.27=F，0012.0=p，0883.42=s。可知模型（41）成立。观察命令 rcoplot(r,rint)所画的残差分布，除第 8 个数据外其余残差的置信区间均包含零点，第 8 个点应视为异常点，将其剔除后重新计算，可得 b=30.7820 109.3985 bint=26.2805 35.2834 76.9014 141.8955 stats=0.9188 67.8534 0.0002 0.8797 应该用修改后的这个结果。表 2 1x元 120 140 190 130 155 175 125 145 180 150 2x元 100 110 90 150 210 150 250 270 300 250 y个 102 100 120 77 46 93 26 69 65 85 -244-例 2 某厂生产的一种电器的销售量y与竞争对手的价格1x和本厂的价格2x有关。表 2 是该商品在 10 个城市的销售记录。试根据这些数据建立y与1x和2x的关系式，对得到的模型和系数进行检验。若某市本厂产品售价 160（元），竞争对手售价 170（元），预测商品在该市的销售量。解 分别画出y关于1x和y关于2x的散点图，可以看出y与2x有较明显的线性关系，而y与1x之间的关系则难以确定，我们将作几种尝试，用统计分析决定优劣。设回归模型为 22110 xxy+=（42）编写如下程序：x1=120 140 190 130 155 175 125 145 180 150;x2=100 110 90 150 210 150 250 270 300 250;y=102 100 120 77 46 93 26 69 65 85;x=ones(10,1),x1,x2;b,bint,r,rint,stats=regress(y,x);b,bint,stats 得到 b=66.5176 0.4139 -0.2698 bint=-32.5060 165.5411 -0.2018 1.0296 -0.4611 -0.0785 stats=0.6527 6.5786 0.0247 351.0445 可以看出结果不是太好：0247.0=p，取05.0=时回归模型（42）可用，但取01.0=则模型不能用；6527.02=R较小；10,的置信区间包含了零点。下面将试图用21,xx的二次函数改进它。4.2 多项式回归 如果从数据的散点图上发现y与x呈较明显的二次（或高次）函数关系，或者用线性模型（20）的效果不太好，就可以选用多项式回归。4.2.1 一元多项式回归 一元多项式回归可用命令 polyfit 实现。例 3 将 17 至 29 岁的运动员每两岁一组分为 7 组，每组两人测量其旋转定向能力，以考察年龄对这种运动能力的影响。现得到一组数据如表 3。表 3 年 龄 17 19 21 23 25 27 29 第一人 20.48 25.13 26.15 30.0 26.1 20.3 19.35 第二人 24.35 28.11 26.3 31.4 26.92 25.7 21.3 试建立二者之间的关系。解 数据的散点图明显地呈现两端低中间高的形状，所以应拟合一条二次曲线。选用二次模型 0122axaxay+=（43）编写如下程序：x0=17:2:29;x0=x0,x0;y0=20.48 25.13 26.15 30.0 26.1 20.3 19.35.24.35 28.11 26.3 31.4 26.92 25.7 21.3;-245-p,s=polyfit(x0,y0,2);p 得到 p=-0.2003 8.9782 -72.2150 即2003.02=a，9782.81=a，2150.720=a。上面的s是一个数据结构，用于计算函数值，如 y,delta=polyconf(p,x0,s);y 得到y的拟合值，及预测值y的置信区间半径delta。182022242628101520253035 图1 拟合的交互式画面 用polytool(x0,y0,2)，可以得到一个如1图的交互式画面，在画面中绿色曲线为拟合曲线，它两侧的红线是y的置信区间。你可以用鼠标移动图中的十字线来改变图下方的x值，也可以在窗口内输入，左边就给出y的预测值及其置信区间。通过左下方的Export下拉式菜单，可以输出回归系数等。这个命令的用法与下面将介绍的rstool相似。4.2.2 多元二项式回归 统计工具箱提供了一个作多元二项式回归的命令rstool，它也产生一个交互式画面，并输出有关信息，用法是 rstool(x,y,model,alpha)其中输入数据x,y分别为mn矩阵和n维向量，alpha为显著性水平（缺省时设定为0.05），model由下列4个模型中选择1个（用字符串输入，缺省时设定为线性模型）：linear(线性)：mmxxy+=L110 purequadratic(纯二次)：=+=mjjjjmmxxxy12110L interaction（交叉）：+=mkjkjjkmmxxxxy1110L quadratic(完全二次)：+=mkjkjjkmmxxxxy,1110L 我们再作一遍例2 商品销售量与价格问题，选择纯二次模型，即 2222211122110 xxxxy+=（44）编程如下：x1=120 140 190 130 155 175 125 145 180 150;x2=100 110 90 150 210 150 250 270 300 250;y=102 100 120 77 46 93 26 69 65 85;-246-x=x1 x2;rstool(x,y,purequadratic)140160180-1000100200300150200250 图2 拟合的交互式画面 得到一个如图2所示的交互式画面，左边是1x（=151）固定时的曲线)(1xy及其置信区间，右边是2x（=188）固定时的曲线)(2xy及其置信区间。用鼠标移动图中的十字线，或在图下方窗口内输入，可改变21,xx。图左边给出y的预测值及其置信区间，就用这种画面可以回答例2提出的“若某市本厂产品售价160（元），竞争对手售价170（元），预测该市的销售量”问题。图的左下方有两个下拉式菜单，一个菜单Export用以向Matlab工作区传送数据，包括beta(回归系数)，rmse（剩余标准差），residuals(残差)。模型（41）的回归系数和剩余标准差为 beta=-312.5871 7.2701 -1.7337 -0.0228 0.0037 rmse=16.6436 另一个菜单model用以在上述4个模型中选择，你可以比较一下它们的剩余标准差，会发现以模型（24）的rmse=16.6436最小。注意本例子在Matlab中完全二次模型的形式为 22521421322110 xbxbxxbxbxbby+=（45）5 非线性回归和逐步回归 本节介绍怎样用Matlab统计工具箱实现非线性回归和逐步回归。5.1 非线性回归 非线性回归是指因变量y对回归系数m,1L（而不是自变量）是非线性的。Matlab统计工具箱中的nlinfit，nlparci，nlpredci，nlintool，不仅给出拟合的回归系数，而且可以给出它的置信区间，及预测值和置信区间等。下面通过例题说明这些命令的用法。例4 在研究化学动力学反应过程中，建立了一个反应速度和反应物含量的数学模型，形式为 33221153241xxxxxy+=-247-其中51,L是未知的参数，321,xxx是三种反应物（氢，n戊烷，异构戊烷）的含量，y是反应速度。今测得一组数据如表4，试由此确定参数51,L，并给出其置信区间。51,L的参考值为（0.1，0.05，0.02，1，2）。表4 序号 反应速度y 氢1x n戊烷2x 异构戊烷3x1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120 5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.32 285 300 10 13 3.13 285 190 120 解 首先，以回归系数和自变量为输入变量，将要拟合的模型写成函数文件huaxue.m：function yhat=huaxue(beta,x);yhat=(beta(4)*x(:,2)-x(:,3)/beta(5)./(1+beta(1)*x(:,1)+.beta(2)*x(:,2)+beta(3)*x(:,3);然后，用nlinfit计算回归系数，用nlparci计算回归系数的置信区间，用nlpredci计算预测值及其置信区间，编程如下：clc,clear x0=1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120 5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.32 285 300 10 13 3.13 285 190 120;x=x0(:,3:5);y=x0(:,2);beta=0.1,0.05,0.02,1,2;%回归系数的初值,任意取的 betahat,r,j=nlinfit(x,y,huaxue,beta);%r,j是下面命令用的信息 betaci=nlparci(betahat,r,jacobian,j);betaa=betahat,betaci%回归系数及其置信区间 yhat,delta=nlpredci(huaxue,x,betahat,r,jacobian,j)%y的预测值及其置信区间的半径，置信区间为yhatdelta。用nlintool得到一个交互式画面，左下方的Export可向工作区传送数据，如剩余标-248-准差等。使用命令 nlintool(x,y,huaxue,beta)可看到画面，并传出剩余标准差rmse=0.1933。4.2 逐步回归 实际问