系统工程学第3章swr

1、城市系统工程学系统分析方法导论第3章逐步回归分析导言¢¢ 分析一个系统，最好能够了解它的结构。¢¢ 传统观点：结构决定功能，功能反作用于结构。¢¢ 新观点：系统结构以及系统与环境的相互作用决定系统的功能。导言为此，需要建立模型，特别是数学模型。好的数学模型有三个层次：特征尺度反映结构的层次。宏观尺度反映环境的影响。微观尺度反映要素的相互作用（参数）。¢¢¢¢1.2.3.导言有时候，我们根本不知道系统的结构。于是利用黑箱（blackbox）原理。输入（input）输出（output）关系。

2、2;¢¢¢¢¢导言输出(output)y输入(input)x黑箱（blackbox）f(x)导言由于不知道系统的结构，我们就探索输入输出关系，然后“模拟”一个结构。这个结构与系统的真实结构可能有区别，但我们关心的是因果关系，其次才是结构。¢¢¢¢导言系统建模的两种类型：理论建模分析法机理型（结构型）。经验建模实验法参数型。第一类模型：参数意义非常明确。第二类模型：参数意义需要我们去探讨。¢¢1.2.¢¢¢¢导言对于数学模型而言，最常用的经验建模方

3、法就是回归分析。¢¢ìì一元线性回归ììïï线性回归íí多元线性回归îîïïïïìì一元非线性回归回归分析ïï非线性回归íííí多元非线性回归îîïïïï其他回归分析ïïïïîî导言¢¢ 多元线性回归不仅仅是普通的线性代

4、数运算，其间存在伪因果和共线性等诸多问题。解决问题的办法已经开发了不少。目前看来效率较高、最受人们青睐的方法是逐步回归法。§6.1逐步回归的基本原理¢¢ 多元逐步回归方法的基本思路在于，自动地从大量的可供选择的变量中选取最重要的变量，据以建立回归分析的或者解释模型。§1逐步回归的基本原理¢¢ 变量选取的根据是自变量对因变量作用程度的大小：保留作用程度大的变量， 剔除作用小的变量。是否选取一个变量， 定量判据之一就是相关系数。§1逐步回归的基本原理¢¢ 假定有m个自变量，1个因变量（用y表示），则全部变量（

5、包括自变量和因变量）之间的相关系数矩阵可以表作éé R11R12 R22"! #!R1m R2m"R1yùùêê R21R2 y úúêêúúúú""R = êêêêRúúRRRêêmy úúm1m2mmêê Ry1úúRy 2RymRyyëëû&#

6、251;§1逐步回归的基本原理¢¢ 根据相关系数定义一个自变量的“贡献” 系数R 2jyP=jRjj¢¢ 按照贡献系数的大小决定一个自变量的去留。§1逐步回归的基本原理¢¢ 根据相关系数定义一个自变量的“贡献” 系数R2jyP=jRjj¢¢ 按照贡献系数的大小决定一个自变量的去留。§1逐步回归的基本原理¢¢ 式中Pj 表示第j个自变量对因变量的贡献系数，Rjy表示第j个自变量与因变量的相关系数，Rjj 表示相关系数矩阵对角线上第j行第j列元素(j=1,2,m)。&#

7、167;§1逐步回归的基本原理¢¢ 由于逐步回归分析是对自变量逐步进行的，每次计算都有一个贡献系数，第l 步计算的贡献系数表示为R (l -1) 2jyP (l )=jR (l -1)jj§1逐步回归的基本原理¢¢ 在逐步回归分析过程中，我们不仅要引入贡献最大的自变量，同时要考虑剔除贡献最小的自变量。因此，变量的存留与否又涉及到另一个统计判据F检验。§1逐步回归的基本原理¢¢ 设定一个显著性水平，查F检验表， 找到F检验的临界值F。当然，也可以根据对系统的认识和建模的需要，由研究者自行设定一个临界值。&#

8、167;1逐步回归的基本原理¢¢ 在第l步计算中，假如第v个自变量的贡献系数最大，数值为= maxP(l) P(l )vh¢¢ 则需要根据F检验来该自变量是否应该被引入模型，这里h为尚且没有被引入模型的变量序号，v为原始变量的序号。§1逐步回归的基本原理¢¢ 计算变量引入的F值n - (l + 1) × P (l)(n - l - 1)P (l)=v=vFinR (l -1) - P (l )R (l -1) - P (l )yyvyyv¢¢ 式中n为样品个数，l为计算步骤数， Pv(l)为第v

