《化工数值计算》

1整理课件回归分析和曲线拟合生产过程和科学实验中,常用的变量大体可分两

类。一类为确定性变量,另一类为随机变量。确

定性变量是指两个或多个变量之间有确定的关

系.即其中某个变量的每个值,都与一变量的一个

或几个完全确定的值相对应,即它们之间存在着2整理课件函数关系:

例如,理想气体的压力P与摩尔体积V间,存在着确

定的函数关系:

3整理课件但在实际问题中,由于变量之间的关系比较

复杂,或由于生产或实验过程中不可避免地存在

着误差,使变量之间的关系具有不确定性,也就是

说,某个变量对应的,不是一个或几个确定的值,而

是整个集合的值,这时,变量x和y间的关系,就称为

相关关系。例如,流体在圆形直管中做湍流时的

情形,通过量纲分析可知,努塞尔特数Nu、普兰特4整理课件数Pr和雷诺数Re之间存在着如下相关关系:

这种关系的不确定性,表现为式中a和b的数

值,在每次测量中不尽相同。不确定的原因,首先

是影响该过程的因素甚多,有些因素至今尚未弄

清;其次是受到实验过程中的偶然因素影响。这5整理课件种不确定性关系并不说明上述三个量纲为1的数

群之间无规律可循。相反,通过大量试验,人们发

现,a和b的数值总是围绕着某一定值波动,而且随

着试验次数的增多,a、b的数值趋于稳定。a、b

的稳定值,可作为a和b的最佳估计值。在一定条

件下,a=0.023,b=0.8。由此可见,通过大量试验,是

可以找到隐藏在随机性后面的统计规律性的。6整理课件回归分析和曲线拟合是一种处理变量相关关系

的数理统计方法。用它可以寻找隐藏在随机性

后面的统计规律性。函数与相关是两种不同类型的变量关系,它们之

间并无严格界限。一方面,相关的变量之间,并无

确定的关系,但在一定的条件下,从一定的统计意7整理课件义上看,它们之间又可能存在着某种确定的函数

关系。另一方面,由于实际测定的数据中,总存在

着误差,即使是确定性变量,也会出现某些非确定

性结果。6.1

一元线性回归

8整理课件一元线性回归处理的是两个变量之间的线

性关系。所用的数学模型为一元线性代数模型,

其模型方程式是

对这种模型参数的估计,就是根据原始数据点(x1,

y1)、(x2,y2)、…、(xi,yi)、…、(xn,yn),确定式(6-9整理课件1)中a、b的估计值。在实际体系中,自变量x与因变量y之间服从线性

关系的情况虽然不多,但在不少情况下,x、y之间

存在着某种函数组合关系。例如f1(x,y),f2(x,y),

设两个函数之间服从线性关系

10整理课件f1与f2是不含待定系数的已知函数。若把f1(x,y)与

f2(x,y)分别视为自变量与因变量,则仍可以借用线

性模型去估计其参数值。这种方法称为化直

法。它在化学化工的实际问题中是常见的。例

如单分子基元反应A

B的动力学方程式为

11整理课件对上式积分得

式中,cA-t是不呈线性关系的函数。若对方程两边

取对数,上式可化为lncA-t的线性函数:

12整理课件又例如,按照阿仑尼乌斯定律,反应速率常数k与

温度T之间不呈线性关系:

但lnk与1/T则呈线性关系:

13整理课件这些都是属于可化为线性关系的例子。一元线性代数模型中的待定参数a和b,称为“估

计值”。之所以称为“估计”值,是因为a,b的值

是从实验值中通过数理统计方法确定的。图6-1一元线性回归14整理课件6.1.1方法概述设有一组实验数据(x1,y1)、(x2,y2)、…、(xn,yn),自变量x与因变

量y存在着式(6-1)的关系。当x取值为xi时,y的测定值为yi,计

算值为yi*,并有

由于参数a,b为未知值,故yi*也是未知值。若将全部实验15整理课件数据标绘在x-y图中(见图6-1),由于各种因素的影响,它们不会

