化工数值计算-6

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,回归分析和曲线拟合,生产过程和科学实验中,常用的变量大体可分两,类。一类为确定性变量,另一类为随机变量。确,定性变量是指两个或多个变量之间有确定的关,系.即其中某个变量的每个值,都与一变量的一个,或几个完全确定的值相对应,即它们之间存在着,函数关系:,例如,理想气体的压力,P,与摩尔体积,V,间,存在着确,定的函数关系:,但在实际问题中,由于变量之间的关系比较,复杂,或由于生产或实验过程中不可避免地存在,着误差,使变量之间的关系具有不确定性,也就是,说,某个变量对应的,不是一个或几个确定的值,而,是整个集合的值,这时,变量,x,和,y,间的关系,就称为,相关关系。例如,流体在圆形直管中做湍流时的,情形,通过量纲分析可知,努塞尔特数,Nu,、普兰特,数,Pr,和雷诺数,Re,之间存在着如下相关关系:,这种关系的不确定性,表现为式中,a,和,b,的数,值,在每次测量中不尽相同。不确定的原因,首先,是影响该过程的因素甚多,有些因素至今尚未弄,清;其次是受到实验过程中的偶然因素影响。这,种不确定性关系并不说明上述三个量纲为1的数,群之间无规律可循。相反,通过大量试验,人们发,现,a,和,b,的数值总是围绕着某一定值波动,而且随,着试验次数的增多,a,、,b,的数值趋于稳定。,a,、,b,的稳定值,可作为,a,和,b,的最佳估计值。在一定条,件下,a,=0.023,b,=0.8。由此可见,通过大量试验,是,可以找到隐藏在随机性后面的统计规律性的。,回归分析和曲线拟合是一种处理变量相关关系,的数理统计方法。用它可以寻找隐藏在随机性,后面的统计规律性。,函数与相关是两种不同类型的变量关系,它们之,间并无严格界限。一方面,相关的变量之间,并无,确定的关系,但在一定的条件下,从一定的统计意,义上看,它们之间又可能存在着某种确定的函数,关系。另一方面,由于实际测定的数据中,总存在,着误差,即使是确定性变量,也会出现某些非确定,性结果。,6.1一元线性回归,一元线性回归处理的是两个变量之间的线,性关系。所用的数学模型为一元线性代数模型,其模型方程式是,对这种模型参数的估计,就是根据原始数据点(,x,1,y,1,)、(,x,2,y,2,)、,、(,x,i,y,i,)、,、(,x,n,y,n,),确定式(6-,1)中,a,、,b,的估计值。,在实际体系中,自变量,x,与因变量,y,之间服从线性,关系的情况虽然不多,但在不少情况下,x,、,y,之间,存在着某种函数组合关系。例如,f,1,(,x,y,),f,2,(,x,y,),设两个函数之间服从线性关系,f,1,与,f,2,是不含待定系数的已知函数。若把,f,1,(,x,y,)与,f,2,(,x,y,)分别视为自变量与因变量,则仍可以借用线,性模型去估计其参数值。这种方法称为化直,法。它在化学化工的实际问题中是常见的。例,如单分子基元反应A,B的动力学方程式为,对上式积分得,式中,c,A,-,t,是不呈线性关系的函数。若对方程两边,取对数,上式可化为ln,c,A,-,t,的线性函数:,又例如,按照阿仑尼乌斯定律,反应速率常数,k,与,温度,T,之间不呈线性关系:,但ln,k,与1/,T,则呈线性关系:,这些都是属于可化为线性关系的例子。,一元线性代数模型中的待定参数,a,和,b,称为“估,计值”。之所以称为“估计”值,是因为,a,b,的值,是从实验值中通过数理统计方法确定的。,图6-1一元线性回归,6.1.1方法概述,设有一组实验数据(,x,1,y,1,)、(,x,2,y,2,)、,、(,x,n,y,n,),自变量,x,与因变,量,y,存在着式(6-1)的关系。当,x,取值为,x,i,时,y,的测定值为,y,i,计,算值为,y,i,*,并有,由于参数,a,b,为未知值,故,y,i,*,也是未知值。若将全部实验,数据标绘在,x,-,y,图中(见图6-1),由于各种因素的影响,它们不会,全部落在一条直线上,即,n,个,y,i,不会与,n,个,y,i,*,完全重合,它们将,随机地分布在与,x,i,呈线性关系的,y,i,*,的周围。以,i,表示它们之,间的差值,则有,这里,i,就是误差。它反映了,x,i,使,y,i,偏离直线的各种影响因素的,总和。,现在,要寻找一条最靠近各个数据点的直线,这条直线称为回,归直线。由于回归直线是一切直线中最接近各数据点(,x,i,y,i,),的,用它代表,x,与,y,之间的线性关系,比任何其他直线更为可,靠。究竟如何确定回归曲线中的参数,a,和,b,呢?目前最常用的,方法就是最小二乘法,即残差平方和最小法。,式(6-3)中的误差,i,又称为残差,表示第,i,个数据与回归直线的,偏离程度,则残差平方和,Q,表示全部数据与回归直线的总偏离程度。显然,Q,是,a,和,b,的,函数:,不用残差和,i,的原因是,i,有正有负,相加时可能彼此抵,消,从而不能反映总的偏离程度,而用残差的平方和不会发生,这种现象。,由多元函数的极值理论可知,要使,Q,值最小,a,、,b,必须满足下,列条件:,即得,式(6-6)称为一元线性回归的正规方程组,通过求解该方程组,可得:,式(6-7)中等号右侧的量全部取自原始数据。因此,就可以确,定回归系数,a,和,b,完成参数估计。,为了简化,a,和,b,的表达式,定义:,式中,、,分别为,x,i,和,y,i,的平均值。,x,i,与,之差(,x,i,-,),称为,x,i,的离差;全部,x,i,的离差平方和,称为,x,的,离差平方和,记为,L,xx,:,y,i,与,之差(,y,i,-,),称为,y,i,的离差;全部,y,i,的离差平方和,称为,y,的,离差平方和,记为,L,yy,同理,再令,L,xy,为全部,x,i,的离差与,y,i,的离差乘积的总和:,将以上关系式代入式(6-7),得,由式(6-12)第二式可以看出,回归直线是通过点(,)的。从,力学观点看,(,)相当于,n,个实验点(,x,i,y,i,)的重心,回归直线是,通过重心的。,应当指出:, 残差,i,只用,y,i,-,y,*,i,表示时,表明,y,i,有测量误差,而,x,i,无测量误差;,或表示与,y,i,相比,x,i,的误差很小。因此,测量误差使实验点偏离,回归直线,都表现为,y,i,偏离,y,*,i,。如果,x,i,的误差与,y,i,的误差相比,不可忽略,则两者都必须考虑。这种情况比较复杂,此处不予,介绍。, 求回归方程的计算过程中,不需要事先假定两个变量之间,必须有相关关系。即使是一组杂乱无章的数据,也可以用最,小二乘法绘制一条直线,以表示,x,与,y,的关系。显然,这种情况,下,绘制的直线并无实际意义。,为了判断两个变量间线性关系的优劣程度,引入一个新的指,标,R,称为简单相关系数,它的定义为,R,值不同时,数据点的分布情况如下。,(1),R,= 0,图6-2R = 0的数据点分布,此时,L,xy,=,0,b,= 0。即回归直线平行于,x,轴,y,的变化与,x,无关,表,示数据点的分布是无规则的,如图6-2所示。,但亦有当,R,= 0时,x,与,y,确实存在明显相关性的情况。这种情,形,不能应用线性回归方法,只能用化直线法或曲线拟合法处,理。,(2) 0 |,R,| 0时,b, 0。数据点的,y,值随着,x,增加而增加,这种情况称为,x,与,y,正相关。,R, 0时,b, 0。数据点的,y,值随着,x,增加而减小,这种情况称为,x,与,y,负相关。,R,的绝对值越小,数据点沿回归直线越分散。,图6-301的数据点分布,1的数据点分布,(3) |,R,| = 1,x,与,y,完全相关。全部数据点均落在回归直线上。,若,x,与,y,为非线性相关,但经变量变换后,用回归直线的方法处,理,所求得的回归系数仅对变换后的变量是最佳的,而对原变,量来说则并非最佳,但通常还能令人满意,此时应注意原变量,的残差平方和并非最小。,由以上讨论可知,相关系数,R,的绝对值在0与1之间,而且越接,近于1,其线性关系越密切,那么|,R,|与1接近到什么程度,才能说,明,x,与,y,之间存在线性相关关系呢?要回答这个问题,就要对相,关系数进行显著性检验。由于篇幅所限,有关相关系数的显,著性检验和回归方程的方差分析等问题将不在此讨论。如,有需要,可参考有关数理统计方面的书籍。,6.1.2程序框图,图6-4是一元线性回归的通用计算程序框图。,程序框图中的主要变量:,N,数据点数,X,、,Y,一维数组,用于存放原始数据中的,x,和,y,值,XXL,x,离差平方和,L,xx,YYL,y,离差平方和,L,yy,XYL,x,离差与,y,离差乘积总和,L,xy,A,回归直线截距,a,B,回归直线斜率,b,R,简单相关系数,6.1.3计算实例,例6-1,已知某反应的速率常数,k,与热力学,温度,T,的实验数据如下:,试求,k,-,T,的关系式。,解通常反应速率常数与热力学温度的关系,服从阿仑尼乌斯定律:,k,=,A,e,-,E,/,RT,式中,E,为反应活化能;,T,为热力学温度;,R,为气体通用常数。,上式取对数,且令,y,= ln,k,x,= -1/,T,可得,y,=ln,A,+,x,。,按图6-4,用一元线性回归求得,A,= 1.966,10,9,/min,E,= 79.571kJ/mol。将实,验数据点和利用关系式获得的计算点一起绘制在图6-5中。,源程序:,*,Example 6-1-Eg6-1.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(50), Y(50),Dim XYA As Variant,Cls,N = 5,XYA = Array(363, 0.00718, 373, 0.01376, 383, 0.02701, _,393, 0.05221, 403, 0.09718),K = 0,For I = 1 To N,X(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2,X(I) = -1 / X(I),Y(I) = Log(Y(I),Next I,CallLINEAR1(N, X(), Y(), A, B, R),A = Exp(A): E = B* 8.314,Print A = ; A: Print E = ; E,Print R = ; R,End Sub,*,Sub LINEAR1(N, X(), Y(), A, B, R),*,XT = 0: YT = 0: XX = 0: YY = 0: XY = 0,For I = 1 To N,XT = XT + X(I): YT = YT + Y(I),XX = XX + X(I) * X(I): YY = YY + Y(I) * Y(I),XY = XY + X(I) * Y(I),Next I,XXL = XX - XT * XT / N,YYL = YY - YT * YT / N,XYL = XY - XT * YT / N,B = XYL / XXL,A = (YT - B * XT) / N,R = XYL / Sqr(XXL * YYL),End Sub,执行结果:,A =1966349283.054212 (指前因子),E =79570.97618674007 (活化能),R =.999718315533107(相关系数),源程序中将一元线性回归计算安排在子程序LINEAR1中。,例6-2,某水样BOD测定数据如下:,试确定该水样中有机物生物氧化降解反应的经验速率方程表达式。,解通常水体中有机物生物氧化降解反应的经验速率方程服从下列方程:,式中,BOD、,L,0,为分别为,t,时和初始时刻的生化需氧量(mg/L);,k,为BOD的降,解系数,即耗氧系数/d,-1,托马斯(Thomas)提出将1-e,-,kt,按幂级数展开如下:,与此展开相似的表达式有:,两展开式仅第四项出现微小差别,故可以近似地取,即 BOD=,L,0,整理得,=,+,t,取,y,=,、,x,=,t,按一元线性回归计算,a,=,= 0.25778、,b,=,= 0.01371,从而解得,k,= 6,b,/,a,= 0.31907d,-1,L,0,=,= 182.963mg/L。,源程序(仅列出主程序,回归子程序LINEAR1同例6-1):,*,Example 6-2-Eg6-2.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(50), Y(50),Dim XYA As Variant,Cls,N = 10,XYA = Array(1, 58, 2, 85, 3, 107, 4, 125, 5, 138, _,6, 147, 7, 155, 8, 161, 9, 167, 10, 170),K = 0,For I = 1 To N,X(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2,Y(I) = (X(I) / Y(I) (1 / 3),Next I,Call LINEAR1(N, X(), Y(), A, B, R),Print A = ; A: Print B = ; B,BK = 6* B / A,BL0 = 1 / BK / A / A / A,Print L0 = ; BL0: Print K = ; BK,Print R = ; R,End Sub,执行结果:,A=.2577799110470602 B=1.370838615288584D-02,L0=182.9634248033556 K=.3190718647672257,R=.9900105997750359,此外,塞里奥特(Theriaut)提出了BOD公式的另一种解法:,式中,k,为待估的,k,的近似值;,h,为,k,的允许偏差量。从而有,因,h,甚小,e,-,ht,1-,ht,故上式变为:,式中,a,=,L,0,d,=,L,0,h,x,1,=1-e,-,k,t,x,2,=,t,e,-,k,t,。这样可以首先假设,k, 的初始值为,k,0,利,用实验数据通过二元线性回归(见下节例6-5)确定出,a,和,d,并求出,L,0,=,a,、,h,=,d,/,L,0,;若|,h,|,(误差允许值),则令,k,=,k,0,+,h,重新进行线性回归计算,直至|,h,|,m,+ 1才能求出上式中的,m,+ 1个回归系数。,同样由多元函数的极值理论可知,要使,Q,值最小,a,0,和,a,j,必须满,足下列条件:,式(6-15)经整理可得:,式(6-16)称为多元线性回归模型的正规方程组。它是一个,m,+,1元的线性代数方程组。由于,x,ij,和,y,i,已知,故可求得,m,+1个待定,系数,a,0,a,1,a,m,。,实际计算时,一般作如下处理:先将式(6-16)的第一式写成,然后将式(6-17)代入方程组(6-16)的第2至第,m,+1式,重新组成,一个,m,元线性方程组,其中有,a,1,a,2,a,m,等,m,个待定系数。通过,求解此,m,元线性方程组,获得系数,a,1,a,2,a,m,再代回式(6-17),求得,a,0,。,为简化计算,用,表示第,j,个,x,的平均值,表示,y,的平均值,则,用,L,jk,表示第,j,个,x,离差与第,k,个,x,离差乘积之和,则,用,L,yy,表示,y,离差的平方和,则,用,L,jy,表示第,j,个,x,离差与,y,离差乘积之和,则,将式(6-17)分别代入式(6-16)的第2至,m,+1式,经简化整理可得,如下,m,元线性方程组:,可用主元素消去法求解此式,然后将求得的,a,1,a,2,a,m,代入式,(6-17),求出,a,0,从而完成对多元线性回归模型的参数估计。,多元线性回归的计算中,常用复相关系数衡量数据点之间的,线性优劣。复相关系数定义如下:,式中,U,称为回归平方和:,应当指出,并非所有曲线都可以按这种方法处理。例如,抛物线,就不能通过变量变换把它化为直线。但是如果令,x,1,=,x,x,2,=,x,2,则上式就化成一个包含两个自变量的线性方程,从而将抛物线按二元线性回归计算。对于含多变量的任意,多项式,也可以通过类似的变换,把它们转化成多元线性回归计算。,6.2.2程序框图,图6-6是多元线性回归的通用计算程序框图。,图6-6(a) 多元线性回归的通用计算程序框图(1),图6-6(b) 多元线性回归的通用计算程序框图(2),程序框图中的主要变量:,N,数据点数,M,多元线性模型元数,X,二维数组,用于存放原始数据的,x,值,Y,一维数组,用于存放原始数据的,y,值,YP,值,YYL,L,yy,值,XP,一维数组,用于存放,值,A,二维数组,用于存放,m,元线性方程组的系数,L,jk,B,一维数组,用于存放,m,元线性方程组的常数项,L,jy,C,一维数组,用于存放多元线性模型的系数,a,j,(,j,= 0,1,M,),R,复相关系数,R,0,U,回归平方和,Q,残差平方和,子程序XYF为列主元消去法求解线性方程组的程序,可参见,图5-2和图5-3。,6.2.3计算实例,例6-4,已知某溶液由两种物质组成,c,A,为物,质A的浓度(g/L),c,B,为物质B的浓度(g,/L),为溶液的黏度(mPas)。设数学模型为,=,a,0,+,a,1,c,A,+,a,2,c,B,试根据下列实验数据,确定,a,0,、,a,1,、,a,2,的值。,解按图6-6编写计算源程序。,源程序:,*,Example 6-4-Eg6-4.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(100, 20), Y(100), C(20),Dim XYA As Variant,Cls,N = 15: M = 2,XYA = Array(25.8, 98, 14.5, 15.8, 116, 9.7, 18.1, 104, 11.3, _,13.3, 99, 26, 20.1, 153, 44.7, 10.1, 98, 21, _,17.1, 103, 25.2, 21, 112, 13.7, 23.7, 113, 38.5, _,11.2, 80, 5.8, 10.2, 87, 17.7, 16.4, 138, 40, _,15.9, 98, 17.1, 8, 102, 3, 26, 155, 37.3),K = 0,For I = 1 To N,For J = 1 To M:X(I, J) = XYA(K):K = K + 1: Next J,Y(I) = XYA(K): K = K + 1,Next I,Call LINEAR2(N, M, X(), Y(), C(), R),Print Tab(4); * Results *,For J = 0 To M,Print A(; J; ) = ; Format$(C(J), #.#),Next J,Print R= ; Format$(R, #.#),End Sub, *,Sub LINEAR2(N, M, X(), Y(), C(), R), *,Dim A(20, 20), B(20), XP(20),YP = 0 = = y Average,For I = 1 To N: YP = YP + Y(I): Next I,YP = YP / N,YYL = 0 = = Lyy Average,For I = 1 To N: YYL = YYL + (Y(I) - YP) * (Y(I) - YP): Next I,For J = 1 To M = = Xj Average,XP(J) = 0,For I = 1 To N: XP(J) = XP(J) + X(I, J): Next I,XP(J) = XP(J) / N,Next J,For J = 1 To M = = Ljy,XYL = 0,For I = 1 To N,XYL = XYL + (X(I, J) - XP(J) * (Y(I) - YP),Next I,B(J) = XYL,Next J,For J = 1 To M = = Ljk,For K = 1 To M,XXL = 0,For I = 1 To N,XXL = XXL + (X(I, J) - XP(J) * (X(I, K) - XP(K),Next I,A(J, K) = XXL,Next K,Next J,Call XYF(A(), B(), M, C(),C(0) = YP =a0,For J = 1 To M,C(0) = C(0) - C(J) * XP(J),Next J,U = 0: Q = 0 = R,For I = 1 To N,YI = C(0),For J = 1 To M,YI = YI + C(J) * X(I, J),Next J,U = U + (YI - YP) * (YI - YP),Q = Q + (YI - Y(I) * (YI -Y(I),Next I,R = Sqr(U / (U + Q) -or R = SQR(U/YYL),End Sub,执行结果:,* Results *,A(0)=-27.4324958,A(1)=0.2327103,A(2)=0.4095299,R= 0.7472224,源程序中将多元线性回归安排在子程序LINEAR2中,子程序XYF和XLZY,见例5-3。,例6-5,按塞里奥特方法确定例6-2的BOD,公式中的参数。,解由前叙述可知,塞里奥特方法计算机求解时要应用二元线性回,归,且无常数项,即,a,0,0。对,a,0,0时的多元线性回归,求解回归方程系数的,m,元线性方程组同样具有方程组(6-18)的形式,不过系数矩阵和常数项必须,用下列式子进行计算:,将例6-4的源程序加以修改即可得到本例的计算源程序。,源程序(子程序XYF和XLZY同例6-4):,*,Example 6-5-Eg6-5.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(100, 20), Y(100), C(20), T(100),Dim XYA As Variant,Cls,N = 10: M = 2,XYA = Array(1, 58, 2, 85, 3, 107, 4, 125, 5, 138, _,6, 147, 7, 155, 8, 161, 9, 167, 10, 170),K = 0,For I = 1 To N,T(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2,Next I,E = 0.000001: X0 = 0.5: H = 100 * E,Do While Abs(H) E,For I = 1 To N,X(I, 1) = 1 - Exp(-X0 * T(I),X(I, 2) = T(I) * Exp(-X0 * T(I),Next I,Call LINEAR20(N, M, X(), Y(), C(), R),CL0 = C(1): H = C(2) / C(1),X0 = X0 + H,Loop,Print Tab(4); * Results *,Print L0 = ; CL0: Print K = ; X0: Print R= ; R,End Sub,*,Sub LINEAR20(N, M, X(), Y(), C(), R), *,Dim A(20, 20), B(20), XP(20),YP = 0 = y Average,For I = 1 To N: YP = YP + Y(I): Next I,YP = YP / N,For J = 1 To M = Sjy,XYL = 0,For I = 1 To N,XYL = XYL + X(I, J) * Y(I),Next I,B(J) = XYL,Next J,For J = 1 To M = Sjk,For K = 1 To M,XXL = 0,For I = 1 To N,XXL = XXL + X(I, J) * X(I, K),Next I,A(J, K) = XXL,Next K,Next J,Call XYF(A(), B(), M, C(),U = 0: Q = 0 = R,For I = 1 To N,YI = 0,For J = 1 To M,YI = YI + C(J) * X(I, J),Next J,U = U + (YI - YP) * (YI - YP),Q = Q + (YI - Y(I) * (YI - Y(I),Next I,R = Sqr(U / (U + Q),End Sub,执行结果:,* Results *,L0= 173.428617756483(mg/L),K=0.3306173539389562 (d,-1,),R= 0.9956168,注意按这种方法进行计算时,k,的初始值必须在真实值附近选取,否则将得,出错误的结果。,注意本例源程序中的多元线性回归子程序LINEAR20,只适用于无常数项,的多元线性方程的回归。,6.3剔除可疑数据及其计算程序,6.3.1剔除可疑数据的方法,在线性回归计算中,假定每个测定数据与回归结果之间的误,差均在随机误差允许的范围之内。然而,由于测量误差或过,失误差等多种原因,在一组实验值中,误差往往会超出随机误,差的允许范围。这些数据,称为可疑数据。为保证回归结果,的可靠性,必须剔除这些可疑的数据。,剔除可疑数据,应当有一个科学的标准。这个标准就是统计,判据,属于统计判据的剔除准则有多种。以一元线性回归为,例,其代数模型为,y,=,a,+,bx,。若自变量,x,无测量误差,则,y,的标,准偏差为,式中,n,为原始数据点数;,m,为回归模型中自变量的个数,对一,元线性回归,m,= 1;,i,为残差,即,i,=,y,i,-,a,、,b,是按最小二乘法求出的最佳估计值。根据数理统计分,析,合理的数据,其残差不应超出,的,k,倍。若取,k,= 3,便是常,用的3,准则。据此,可以把残差绝对值超过3,的个别数据(,x,i,y,i,),判为可疑数据而加以剔除。必须指出,3,准,则是以数据点数,n,􀅹,为前提的,当,n,为有限值时,3,判据并不,十分可靠。下面介绍一种广泛采用的判据,即所谓肖维奈特,准则。,按肖维奈特准则,若,n,次等精度测量中,有某个测量值,y,i,其残差,的绝对值超出,k,就可以认为是可疑数据而予以剔除。表6-1,列出了肖维奈特准则中与,n,相对应的,k,值。,表6-1肖维奈特准则的n和k值,n,k,n,k,n,k,n,k,5,1.65,15,2.13,25,2.33,80,2.74,6,1.73,16,2.16,26,2.34,100,2.81,7,1.79,17,2.18,27,2.35,150,2.93,8,1.86,18,2.20,28,2.37,185,3.00,9,1.92,19,2.22,29,2.38,200,3.02,10,1.96,20,2.24,30,2.39,250,3.11,11,2.00,21,2.26,35,2.45,500,3.29,12,2.04,22,2.28,40,2.50,1000,3.48,13,2.07,23,2.30,50,2.58,2000,3.66,14,2.10,24,2.32,60,2.64,5000,3.89,使用这个准则时,可根据回归结果,对全部实验值进行逐级检,查,把属于可疑数据的实验值选出。若发现不止一个可疑数,据,则应把其中残差绝对值最大者剔除,然后重新计算,值。,根据新的,值,再次用肖维奈特准则进行检查。每次只剔除,一个可疑数据,其余数据重新进行回归,直至回归所用的数据,中不再含有可疑数据为止。,6.3.2剔除可疑数据的计算程序框图,图6-7是具有剔除可疑数据功能的一元线性回归通用计算程,序框图,整个计算过程分为输入原始数据、一元线性回归计,算、确定肖维奈特准则的,k,值、确定残差绝对值最大的数据,点、剔除最可疑数据点(即残差绝对值最大的数据点)。,图6-7具有剔除可疑数据功能的一元线性回归通用计算程序框图,程序框图中的主要变量:,N,原始数据点数或剔除可疑数据后的合格数据点数,N,1可疑数据点数,X,一维数组,用于存放原始数据及合格数据中的,x,值,Y,一维数组,用于存放原始数据或合格数据中的,y,值,X,1一维数组,用于存放可疑数据点的,x,值,Y,1一维数组,用于存放可疑数据点的,y,值,A,回归直线截距,B,回归直线斜率,R,简单相关系数,SD,标准偏差,ER,平均相对误差,DALTA,绝对值最大的残差,ID,残差绝对值最大的数据点序号,U,肖维奈特准则的,k,值,子程序LINEAR1A为一元线性回归计算子程序,比例6-1,中的子程序LINEAR1增加了标准偏差和平均相对误差的计,算;子程序RULES为肖维奈特准则中,k,值的计算程序。,采用类似的方法,可以编写能剔除可疑数据的多元线性回归,计算程序框图。,6.3.3计算实例,例6-6,在定温下某反应的活化能(,E,)与压力,(,P,)呈直线关系:,E,=,a,+,bP,由实验测得如下数据:,试用肖维奈特准则剔除其中的可疑数据,并确定,a,和,b,的值,计算相关系,数、标准偏差和平均相对偏差。,解按图6-7的程序框图编写源程序。,源程序:,*,Example 6-6-Eg6-6.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(500), Y(500), X1(500), Y1(500),Dim XYA As Variant,Cls,N = 10,XYA = Array(1, 40.2, 2, 80, 3, 40.9, 4, 41.6, 5, 41.8, _,6, 42.6, 7, 42.6, 8, 70, 9, 43.7, 10, 43.8),K = 0,For I = 1 To N,X(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2,Next I,-,IC = 1: N1 = 0,Do While IC = 1,Call LINEAR1A(N, X(), Y(), A, B, R, SD, ER),Call RULES(N, U),DALTA = 0: ID = 0,For I = 1 To N,DT = Abs(Y(I) - A - B * X(I),If DT DALTA Then,DALTA = DT,ID = I,End If,Next I,If DALTA U * SD Then,INN = 0,For I = 1 To N,If I = ID Then,N1 = N1 + 1,X1(N1) = X(I): Y1(N1) = Y(I),Else,INN = INN + 1,X(INN) = X(I): Y(INN) = Y(I),End If,Next I,N = N - 1: IC = 1,Else,IC = 0,End If,Loop,-,Print N = ; N: Print N1 = ; N1,Print I, X1, Y1,For I = 1 To N1,Print I, X1(I), Y1(I),Next I,Print A = ; A: Print B = ; B: Print R = ; R,Print SD = ; SD: Print ER = ; ER; %,End Sub,*,Sub LINEAR1A(N, X(), Y(), A, B, R, SD, ER),*,XT = 0: YT = 0: XX = 0: YY = 0: XY = 0,For I = 1 To N,XT = XT + X(I): YT = YT + Y(I),XX = XX + X(I) * X(I): YY = YY + Y(I) * Y(I),XY = XY + X(I) * Y(I),Next I,XXL = XX - XT * XT / N,YYL = YY - YT * YT / N,XYL = XY - XT * YT / N,B = XYL / XXL,A = (YT - B * XT) / N,R = XYL / Sqr(XXL * YYL),SD = 0: ER = 0,For I = 1 To N,SD = SD + (Y(I) - A - B * X(I) * (Y(I) - A - B * X(I),ER = ER + Abs(Y(I) - A - B * X(I) / (A + B * X(I),Next I,SD = Sqr(SD / (N - 2): ER = (ER / N) * 100,End Sub,Sub RULES(N, U),Select Case N,执行结果:,N =8,N1 =2,IX1Y1,1 2 80,2 8 70,A =39.80313111545988,B =0.4172211350293531,R =0.9910796148007959,SD =0.1830558892195286,ER =0.3321330871296535%,从源程序的执行结果看,按肖维奈特准则,10个数据点中有两个属于可疑,数据,可疑数据的活化能测定值分别为80和70。,6.4多项式拟合,在化学化工的实验或科研中,经常需要从一,组测定数据,例如从,n,对(,x,i,y,i,)数据,去求自变量,x,和,因变量,y,的近似函数关系式,y,=,p,(,x,)。从图形上,看,这是由给定的,n,个点(,x,i,y,i,)(,i,= 1,2,n,)作曲线,拟合。,在曲线拟合中,多项式拟合问题占特殊的地位。,任何函数在一个比较小的范围内,可以用多项式,任意逼近。因此,在比较复杂的实际问题中,可以,不问,y,与各因素的确切关系,而用多项式拟合进行,分析和计算。,下面以多项式拟合为例,说明曲线拟合的方法和,计算程序。,6.4.1方法概述,设用下列,m,次多项式 :,拟合一组数据(,x,i,y,i,)(,i,= 1,2,n,),即曲线,y,=,f,(,x,)上已给定,n,个,点,用多项式求作该曲线的近似图形。这一问题与前述的插,值问题有类似之处。但插值问题要求近似曲线,y,=,p,(,x,)严格,地通过所给的,n,个点,这一要求将会使近似曲线,y,=,p,(,x,)保留,数据的全部测试点的测量误差。如果个别数据的误差很大,那么插值的效果显然是不够理想的。鉴于这种情况,考虑放,弃严格通过所有结点(,x,i,y,i,)这一要求,而采用别的方法去构造,近似曲线,以尽可能反映所给数据的总趋势。曲线拟合的常,用方法仍然是最小二乘法,即残差平方和最小法。,若以,i,代表结点处的残差,则,残差的平方和为,由于,x,i,与,y,i,为已知值,故,Q,是,a,j,(,j,= 0,1,2,m,)的函数。由,多元函数的极值理论可知,要使,Q,最小,则系数,a,k,必须满足下,式:,即,(,a,j,)=,y,i,上式变换后得:,令,S,l,=,将上两式代入式(6-24)可得,m,+1元线性方程组:,式(6-25)称为正规方程组。这是,m,+1元线性方程组,其系数矩,阵是对称矩阵,可以证明,上述正规方程组有唯一解。所得的,m,次多项式,P,m,(,x,)确能使残差平方和,Q,最小,故,P,m,(,x,)即为所求的,拟合多项式。它是函数,y,=,f,(,x,)的近似表达式,也可用它代替,f,(,x,)作微分、积分计算。此外,必须注意的是上式计算中,1,。,6.4.2程序框图,图6-8是多项式拟合的通用计算程序框图。,图6-8多项式拟合的通用计算程序框图,程序框图中的主要变量:,M,多项式次数,N,原始数据组数,X,、,Y,一维数组,用于存放原始数据的,x,、,y,值,S,一维数组,用于存放正规方程组的,S,l,数值,T,一维数组,用于存放正规方程组各常数项,T,k,的值,B,二维数组,用于存放正规方程组的系数矩阵,A,一维数组,用于存放正规方程组的变量值,即多项式,系数的值,对图6-8需作如下几点说明。, 由于计算机在做乘幂运算时,是将它转化为对数而进行运,算的,所以要考虑当自变量,x,i,取为零时,对S,l,及,T,k,的影响,即应考,虑到,否则计算机将指出发生运算逻辑错误。, 子程序XYF0为列主元消去法求解多元线性方程组的计算,程序,注意与例6-4中子程序XYF稍有不同,本例中为,m,+ 1元,线性方程组,而例6-4中为,m,元线性方程组。, 在实际计算工作中,当多项式的次数,m,7时,方程组容易,出现病态,所拟合的曲线出现较大的起伏。为了避免此种现,象,m,的数值不宜取得太大。,6.4.3计算实例,例6-7,已知某化工厂产品的废品率,y,与该,产品中的某种物质百分含量,x,有,关。实际测定的,数据如下:,试用三次多项式拟合这组数据,并绘制图形。,解按图6-8的程序框图,编写拟合的计算程序。,源程序:(子程序XYF0见例4-6),*,Example 6-7-Eg6-7.frm,*,DefDbl A-H, O-Z,Private Sub Command1_Click(),Dim X(100), Y(100), A(10),Dim XYA As Variant,Cls,N = 16: M = 3,XYA = Array(0.34, 1.3, 0.36, 1.0, 0.37, 0.73, 0.38, 0.9, _,0.39, 0.81, 0.39, 0.7, 0.39, 0.6, 0.4, 0.5, _,0.4, 0.44, 0.41, 0.56, 0.42, 0.3, 0.43, 0.42, _,0.43, 0.35, 0.45, 0.4,