9、个变量第l步的贡献系数， Ryy为因变量的自相关系数。§1逐步回归的基本原理¢¢ 如果Fin>F，则在这个显著性水平下，该变量可以被引入模型，否则不要引入。§1逐步回归的基本原理¢¢ 在第l步计算中，如果第v个自变量的贡献系数为P(l )= minP(l) vh¢¢ 则可以根据F检验来应该被剔除。该自变量是否§1逐步回归的基本原理¢¢ 计算变量剔除的F值n - (l + 1) - 1 × P (l)(n - l - 2)P (l)=v=vFoutR (l -1) - P

10、 (l )R (l -1) - P (l )yyvyyv¢¢ 如果Fout F，则在这个显著性水平下， 该变量应该被剔除，否则就要保留。§2 计算方法¢¢ 1）数据准备¢¢ 逐步回归过程的抽象描述过程常常让读者特别是初学者把握不住要领。¢¢ 为了使得讨论具体化，下面借助一个简单的实例说明逐步回归分析的方法。§2 计算方法¢¢ 问题是山东省淄博市旅游业的发展分析， 我们想搞清楚哪些因素影响淄博市的旅游总收入。¢¢ 所能考虑的因素包括：国内游客数量、海外游客数

11、量、第三产业的发展和人均GDP数量（m=4）。¢¢ 从1995年到2004年一共10个年份的数据（n=10）。§2 计算方法§2 计算方法¢¢ 初步分析表明，这些因素都与旅游业总收入具有明确的关系。而且，作为自变量，它们彼此之间也有很强的关系。换言之，如果将这四个变量全部引入模型， 就会导致多重共线性的问题。为了得到简约、可靠的模型，需要借助逐步回归分析技术。§2 计算方法¢¢ 自变量的先后排列顺序对多元逐步回归没有实质性的影响。但是，作为一个教学实例，为了更为有效地说明问题，我们对上表的变量排列顺序稍作

12、调整，结果如下。§2 计算方法§2 计算方法¢¢ 利用上表的数据，容易计算相关系数， 得到矩阵如下。逐步回归计算就是从这种相关系数矩阵出发的。为了表述方便，这个矩阵记为R(0)。§2 计算方法§2 计算方法¢¢ 首先设定F统计量的临界值。取显著性 水平=0.05 ， 我们有m=4 个自变量， n=10个观测值。不妨取回归自由度为4、剩余自由度为n-m-1=10-4-1=5的临界值为我们引入变量的F 值下限， 即取Fc(in)=5.192。§2 计算方法¢¢ 另一方面，假定一个变量被淘汰

13、，则有m=3。我们取显著性水平=0.05 、回归自由度为3 、剩余自由度为n-m- 1=10-3-1=6的F临界值为剔除一个变量的上限，即取Fc(out)=4.757。§2 计算方法¢¢ 上述假定的意思是说， 只要F 值大于5.192，变量引入模型；如果F值小于4.757，则剔除变量不将变量引入模型。当然，这个临界值只是初步的设定，或者是一种临时的设定，后面可以根据研究者对置信度的具体需要调整临界值。§2 计算方法¢¢ 2）第一轮计算¢¢ 这一步的计算可以分解为如下几个步骤。¢¢ （1）计算自变

14、量的贡献系数§2 计算方法¢¢ 利用前述公式计算各个自变量对因变量的贡献系数，结果如下P(1)=R(0) 2/ R(0)= 0.991192= 0.969932= 0.978972/1 = 0.98246/1 = 0.94077/1 = 0.9583811y11P(1)=R(0) 2/ R(0)22 y22P(1)=R(0) 2/ R(0)33 y33P(1) =R(0) / R(0)= 0.962152 /1 = 0.9257441y44§2 计算方法¢¢ （2）找出最大贡献系数及其对应的变量序号¢¢ 显然，P1