全部落在一条直线上,即n个yi不会与n个yi*完全重合,它们将

随机地分布在与xi呈线性关系的yi*的周围。以δi表示它们之

间的差值,则有

这里δi就是误差。它反映了xi使yi偏离直线的各种影响因素的

总和。16整理课件现在,要寻找一条最靠近各个数据点的直线,这条直线称为回

归直线。由于回归直线是一切直线中最接近各数据点(xi,yi)

的,用它代表x与y之间的线性关系,比任何其他直线更为可

靠。究竟如何确定回归曲线中的参数a和b呢?目前最常用的

方法就是最小二乘法,即残差平方和最小法。式(6-3)中的误差δi又称为残差,表示第i个数据与回归直线的

偏离程度,则残差平方和17整理课件

Q表示全部数据与回归直线的总偏离程度。显然Q是a和b的

函数:

不用残差和

δi的原因是δi有正有负,相加时可能彼此抵消,从而不能反映总的偏离程度,而用残差的平方和不会发生18整理课件这种现象。由多元函数的极值理论可知,要使Q值最小,a、b必须满足下

列条件:

19整理课件即得

式(6-6)称为一元线性回归的正规方程组,通过求解该方程组,

可得:20整理课件

式(6-7)中等号右侧的量全部取自原始数据。因此,就可以确

定回归系数a和b,完成参数估计。21整理课件为了简化a和b的表达式,定义:

式中,

、

分别为xi和yi的平均值。xi与

之差(xi-

),称为xi的离差;全部xi的离差平方和,称为x的22整理课件离差平方和,记为Lxx:

23整理课件yi与

之差(yi-

),称为yi的离差;全部yi的离差平方和,称为y的离差平方和,记为Lyy,同理

再令Lxy为全部xi的离差与yi的离差乘积的总和:24整理课件

将以上关系式代入式(6-7),得25整理课件

由式(6-12)第二式可以看出,回归直线是通过点(

,

)的。从力学观点看,(

,

)相当于n个实验点(xi,yi)的重心,回归直线是通过重心的。应当指出:①残差δi只用yi-y*i表示时,表明yi有测量误差,而xi无测量误差;26整理课件或表示与yi相比,xi的误差很小。因此,测量误差使实验点偏离

回归直线,都表现为yi偏离y*i。如果xi的误差与yi的误差相比,

不可忽略,则两者都必须考虑。这种情况比较复杂,此处不予

介绍。②求回归方程的计算过程中,不需要事先假定两个变量之间

必须有相关关系。即使是一组杂乱无章的数据,也可以用最

小二乘法绘制一条直线,以表示x与y的关系。显然,这种情况

下,绘制的直线并无实际意义。27整理课件为了判断两个变量间线性关系的优劣程度,引入一个新的指

标R,称为简单相关系数,它的定义为

R值不同时,数据点的分布情况如下。(1)R=028整理课件图6-2R=0的数据点分布29整理课件此时Lxy

=

0,b=0。即回归直线平行于x轴,y的变化与x无关,表

示数据点的分布是无规则的,如图6-2所示。但亦有当R=0时,x与y确实存在明显相关性的情况。这种情

形,不能应用线性回归方法,只能用化直线法或曲线拟合法处

理。(2)0<|R|<130整理课件绝大多数属于这种情况,它表示x与y之间存在着一定的线性

关系,如图6-3所示。R>0时,b>0。数据点的y值随着x增加而增加,这种情况称为

x与y正相关。R<0时,b<0。数据点的y值随着x增加而减小,这种情况称为

x与y负相关。31整理课件R的绝对值越小,数据点沿回归直线越分散。图6-30<<1的数据点分布32整理课件

<1的数据点分布(3)|R|=1x与y完全相关。全部数据点均落在回归直线上。若x与y为非线性相关,但经变量变换后,用回归直线的方法处

理,所求得的回归系数仅对变换后的变量是最佳的,而对原变