15、(1)等于0.98246最大，对应的变量序号v=1。因此，首先考虑引入的变量是国内游客数量x1。§2 计算方法¢¢ （3）计算变量引入的F统计量¢¢ 国内游客数量这个变量是否能被引入模型，还要进行一次F检验。对于我们的问题，n=10，现在计算第l=1步。因此可知(10 - 1 - 1) * 0.98246= 448.035Fin1 - 0.98246§2 计算方法¢¢ 这个数值远远大于我们设定的临界值Fc(in)=5.192，因此变量x1可以被引入模型。§2 计算方法¢¢ 接下来考虑排

16、除贡献系数最小的变量。但是否排除，要视Fout值而定。根据上述计算结果，P4(1)=0.92574最小，其Fout值为87.26(10 - 2 - 1) * 0.92574= 99.725Fout1 - 0.92574¢¢ 这个数值高于剔除变量的F临界值4.757， 因此第一步不能剔除。§2 计算方法¢¢ 作为对比，可以计算出所有变量的F变化值。例如，对于第二个变量“第三产业产值”，变量引入和剔除的F值分别为(10 - 1 - 1) * 0.94077=F= 127.075in1 - 0.94077(10 - 2 - 1) * 0.94077=

17、F= 111.190out1 - 0.94077§2 计算方法¢¢ 其余的Fin值和Fout值都按这种方法计算。结果表明，它们都大于临界值Fc(out)= 4.757。因此这一阶段没有变量被剔除。并且，由于每一轮只能引入一个变量， 其他变量的Fin值即便大于临界值，也 只具有参考价值，不能据此同时引入几 个变量。§2 计算方法§2 计算方法¢¢ （4）相关系数矩阵变换，将R(0)化为R(1)¢¢ 假定第v个变量在第l步被引入，则相关系数矩阵的第v个元素称为主元。矩阵变换是主元进行的。相关系数矩阵的变换公式

18、如下§2 计算方法¢¢ 相关系数矩阵的变换公式如下ììR (l )× R (l )- R (l )/ R (l )( j ¹ v, k ¹ v)( j = v, k ¹ v)( j ¹ v, k = v)( j = v, k = v)jkjvvkvvïïïïR (l )/ R (l )(l +1)jkvvvkR= ííïï- R (l )/ R (l )kvvvïï(l )vv1/ Rî

19、î§2 计算方法¢¢ 式中j、k分别为相关系数矩阵的行列编号。¢¢ 根据这个公式，第一步应该改变非主元所在的行、列的元素（jv, k v ）;¢¢ 第二步改变主要所在的行的元素（v ）;¢¢ 第三步改变主元所在的列的元素（=v ）;¢¢ 第四步改变主元本身（ j=v, k =v ）。j=v,kjv,k§2 计算方法¢¢ 首先变换非主元所在的行和列的元素。我们的主元在第j=1行、第k=1列，故非主元所在的元素为1行、1列以外的元素。举例如下。R(1

20、) = R(0) - R(0) * R(0) / R(0)= 1 - 0.95090 * 0.95090 /1 = 0.095792222211211R(1) = R(0) - R(0) * R(0) / R(0)= 0.95363 - 0.95090 * 0.95261/1 = 0.047802323211311§2 计算方法R(1) = R(0) - R(0) * R(0) / R(0)= 0.95363 - 0.95261* 0.95090 /1 = 0.047803232311211R(1) = R(0) - R(0) * R(0) / R(0)=1 - 0.95261* 0

21、.95261/1 = 0.092543333311311R(1) = R(0) - R(0) * R(0) / R(0)= 0.98102 - 0.96525* 0.95090 /1 = 0.063174242411211¢¢ 其余计算依此类推。§2 计算方法¢¢ 然后改变主元所在行的元素。我们的主元在第j=1行，故改变第1行的元素。例如R(1) = R(0) / R(0)= 0.95090 /1 = 0.95090121211R(1) = R(0) / R(0)= 0.95261/1 = 0.95261131311§2 计算方法&#

22、162;¢ 其余计算依此类推。可以看出，对于本轮计算，主元所在的第1行元素数值不变。§2 计算方法¢¢ 再然后改变主元所在列的元素。我们的主元在第k=1列，故改变第1列的元素。例如R(1) = -R(0) / R(0)= -0.95090 /1 = -0.95090122111R(1) = -R(0) / R(0)= -0.95261/1 = -0.95261133111¢¢ 其余计算依此类推。§2 计算方法¢¢ 可以看出，对于本轮计算，主元所在的第1列元素仅仅改变符号（正号变负号）。§2 计算

23、方法¢¢ 最后改变主元所在的元素。对于本轮计算，主元实际不变：R(1)= 1/ R(0)= 1/1 = 11111¢¢ 于是整个矩阵变换完成工作。这样，我们得到矩阵R(1)。§2 计算方法§2 计算方法¢¢ 在这个矩阵中，第1行最后一列的元素可以用于建立一元线性回归模型。如果我们只打算引入一个关系最密切的变量， 则在数据标准化的情况下，可以建立如下模型y*= 0.99119x*1¢¢ 式中符号带星号表示数据标准化。§2 计算方法¢¢ 3）第二轮计算¢

24、62; （1）计算自变量的贡献系数¢¢ 基于第一个相关系数矩阵的变换结果R(1)计算各个自变量对因变量的贡献系数如下。§2 计算方法P(2)=R(1) 2/ R(1)= 0.991192/1 = 0.9824611y11P(2)=R(1) 2/ R(1)= 0.027412/ 0.09579 = 0.0078522 y22P(2)=R(1) 2/ R(1)= 0.034762/ 0.09254 = 0.0130533 y33P(2)=R(1) 2/ R(1)= 0.005412/ 0.06830 = 0.0004341y44§2 计算方法¢&#

25、162; （2）找出最大贡献系数及其对应的变量序号¢¢ 从上面的计算结果可以看出，不考虑已经被引入模型的第一个变量，在剩余变量中P3(2)=0.01305为最大，对应的变量序号v=3。因此，第二次可能引入的变量是海外游客数量x3。¢¢ 注意，我们的主元已经转移到第三个变量对应的对角线上。§2 计算方法¢¢ （3）计算变量引入和剔除的F统计量¢¢ 海外游客数量能否被引入模型，依然需要借助F检验因此应有。现在计算第l=2步，(10 - 2 - 1) * 0.01305= 20.359Fin0.01754 -

26、0.01305§2 计算方法¢¢ 这个数值大于我们设定的临界值Fc(in)=5.192，因此变量x3可以被引入模型。¢¢ 当我们引入x1的时候，F值为448.035； 现在引入x3，F值在原来的基础上增加了20.359。§2 计算方法¢¢ 在没有被引入也没有被排除的变量中， 找到最小贡献系数，考虑剔除相应的变量。但是是否剔除，依然要视Fout值而定。§2 计算方法¢¢ 根据上面的计算结果，第四个变量“人均GDP”的贡献系数P4(2)=0.00043最小， 其Fout值为(10 - 2

27、- 2) * 0.00043= 0.150Fout0.01754 - 0.00043¢¢ 小于 Fc(out)=4.757,因此，这个变量可以被剔除，不再考虑它的引入。§2 计算方法§2 计算方法¢¢ 特别指出两点：其一，对于已经引入模型的变量，其第l步的F值是采用上一个相关系数矩阵R(l-1)值及其相应的贡献系数计算的，不同之处在于l值增大了（l变成l+1）。其二，非主元的变量对应的F 值是用它们各自的贡献系数计算的，计算公式完全一样。§2 计算方法¢¢ （4）相关系数矩阵变换，将R(1)化为R(2)&

28、#162;¢ 首先变换非主元所在的行和列的元素。我们的主元在第j=3行、第k=3列，故非 主元所在的元素为3行、3列以外的元素。例如R(2) = R(1) - R(1) * R(1) / R(1)1111133133= 1 - 0.95261* (-0.95261) / 0.09254 = 10.80601§2 计算方法R (2) = R (1) - R (1) * R (1) / R (1)1212133233= 0.95090 - 0.95261 * 0.04780 / 0.09254 = 0.45890R(2) = R(1) - R(1) * R(1) / R(1)2