量来说则并非最佳,但通常还能令人满意,此时应注意原变量33整理课件的残差平方和并非最小。由以上讨论可知,相关系数R的绝对值在0与1之间,而且越接

近于1,其线性关系越密切,那么|R|与1接近到什么程度,才能说

明x与y之间存在线性相关关系呢?要回答这个问题,就要对相

关系数进行显著性检验。由于篇幅所限,有关相关系数的显

著性检验和回归方程的方差分析等问题将不在此讨论。如

有需要,可参考有关数理统计方面的书籍。34整理课件6.1.2程序框图图6-4是一元线性回归的通用计算程序框图。程序框图中的主要变量:

N数据点数X、Y一维数组,用于存放原始数据中的x和y值35整理课件XXL

x离差平方和LxxYYL

y离差平方和LyyXYL

x离差与y离差乘积总和LxyA回归直线截距aB回归直线斜率b36整理课件R简单相关系数6.1.3计算实例37整理课件例6-1已知某反应的速率常数k与热力学

温度T的实验数据如下:试求k-T的关系式。解通常反应速率常数与热力学温度的关系,服从阿仑尼乌斯定律:k=Ae-E/

RT式中,E为反应活化能;T为热力学温度;R为气体通用常数。上式取对数,且令y=lnk,x=-1/T,可得y=lnA+

x。按图6-4,用一元线性回归求得A=1.966×109/min,E=79.571kJ/mol。将实

验数据点和利用关系式获得的计算点一起绘制在图6-5中。38整理课件

39整理课件源程序:’********************************************’Example6-1------Eg6-1.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(50),Y(50)DimXYAAsVariantClsN=5XYA=Array(363,0.00718,373,0.01376,383,0.02701,_393,0.05221,403,0.09718)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2X(I)=-1/X(I)Y(I)=Log(Y(I))NextICallLINEAR1(N,X(),Y(),A,B,R)A=Exp(A):E=B*8.314Print"A=";A:Print"E=";EPrint"R=";REndSub’******************SubLINEAR1(N,X(),Y(),A,B,R)’******************XT=0:YT=0:XX=0:YY=0:XY=0ForI=1ToNXT=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)NextIXXL=XX-XT*XT/NYYL=YY-YT*YT/NXYL=XY-XT*YT/NB=XYL/XXLA=(YT-B*XT)/NR=XYL/Sqr(XXL*YYL)EndSub执行结果:A=1966349283.054212(指前因子)E=79570.97618674007(活化能)R=.999718315533107(相关系数)源程序中将一元线性回归计算安排在子程序LINEAR1中。40整理课件例6-2某水样BOD测定数据如下:41整理课件试确定该水样中有机物生物氧化降解反应的经验速率方程表达式。解通常水体中有机物生物氧化降解反应的经验速率方程服从下列方程:

式中,BOD、L0为分别为t时和初始时刻的生化需氧量(mg/L);k为BOD的降

解系数,即耗氧系数/d-1托马斯(Thomas)提出将1-e-kt按幂级数展开如下:

与此展开相似的表达式有:

两展开式仅第四项出现微小差别,故可以近似地取

即BOD=L0

整理得

=

+

t取y=

、x=t按一元线性回归计算a=

=0.25778、b=

=0.01371,从而解得k=6b/a=0.31907d-1,L0=

=182.963mg/L。源程序(仅列出主程序,回归子程序LINEAR1同例6-1):’********************************************’Example6-2------Eg6-2.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(50),Y(50)DimXYAAsVariantClsN=10XYA=Array(1,58,2,85,3,107,4,125,5,138,_6,147,7,155,8,161,9,167,10,170)K=0ForI=1ToN42整理课件X(I)=XYA(K):Y(I)=XYA(K+1):K=K+2Y(I)=(X(I)/Y(I))^(1/3)NextICallLINEAR1(N,X(),Y(),A,B,R)Print"A=";A:Print"B=";BBK=6*B/ABL0=1/BK/A/A/APrint"L0=";BL0:Print"K=";BKPrint"R=";REndSub执行结果:A=.2577799110470602B=1.370838615288584D-02L0=182.9634248033556K=.3190718647672257R=.9900105997750359此外,塞里奥特(Theriaut)提出了BOD公式的另一种解法:

式中,k'为待估的k的近似值;h为k的允许偏差量。从而有

因h甚小,e-ht≈1-ht,故上式变为:

式中,a=L0,d=L0h,x1=1-e-k't,x2=te-k't。这样可以首先假设k'的初始值为k'0,利

用实验数据通过二元线性回归(见下节例6-5)确定出a和d,并求出L0=a、h

=d/L0;若|h|≥

ε(误差允许值),则令k'

=k'0+h,重新进行线性回归计算,直至|h

|<ε。43整理课件例6-3尿中胆色素经处理后,在550nm处有

很强的吸光性,现测得配制好的不同

胆色素质量浓度的标准溶液的吸光率数据如下。假定标定曲线用y=bx表示,试

计算该方程中参数b的值。解就一元线性回归计算而言,回归直线方程y=bx对有些实际问题更加

合理,本例就是这样。此时,对式(6-6),a=0,则

将例6-1源程序加以修改就可用于本例计算。44整理课件源程序:’********************************************’Example6-3------Eg6-3.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(50),Y(50)DimXYAAsVariantClsN=9XYA=Array(0.5,0.039,0.75,0.061,1,0.087,1.25,0.107,1.5,_0.119,1.75,0.163,2,0.179,2.25,0.194,2.5,0.213)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2NextICallLINEAR10(N,X(),Y(),B,R)Print"B=";B:Print"R=";REndSub’******************SubLINEAR10(N,X(),Y(),B,R)’******************XT=0:YT=0:XX=0:YY=0:XY=0ForI=1ToNXT=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)NextIXXL=XX-XT*XT/NYYL=YY-YT*YT/NXYL=XY-XT*YT/NB=XY/XXR=XYL/Sqr(XXL*YYL)EndSub执行结果:B=0.08653125R=0.9949028678863266源程序中将一元线性回归计算安排在子程序LINEAR10中,该子程序适用

于规定截距aº0的场合。45整理课件6.2

多元线性回归

一元线性代数模型中,只有一个自变量。若

有多个因素影响体系的性质时,必须考虑因变量y

与多个自变量xl,x2,…,xn之间的关系。例如,化学

反应速率要受到温度、压力和浓度的影响。在46整理课件气相反应动力学中,反应动力学方程可表示为

式中,r为反应速率,pA、pB、pC分别为反应物A、

B、C组分的分压;a、b、c是方程式中的待定指

数;k为反应速率常数。47整理课件若将上式取对数得到

再令y=lnr,d=lnk,x1=lnpA,x2=lnpB,x3=lnpC,则

得

48整理课件可见该式具有多元线性方程式的特征,a、b、c、

d为系数,x1、x2、x3为自变量。多元线性方程的普遍式为

它是一个含有m+1个系数的m元线性方程式,下49整理课件面介绍多元线性回归的最小二乘法。6.2.1方法概述设x取值为xi1,xi2,…,xim时,实验测定的y值为yi(i=1,2,…,n),由于

测定值yi存在着误差,所以会偏离线性关系。现在要寻找一组

aj的估计值以构成回归方程。确定aj的原则,仍然是使yi的实

验值与回归方程计算值的残差平方和最小,即使50整理课件

最小。式中i表示实验点序号(i=1,2,…,n);j表示自变量标号(j

=1,2,…,m);自变量xij为第j个自变量的第i次测定值。此外注

意仅n>m+1才能求出上式中的m+1个回归系数。同样由多元函数的极值理论可知,要使Q值最小,a0和aj必须满51整理课件足下列条件:

式(6-15)经整理可得:52整理课件

式(6-16)称为多元线性回归模型的正规方程组。它是一个m+

1元的线性代数方程组。由于xij和yi已知,故可求得m+1个待定53整理课件系数a0,a1,…,am。实际计算时,一般作如下处理:先将式(6-16)的第一式写成