29、121233133= -0.95090 - 0.04780 *(-0.95261) / 0.09254 = -0.45890¢¢ 其余计算依此类推。§2 计算方法¢¢ 然后改变主元所在行的元素。我们的主元现在在第j=3行，故改变第3行的元素。例如R(2) = R(1) / R(1)= -0.95261/ 0.09254 = -10.29388313133R(2) = R(1) / R(1)= 0.04780 / 0.09254 = 0.51648323233¢¢ 其余计算依此类推。§2 计算方法¢

30、2; 再然后改变主元所在列的元素。我们的主元在第k=3列，故改变第1列的元素。例如R(2) = -R(1) / R(1)= -0.95261/ 0.09254 = -10.29388313133R(2) = -R(1) / R(1)= -0.04780 / 0.09254 = -0.51648323233¢¢ 其余计算依此类推。§2 计算方法¢¢ 最后改变主元所在的元素R(2)=1/ R(1)=1/ 0.09254 = 10.806013333¢¢ 这样，我们得到相关矩阵R(2) 。§2 计算方法§2 计

31、算方法¢¢ 假如我们只准备引入两个变量，则基于标准化数据可以建立如下模型y*= 0.63341x* + 0.37558x*13¢¢ 也就是说，在矩阵R(2)的最后一列，第一个元素表示第一个标准化回归系数（对应于x1），第三个元素表示第二个标准化回归系数（对应于x3）。§2 计算方法¢¢ 4）第三轮计算¢¢ 基于第二个相关系数矩阵的变换结果R(2)计算各个自变量对因变量的贡献系数，方法与前面两轮完全一样，接下来的各种计算也是一样。不同的是，每一步计算都是结果进行的。新的相关系数矩阵变换§2 计算方

32、法¢¢ 这一次第二个变量“第三产业产值”的贡献系数最大，为0.00126。但是，Fin值没有达到被引入的标准，而Fout值则达到被剔除的标准。¢¢ 如果我们继续引入新的变量，F值的变化将会很小，或者说F值的增加量很不显著。因此，可以考虑中止引入变量的计算，不再在模型中添加其他变量。§2 计算方法§2 计算方法¢¢ 现在我们已因此，这时可以考虑用两个变量作为剔除标准，三个标量作为引入标准。利用finv函数在Excel中查到临界值如下F(in)0.05,3,10-3-1= 4.757F(out)0.05,2,10-2-

33、1= 4.737§2 计算方法¢¢ 根据调整后的标准，依然是两个变量（x1和x3）引入，两个变量（x2和x4）剔除。如果将显著性水平降低到0.01，则临界值变为F(in)0.01,3,10-3-1= 9.780F(out)0.01,2,10-2-1= 9.547¢¢ 可见在=0.01的显著性水平下，模型引入的变量不变。§2 计算方法¢¢ 到此为止，根据我们的选择标准，变量的引入和剔除计算过程可以结束。§2 计算方法§2 计算方法¢¢ 5）参数估计和模型建设¢

34、2; 接下来我们需要计算模型的回归系数， 然后建立回归分析模型。前面的第一个 相关系数变换矩阵最后一列给出了引入 一个变量时的标准化回归系数：0.99119， 这个数值就是第一个自变量与因变量的 相关系数。§2 计算方法¢¢ 第二个相关系数变换矩阵给出了引入两个变量时的标准化回归系数：0.63341和0.37558，二者之和接近于1。¢¢ 如果我们需要的仅仅是解释模型而非预测模型，则我们的建模工作可以到此为止，得到模型如前所示。§2 计算方法¢¢ 但是，如果我们需开展计算工作，将标准化回归参数转换为非标准化的回归系

35、数。¢¢ 计算过程如下。§2 计算方法¢¢ 第一步，计算原始数据的协方差。¢¢ 利用原始数据表中的数据，不难计算出协方差，并矩阵如下。§2 计算方法¢¢ 第二步，计算非标准化回归系数。¢¢ 有了协方差矩阵，结合前面的相关系数矩阵第二步变换结果R(2)，就可以计算非标准化回归系数。§2 计算方法¢¢ 非标准化回归系数计算公式为c yy× R (l )b=jjycjjb0= y - åbj x jj§2 计算方法