然后将式(6-17)代入方程组(6-16)的第2至第m+1式,重新组成

一个m元线性方程组,其中有a1,a2,…,am等m个待定系数。通过

求解此m元线性方程组,获得系数a1,a2,…,am,再代回式(6-17),54整理课件求得a0。为简化计算,用

表示第j个x的平均值,

表示y的平均值,则

用Ljk表示第j个x离差与第k个x离差乘积之和,则55整理课件

用Lyy表示y离差的平方和,则

56整理课件用Ljy表示第j个x离差与y离差乘积之和,则

将式(6-17)分别代入式(6-16)的第2至m+1式,经简化整理可得

如下m元线性方程组:57整理课件

可用主元素消去法求解此式,然后将求得的a1,a2,…,am代入式

(6-17),求出a0,从而完成对多元线性回归模型的参数估计。多元线性回归的计算中,常用复相关系数衡量数据点之间的

线性优劣。复相关系数定义如下:58整理课件

式中,U称为回归平方和:

59整理课件应当指出,并非所有曲线都可以按这种方法处理。例如

抛物线

就不能通过变量变换把它化为直线。但是如果令x1=x,x2=x

2,则上式就化成一个包含两个自变量的线性方程

60整理课件从而将抛物线按二元线性回归计算。对于含多变量的任意

多项式

也可以通过类似的变换,把它们转化成多元线性回归计算。6.2.2程序框图图6-6是多元线性回归的通用计算程序框图。61整理课件图6-6(a)多元线性回归的通用计算程序框图(1)62整理课件图6-6(b)多元线性回归的通用计算程序框图(2)63整理课件程序框图中的主要变量:

N数据点数M多元线性模型元数X二维数组,用于存放原始数据的x值Y一维数组,用于存放原始数据的y值64整理课件YP

值YYL

Lyy值XP一维数组,用于存放

值A二维数组,用于存放m元线性方程组的系数LjkB一维数组,用于存放m元线性方程组的常数项Ljy65整理课件C一维数组,用于存放多元线性模型的系数aj(j=0,1,

…,M)R复相关系数R0U回归平方和Q残差平方和子程序XYF为列主元消去法求解线性方程组的程序,可参见66整理课件图5-2和图5-3。6.2.3计算实例67整理课件例6-4已知某溶液由两种物质组成,cA为物

质A的浓度(g/L),cB为物质B的浓度(g

/L),68整理课件μ为溶液的黏度(mPa·s)。设数学模型为μ=a0+a1cA+a2cB试根据下列实验数据,确定a0、a1、a2的值。解按图6-6编写计算源程序。源程序:’********************************************’Example6-4------Eg6-4.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(100,20),Y(100),C(20)DimXYAAsVariantClsN=15:M=2XYA=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=0ForI=1ToNForJ=1ToM:X(I,J)=XYA(K):K=K+1:NextJY(I)=XYA(K):K=K+1NextICallLINEAR2(N,M,X(),Y(),C(),R)PrintTab(4);"***Results***"ForJ=0ToMPrint"A(";J;")=";Format$(C(J),"#####.#######")69整理课件NextJPrint"R=";Format$(R,"###.#######")EndSub’***************************SubLINEAR2(N,M,X(),Y(),C(),R)’***************************DimA(20,20),B(20),XP(20)YP=0'==yAverageForI=1ToN:YP=YP+Y(I):NextIYP=YP/NYYL=0'==LyyAverageForI=1ToN:YYL=YYL+(Y(I)-YP)*(Y(I)-YP):NextIForJ=1ToM'==XjAverageXP(J)=0ForI=1ToN:XP(J)=XP(J)+X(I,J):NextIXP(J)=XP(J)/NNextJForJ=1ToM'==LjyXYL=0ForI=1ToNXYL=XYL+(X(I,J)-XP(J))*(Y(I)-YP)NextIB(J)=XYLNextJForJ=1ToM'==LjkForK=1ToMXXL=0ForI=1ToNXXL=XXL+(X(I,J)-XP(J))*(X(I,K)-XP(K))NextIA(J,K)=XXLNextKNextJCallXYF(A(),B(),M,C())C(0)=YP'==a0ForJ=1ToMC(0)=C(0)-C(J)*XP(J)NextJU=0:Q=0'==RForI=1ToNYI=C(0)ForJ=1ToMYI=YI+C(J)*X(I,J)NextJU=U+(YI-YP)*(YI-YP)Q=Q+(YI-Y(I))*(YI-Y(I))NextIR=Sqr(U/(U+Q))'--orR=SQR(U/YYL)EndSub70整理课件执行结果:***Results***A(0)=-27.4324958A(1)=0.2327103A(2)=0.4095299R=0.7472224源程序中将多元线性回归安排在子程序LINEAR2中,子程序XYF和XLZY

见例5-3。71整理课件例6-5按塞里奥特方法确定例6-2的BOD

公式中的参数。72整理课件解由前叙述可知,塞里奥特方法计算机求解时要应用二元线性回

归,且无常数项,即a0º0。对a0º0时的多元线性回归,求解回归方程系数的m

元线性方程组同样具有方程组(6-18)的形式,不过系数矩阵和常数项必须

用下列式子进行计算:

将例6-4的源程序加以修改即可得到本例的计算源程序。源程序(子程序XYF和XLZY同例6-4):’********************************************’Example6-5------Eg6-5.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(100,20),Y(100),C(20),T(100)DimXYAAsVariantClsN=10:M=2XYA=Array(1,58,2,85,3,107,4,125,5,138,_6,147,7,155,8,161,9,167,10,170)K=0ForI=1ToNT(I)=XYA(K):Y(I)=XYA(K+1):K=K+2NextIE=0.000001:X0=0.5:H=100*EDoWhileAbs(H)>EForI=1ToNX(I,1)=1-Exp(-X0*T(I))X(I,2)=T(I)*Exp(-X0*T(I))NextICallLINEAR20(N,M,X(),Y(),C(),R)CL0=C(1):H=C(2)/C(1)X0=X0+HLoopPrintTab(4);"***Results***"Print"L0=";CL0:Print"K=";X0:Print"R=";REndSub73整理课件’***************************SubLINEAR20(N,M,X(),Y(),C(),R)’***************************DimA(20,20),B(20),XP(20)YP=0'==yAverageForI=1ToN:YP=YP+Y(I):NextIYP=YP/NForJ=1ToM'==SjyXYL=0ForI=1ToNXYL=XYL+X(I,J)*Y(I)NextIB(J)=XYLNextJForJ=1ToM'==SjkForK=1ToMXXL=0ForI=1ToNXXL=XXL+X(I,J)*X(I,K)NextIA(J,K)=XXLNextKNextJCallXYF(A(),B(),M,C())U=0:Q=0'==RForI=1ToNYI=0ForJ=1ToMYI=YI+C(J)*X(I,J)NextJU=U+(YI-YP)*(YI-YP)Q=Q+(YI-Y(I))*(YI-Y(I))NextIR=Sqr(U/(U+Q))EndSub执行结果:***Results***L0=173.428617756483(mg/L)K=0.3306173539389562(d-1)R=0.9956168注意按这种方法进行计算时,k'的初始值必须在真实值附近选取,否则将得

出错误的结果。注意本例源程序中的多元线性回归子程序LINEAR20,只适用于无常数项

的多元线性方程的回归。74整理课件6.3

剔除可疑数据及其计算程序

6.3.1剔除可疑数据的方法在线性回归计算中,假定每个测定数据与回归结果之间的误

差均在随机误差允许的范围之内。然而,由于测量误差或过

失误差等多种原因,在一组实验值中,误差往往会超出随机误75整理课件差的允许范围。这些数据,称为可疑数据。为保证回归结果

的可靠性,必须剔除这些可疑的数据。剔除可疑数据,应当有一个科学的标准。这个标准就是统计