36、2;¢ 这里b0为截距，bj为第j个回归系数， l=2 为计算步骤的编号数，Rjy(l)为相关系数矩阵第l=2步变换结果的最后一列的第 j个元素，cyy为协方差矩阵对角线上的最后一个元素，cjj 为协方差矩阵对角线上对应于第j个被引入变量的元素， 、x j 为y 自变量和因变量的平均值。至于未被引入的变量，回归系数以0计算。§2 计算方法¢¢ 对于上述问题，我们引进了两个变量x1= 国内游客数，x3=海外游客数。可见， j=1对应于国内游客数，j=3对应于海外游客数。§2 计算方法¢¢ 容易算出c yy152.86810&#

37、215; R (2)b× 0.63341 = 0.05842=11yc17969.8340211c yy152.86810× R (2)b× 0.37558 = 5.53314=33 yc0.7043433b0 = y - b1 x1 - b3 x3= 22.768 - 0.05842 * 473.970 - 5.53314 *1.574 = -13.631§2 计算方法¢¢ 于是非标准化数据的回归模型可以表作y = b0 + b1 x1 + b3¢¢ 至此，模型建设工作结束。¢¢ 需要说明的是

38、，上述数据全部是连续计算过程中截取小数点后面5位得到的结果。3§2 计算方法§3 相关矩阵的消元变换¢¢ 上面的相关系数矩阵变换过程非常繁琐， 特别是自变量较多且操作不够熟练的时候，一就会出错。下面介绍一种比较方便的矩阵变换方法：利用Gauss消元法对增广矩阵进行消元变换。§3 相关矩阵的消元变换¢¢ 在相关系数矩阵R(0)旁边增加一个并排的(m+1)×(m+1)=5×56-2-10）。因为第一轮引入了第一个变量x1， 主元在第1行第1列，且主元为1，我们就用第1行的元素对其他行进行消元。§3

39、相关矩阵的消元变换§3 相关矩阵的消元变换¢¢ 消除的办法很简单，用第2行减去第1行乘以第2行的第一个元素0.95090，用第3 行减去第1行乘以第3行的第一个元素0.95261。其余依此类推。这样，很快将第1列的非主元元素全部化为0 。§3 相关矩阵的消元变换§3 相关矩阵的消元变换¢¢ 这个时候，增广矩阵右边的有对应于主元所在列的第1列发生了变化， 其他列没有改变。将这右边对应于x1的第 1列剪贴到左边对应于x1的第1列即主元所在列，就可以完成相关系数矩阵的第一次变换工作。§3 相关矩阵的消元变换¢&

40、#162; 去掉矩阵的“增广”部分，R(0)变成了R(1) 。最终结果与前面的变换结果一样， 但操作过程更为简捷，不容易出错。§3 相关矩阵的消元变换§3 相关矩阵的消元变换¢¢ 第二轮引入了第三个变量x3，主元在矩阵的第3行第3列，且主元不为1。因此首先考虑将第3行第3列化为1，用第3行的所有元素除以第3行第3列的元素0.09254即可。§3 相关矩阵的消元变换¢¢ 接下来用第3行对其余各行消元，将第3 列的其余元素全部变成0，方法同第一步一样。不同之处在于，我们这一次是借助第3行第3列实施变换。§3 相关矩阵的

41、消元变换¢¢ 举例说来，用第1行减去第3行乘以第1行的第三个元素0.95261，用第2行减去第3 行乘以第2行的第三个元素0.04780。其余依此类推。这样，很快将第3列的其他元素全部化为0。§3 相关矩阵的消元变换§3 相关矩阵的消元变换¢¢ 待到第3列非主元元素全部变为0之后我们可以看到，右边矩阵只有对应于主元位置的第3列（总第8列）元素发生了变化，其余元素均无改变。将右边的第3列（对应于右边的x3）剪切并粘贴到左边第3列（对应于左边的x3），即可完成全部变换过程。这样，我们得到相关矩阵R(2)。§3 相关矩阵的消元变换§4 回归结果检验¢¢ 检验过程与多元线性回归分析一样，包括相关系数检验、标准误差检验、F检验、t检验和DW检验。§4 回归结果检验相关系数用于检验模型的拟合优度； 标准误差用于检验模型的F检验用于检验变量之间的线性关系是否可以接受；¢¢¢¢¢¢t检验回归系数是否与0有显著差异；¢¢¢¢DW检验主要用于残差序列是否存在自相关，从而检验