判据,属于统计判据的剔除准则有多种。以一元线性回归为

例,其代数模型为y=a+bx。若自变量x无测量误差,则y的标

准偏差为76整理课件

式中,n为原始数据点数;m为回归模型中自变量的个数,对一

元线性回归m=1;δi为残差,即

δi=yi-

a、b是按最小二乘法求出的最佳估计值。根据数理统计分

析,合理的数据,其残差不应超出σ的k倍。若取k=3,便是常77整理课件用的3σ准则。据此,可以把残差绝对值超过3σ

的个别数据(xi,yi),判为可疑数据而加以剔除。必须指出,3σ准

则是以数据点数n

¥为前提的,当n为有限值时,3σ判据并不

十分可靠。下面介绍一种广泛采用的判据,即所谓肖维奈特

准则。按肖维奈特准则,若n次等精度测量中,有某个测量值yi,其残差

的绝对值超出kσ,就可以认为是可疑数据而予以剔除。表6-1

列出了肖维奈特准则中与n相对应的k值。78整理课件表6-1肖维奈特准则的n和k值nk

nk

nk

nk51.65152.13252.33802.7461.73162.16262.341002.8171.79172.18272.351502.9381.86182.20282.371853.0091.92192.22292.382003.02101.96202.24302.392503.11112.00212.26352.455003.29122.04222.28402.5010003.48132.07232.30502.5820003.66142.10242.32602.6450003.8979整理课件使用这个准则时,可根据回归结果,对全部实验值进行逐级检

查,把属于可疑数据的实验值选出。若发现不止一个可疑数

据,则应把其中残差绝对值最大者剔除,然后重新计算σ值。

根据新的σ值,再次用肖维奈特准则进行检查。每次只剔除

一个可疑数据,其余数据重新进行回归,直至回归所用的数据

中不再含有可疑数据为止。6.3.2剔除可疑数据的计算程序框图80整理课件图6-7是具有剔除可疑数据功能的一元线性回归通用计算程

序框图,整个计算过程分为输入原始数据、一元线性回归计

算、确定肖维奈特准则的k值、确定残差绝对值最大的数据

点、剔除最可疑数据点(即残差绝对值最大的数据点)。图6-7具有剔除可疑数据功能的一元线性回归通用计算程序框图81整理课件程序框图中的主要变量:N原始数据点数或剔除可疑数据后的合格数据点数N1可疑数据点数X一维数组,用于存放原始数据及合格数据中的x值Y一维数组,用于存放原始数据或合格数据中的y值82整理课件X1

一维数组,用于存放可疑数据点的x值Y1一维数组,用于存放可疑数据点的y值A回归直线截距B回归直线斜率R简单相关系数83整理课件SD标准偏差σER平均相对误差DALTA绝对值最大的残差ID残差绝对值最大的数据点序号U肖维奈特准则的k值84整理课件子程序LINEAR1A为一元线性回归计算子程序,比例6-1

中的子程序LINEAR1增加了标准偏差和平均相对误差的计

算;子程序RULES为肖维奈特准则中k值的计算程序。采用类似的方法,可以编写能剔除可疑数据的多元线性回归

计算程序框图。6.3.3计算实例85整理课件例6-6在定温下某反应的活化能(E)与压力

(P)呈直线关系:E=a+bP86整理课件由实验测得如下数据:试用肖维奈特准则剔除其中的可疑数据,并确定a和b的值,计算相关系

数、标准偏差和平均相对偏差。解按图6-7的程序框图编写源程序。源程序:’********************************************’Example6-6-------Eg6-6.frm’********************************************DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(500),Y(500),X1(500),Y1(500)DimXYAAsVariantClsN=10XYA=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=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+287整理课件NextI’----------IC=1:N1=0DoWhileIC=1CallLINEAR1A(N,X(),Y(),A,B,R,SD,ER)CallRULES(N,U)DALTA=0:ID=0ForI=1ToNDT=Abs(Y(I)-A-B*X(I))IfDT>DALTAThenDALTA=DTID=IEndIfNextIIfDALTA>U*SDThenINN=0ForI=1ToNIfI=IDThenN1=N1+1X1(N1)=X(I):Y1(N1)=Y(I)ElseINN=INN+1X(INN)=X(I):Y(INN)=Y(I)EndIfNextIN=N-1:IC=1ElseIC=0EndIfLoop’---------Print"N=";N:Print"N1=";N1Print"I","X1","Y1"ForI=1ToN1PrintI,X1(I),Y1(I)NextIPrint"A=";A:Print"B=";B:Print"R=";RPrint"SD=";SD:Print"ER=";ER;"%"EndSub’******************SubLINEAR1A(N,X(),Y(),A,B,R,SD,ER)’******************XT=0:YT=0:XX=0:YY=0:XY=0ForI=1ToNXT=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)NextIXXL=XX-XT*XT/N88整理课件YYL=YY-YT*YT/NXYL=XY-XT*YT/NB=XYL/XXLA=(YT-B*XT)/NR=XYL/Sqr(XXL*YYL)SD=0:ER=0ForI=1ToNSD=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)))NextISD=Sqr(SD/(N-2)):ER=(ER/N)*100EndSubSubRULES(N,U)SelectCaseN89整理课件执行结果:N=8N1=2IX1Y112802870A=39.80313111545988B=0.4172211350293531R=0.9910796148007959SD=0.1830558892195286ER=0.3321330871296535%从源程序的执行结果看,按肖维奈特准则,10个数据点中有两个属于可疑

数据,可疑数据的活化能测定值分别为80和70。90整理课件6.4

多项式拟合

在化学化工的实验或科研中,经常需要从一

组测定数据,例如从n对(xi,yi)数据,去求自变量x和

因变量y的近似函数关系式y=p(x)。从图形上

看,这是由给定的n个点(xi,yi)(i=1,2,…,n)作曲线91整理课件拟合。在曲线拟合中,多项式拟合问题占特殊的地位。

任何函数在一个比较小的范围内,可以用多项式

任意逼近。因此,在比较复杂的实际问题中,可以

不问y与各因素的确切关系,而用多项式拟合进行

分析和计算。92整理课件下面以多项式拟合为例,说明曲线拟合的方法和

计算程序。6.4.1方法概述设用下列m次多项式:

93整理课件拟合一组数据(xi,yi)(i=1,2,…,n),即曲线y=f(x)上已给定n个

点,用多项式求作该曲线的近似图形。这一问题与前述的插

值问题有类似之处。但插值问题要求近似曲线y=p(x)严格

地通过所给的n个点,这一要求将会使近似曲线y=p(x)保留

数据的全部测试点的测量误差。如果个别数据的误差很大,

那么插值的效果显然是不够理想的。鉴于这种情况,考虑放

弃严格通过所有结点(xi,yi)这一要求,而采用别的方法去构造

近似曲线,以尽可能反映所给数据的总趋势。曲线拟合的常94整理课件用方法仍然是最小二乘法,即残差平方和最小法。若以δi代表结点处的残差,则

残差的平方和为

95整理课件由于xi与yi为已知值,故Q是aj(j=0,1,2,…,m)的函数。由

多元函数的极值理论可知,要使Q最小,则系数ak必须满足下

式:

即

(

aj

)=

yi

96整理课件

上式变换后得:

令Sl=

97整理课件

将上两式代入式(6-24)可得m+1元线性方程组:

98整理课件式(6-25)称为正规方程组。这是m+1元线性方程组,其系数矩

阵是对称矩阵,可以证明,上述正规方程组有唯一解。所得的

m次多项式Pm(x)确能使残差平方和Q最小,故Pm(x)即为所求的

拟合多项式。它是函数y=f(x)的近似表达式,也可用它代替f

(x)作微分、积分计算。此外,必须注意的是上式计算中

≡1

。6.4.2程序框图99整理课件图6-8是多项式拟合的通用计算程序框图。图6-8多项