化工數(shù)值計(jì)算-6

1、回歸分析和曲線擬合生產(chǎn)過程和科學(xué)實(shí)驗(yàn)中,常用的變量大體可分兩類。一類為確定性變量,另一類為隨機(jī)變量。確定性變量是指兩個(gè)或多個(gè)變量之間有確定的關(guān)系.即其中某個(gè)變量的每個(gè)值,都與一變量的一個(gè)或幾個(gè)完全確定的值相對應(yīng),即它們之間存在著函數(shù)關(guān)系:例如,理想氣體的壓力P與摩爾體積V間,存在著確定的函數(shù)關(guān)系:但在實(shí)際問題中,由于變量之間的關(guān)系比較復(fù)雜,或由于生產(chǎn)或?qū)嶒?yàn)過程中不可避免地存在著誤差,使變量之間的關(guān)系具有不確定性,也就是說,某個(gè)變量對應(yīng)的,不是一個(gè)或幾個(gè)確定的值,而是整個(gè)集合的值,這時(shí),變量x和y間的關(guān)系,就稱為相關(guān)關(guān)系。例如,流體在圓形直管中做湍流時(shí)的情形,通過量綱分析可知,努塞爾特?cái)?shù)Nu、普

2、蘭特?cái)?shù)Pr和雷諾數(shù)Re之間存在著如下相關(guān)關(guān)系:這種關(guān)系的不確定性,表現(xiàn)為式中a和b的數(shù)值,在每次測量中不盡相同。不確定的原因,首先是影響該過程的因素甚多,有些因素至今尚未弄清;其次是受到實(shí)驗(yàn)過程中的偶然因素影響。這種不確定性關(guān)系并不說明上述三個(gè)量綱為1的數(shù)群之間無規(guī)律可循。相反,通過大量試驗(yàn),人們發(fā)現(xiàn),a和b的數(shù)值總是圍繞著某一定值波動,而且隨著試驗(yàn)次數(shù)的增多,a、b的數(shù)值趨于穩(wěn)定。a、b的穩(wěn)定值,可作為a和b的最佳估計(jì)值。在一定條件下,a=0.023,b=0.8。由此可見,通過大量試驗(yàn),是可以找到隱藏在隨機(jī)性后面的統(tǒng)計(jì)規(guī)律性的?；貧w分析和曲線擬合是一種處理變量相關(guān)關(guān)系的數(shù)理統(tǒng)計(jì)方法。用它可以

3、尋找隱藏在隨機(jī)性后面的統(tǒng)計(jì)規(guī)律性。函數(shù)與相關(guān)是兩種不同類型的變量關(guān)系,它們之間并無嚴(yán)格界限。一方面,相關(guān)的變量之間,并無確定的關(guān)系,但在一定的條件下,從一定的統(tǒng)計(jì)意義上看,它們之間又可能存在著某種確定的函數(shù)關(guān)系。另一方面,由于實(shí)際測定的數(shù)據(jù)中,總存在著誤差,即使是確定性變量,也會出現(xiàn)某些非確定性結(jié)果。6.1一元線性回歸一元線性回歸處理的是兩個(gè)變量之間的線性關(guān)系。所用的數(shù)學(xué)模型為一元線性代數(shù)模型,其模型方程式是對這種模型參數(shù)的估計(jì),就是根據(jù)原始數(shù)據(jù)點(diǎn)(x1, y1)、(x2, y2)、(xi, yi)、(xn, yn),確定式(6-1)中a、b的估計(jì)值。在實(shí)際體系中,自變量x與因變量y之間服從線

4、性關(guān)系的情況雖然不多,但在不少情況下,x、y之間存在著某種函數(shù)組合關(guān)系。例如f1 (x, y),f2 (x, y),設(shè)兩個(gè)函數(shù)之間服從線性關(guān)系f1與f2是不含待定系數(shù)的已知函數(shù)。若把f1 (x,y)與f2 (x,y)分別視為自變量與因變量,則仍可以借用線性模型去估計(jì)其參數(shù)值。這種方法稱為化直法。它在化學(xué)化工的實(shí)際問題中是常見的。例如單分子基元反應(yīng)AB的動力學(xué)方程式為對上式積分得式中,cA-t是不呈線性關(guān)系的函數(shù)。若對方程兩邊取對數(shù),上式可化為lncA-t的線性函數(shù):又例如,按照阿侖尼烏斯定律,反應(yīng)速率常數(shù)k與溫度T之間不呈線性關(guān)系:但lnk與1/T則呈線性關(guān)系:這些都是屬于可化為線性關(guān)系的例子

5、。一元線性代數(shù)模型中的待定參數(shù)a和b,稱為“估計(jì)值”。之所以稱為“估計(jì)”值,是因?yàn)閍,b的值是從實(shí)驗(yàn)值中通過數(shù)理統(tǒng)計(jì)方法確定的。圖6-1一元線性回歸6.1.1方法概述設(shè)有一組實(shí)驗(yàn)數(shù)據(jù)(x1,y1)、(x2,y2)、(xn,yn),自變量x與因變量y存在著式(6-1)的關(guān)系。當(dāng)x取值為xi時(shí),y的測定值為yi,計(jì)算值為yi*,并有由于參數(shù)a,b為未知值,故yi*也是未知值。若將全部實(shí)驗(yàn)數(shù)據(jù)標(biāo)繪在x-y圖中(見圖6-1),由于各種因素的影響,它們不會全部落在一條直線上,即n個(gè)yi不會與n個(gè)yi*完全重合,它們將隨機(jī)地分布在與xi呈線性關(guān)系的yi*的周圍。以i表示它們之間的差值,則有這里i就是誤差。

6、它反映了xi使yi偏離直線的各種影響因素的總和?，F(xiàn)在,要尋找一條最靠近各個(gè)數(shù)據(jù)點(diǎn)的直線,這條直線稱為回歸直線。由于回歸直線是一切直線中最接近各數(shù)據(jù)點(diǎn)(xi,yi)的,用它代表x與y之間的線性關(guān)系,比任何其他直線更為可靠。究竟如何確定回歸曲線中的參數(shù)a和b呢?目前最常用的方法就是最小二乘法,即殘差平方和最小法。式(6-3)中的誤差i又稱為殘差,表示第i個(gè)數(shù)據(jù)與回歸直線的偏離程度,則殘差平方和Q表示全部數(shù)據(jù)與回歸直線的總偏離程度。顯然Q是a和b的函數(shù):不用殘差和i的原因是i有正有負(fù),相加時(shí)可能彼此抵消,從而不能反映總的偏離程度,而用殘差的平方和不會發(fā)生這種現(xiàn)象。由多元函數(shù)的極值理論可知,要使Q值最

7、小,a、b必須滿足下列條件:即得式(6-6)稱為一元線性回歸的正規(guī)方程組,通過求解該方程組,可得:式(6-7)中等號右側(cè)的量全部取自原始數(shù)據(jù)。因此,就可以確定回歸系數(shù)a和b,完成參數(shù)估計(jì)。為了簡化a和b的表達(dá)式,定義:式中,、分別為xi和yi的平均值。xi與之差(xi-),稱為xi的離差;全部xi的離差平方和,稱為x的離差平方和,記為Lxx:yi與之差(yi-),稱為yi的離差;全部yi的離差平方和,稱為y的離差平方和,記為Lyy,同理再令Lxy為全部xi的離差與yi的離差乘積的總和:將以上關(guān)系式代入式(6-7),得由式(6-12)第二式可以看出,回歸直線是通過點(diǎn)(,)的。從力學(xué)觀點(diǎn)看,(,)

8、相當(dāng)于n個(gè)實(shí)驗(yàn)點(diǎn)(xi,yi)的重心,回歸直線是通過重心的。應(yīng)當(dāng)指出: 殘差i只用yi-y*i表示時(shí),表明yi有測量誤差,而xi無測量誤差;或表示與yi相比,xi的誤差很小。因此,測量誤差使實(shí)驗(yàn)點(diǎn)偏離回歸直線,都表現(xiàn)為yi偏離y*i。如果xi的誤差與yi的誤差相比,不可忽略,則兩者都必須考慮。這種情況比較復(fù)雜,此處不予介紹。 求回歸方程的計(jì)算過程中,不需要事先假定兩個(gè)變量之間必須有相關(guān)關(guān)系。即使是一組雜亂無章的數(shù)據(jù),也可以用最小二乘法繪制一條直線,以表示x與y的關(guān)系。顯然,這種情況下,繪制的直線并無實(shí)際意義。為了判斷兩個(gè)變量間線性關(guān)系的優(yōu)劣程度,引入一個(gè)新的指標(biāo)R,稱為簡單相關(guān)系數(shù),它的定義為

9、R值不同時(shí),數(shù)據(jù)點(diǎn)的分布情況如下。(1) R = 0圖6-2R = 0的數(shù)據(jù)點(diǎn)分布此時(shí)Lxy = 0,b = 0。即回歸直線平行于x軸,y的變化與x無關(guān),表示數(shù)據(jù)點(diǎn)的分布是無規(guī)則的,如圖6-2所示。但亦有當(dāng)R = 0時(shí),x與y確實(shí)存在明顯相關(guān)性的情況。這種情形,不能應(yīng)用線性回歸方法,只能用化直線法或曲線擬合法處理。(2) 0 |R| 0時(shí),b 0。數(shù)據(jù)點(diǎn)的y值隨著x增加而增加,這種情況稱為x與y正相關(guān)。R 0時(shí),b 0。數(shù)據(jù)點(diǎn)的y值隨著x增加而減小,這種情況稱為x與y負(fù)相關(guān)。R的絕對值越小,數(shù)據(jù)點(diǎn)沿回歸直線越分散。圖6-301的數(shù)據(jù)點(diǎn)分布1的數(shù)據(jù)點(diǎn)分布(3) |R| = 1x與y完全相關(guān)。全部

10、數(shù)據(jù)點(diǎn)均落在回歸直線上。若x與y為非線性相關(guān),但經(jīng)變量變換后,用回歸直線的方法處理,所求得的回歸系數(shù)僅對變換后的變量是最佳的,而對原變量來說則并非最佳,但通常還能令人滿意,此時(shí)應(yīng)注意原變量的殘差平方和并非最小。由以上討論可知,相關(guān)系數(shù)R的絕對值在0與1之間,而且越接近于1,其線性關(guān)系越密切,那么|R|與1接近到什么程度,才能說明x與y之間存在線性相關(guān)關(guān)系呢?要回答這個(gè)問題,就要對相關(guān)系數(shù)進(jìn)行顯著性檢驗(yàn)。由于篇幅所限,有關(guān)相關(guān)系數(shù)的顯著性檢驗(yàn)和回歸方程的方差分析等問題將不在此討論。如有需要,可參考有關(guān)數(shù)理統(tǒng)計(jì)方面的書籍。6.1.2程序框圖圖6-4是一元線性回歸的通用計(jì)算程序框圖。程序框圖中的主要

11、變量:N數(shù)據(jù)點(diǎn)數(shù)X、Y 一維數(shù)組,用于存放原始數(shù)據(jù)中的x和y值XXLx離差平方和LxxYYLy離差平方和LyyXYLx離差與y離差乘積總和LxyA回歸直線截距aB回歸直線斜率bR簡單相關(guān)系數(shù)6.1.3計(jì)算實(shí)例例6-1已知某反應(yīng)的速率常數(shù)k與熱力學(xué)溫度T的實(shí)驗(yàn)數(shù)據(jù)如下:試求k-T的關(guān)系式。解通常反應(yīng)速率常數(shù)與熱力學(xué)溫度的關(guān)系,服從阿侖尼烏斯定律:k=Ae-E/RT式中,E為反應(yīng)活化能;T為熱力學(xué)溫度;R為氣體通用常數(shù)。上式取對數(shù),且令y = lnk,x = -1/T,可得y=lnA+x。按圖6-4,用一元線性回歸求得A = 1.966109/min,E = 79.571kJ/mol。將實(shí)驗(yàn)數(shù)據(jù)點(diǎn)

12、和利用關(guān)系式獲得的計(jì)算點(diǎn)一起繪制在圖6-5中。 源程序:*Example 6-1-Eg6-1.frm*DefDbl A-H, O-ZPrivate Sub Command1_Click()Dim X(50), Y(50)Dim XYA As VariantClsN = 5XYA = Array(363, 0.00718, 373, 0.01376, 383, 0.02701, _393, 0.05221, 403, 0.09718)K = 0For I = 1 To NX(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2X(I) = -1 / X(I)Y(I)

13、 = Log(Y(I)Next ICallLINEAR1(N, X(), Y(), A, B, R)A = Exp(A): E = B* 8.314Print A = ; A: Print E = ; EPrint R = ; REnd Sub*Sub LINEAR1(N, X(), Y(), A, B, R)*XT = 0: YT = 0: XX = 0: YY = 0: XY = 0For I = 1 To NXT = XT + X(I): YT = YT + Y(I)XX = XX + X(I) * X(I): YY = YY + Y(I) * Y(I)XY = XY + X(I) *

14、Y(I)Next IXXL = XX - XT * XT / NYYL = YY - YT * YT / NXYL = XY - XT * YT / NB = XYL / XXLA = (YT - B * XT) / NR = XYL / Sqr(XXL * YYL)End Sub執(zhí)行結(jié)果:A =1966349283.054212 (指前因子)E =79570.97618674007 (活化能)R =.999718315533107(相關(guān)系數(shù))源程序中將一元線性回歸計(jì)算安排在子程序LINEAR1中。例6-2某水樣BOD測定數(shù)據(jù)如下:試確定該水樣中有機(jī)物生物氧化降解反應(yīng)的經(jīng)驗(yàn)速率方程表達(dá)式。解通

15、常水體中有機(jī)物生物氧化降解反應(yīng)的經(jīng)驗(yàn)速率方程服從下列方程:式中,BOD、L0為分別為t時(shí)和初始時(shí)刻的生化需氧量(mg/L);k為BOD的降解系數(shù),即耗氧系數(shù)/d-1托馬斯(Thomas)提出將1-e-kt按冪級數(shù)展開如下:與此展開相似的表達(dá)式有:兩展開式僅第四項(xiàng)出現(xiàn)微小差別,故可以近似地取即 BOD=L0整理得=+t取y =、x = t按一元線性回歸計(jì)算a = 0.25778、b = 0.01371,從而解得k = 6b/a = 0.31907d-1,L0 = 182.963mg/L。源程序(僅列出主程序,回歸子程序LINEAR1同例6-1):*Example 6-2-Eg6-2.frm*De

16、fDbl A-H, O-ZPrivate Sub Command1_Click()Dim X(50), Y(50)Dim XYA As VariantClsN = 10XYA = Array(1, 58, 2, 85, 3, 107, 4, 125, 5, 138, _ 6, 147, 7, 155, 8, 161, 9, 167, 10, 170)K = 0For I = 1 To NX(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2Y(I) = (X(I) / Y(I) (1 / 3)Next ICall LINEAR1(N, X(), Y(), A,

17、B, R)Print A = ; A: Print B = ; BBK = 6* B / ABL0 = 1 / BK / A / A / APrint L0 = ; BL0: Print K = ; BKPrint R = ; REnd Sub執(zhí)行結(jié)果:A=.2577799110470602 B=1.370838615288584D-02L0=182.9634248033556 K=.3190718647672257R=.9900105997750359此外,塞里奧特(Theriaut)提出了BOD公式的另一種解法:式中,k為待估的k的近似值;h為k的允許偏差量。從而有因h甚小,e-ht1-h

18、t,故上式變?yōu)?式中,a = L0,d = L0h,x1=1-e-kt,x2=te-kt。這樣可以首先假設(shè)k 的初始值為k0,利用實(shí)驗(yàn)數(shù)據(jù)通過二元線性回歸(見下節(jié)例6-5)確定出a和d,并求出L0 = a、h = d/L0;若|h| (誤差允許值),則令k = k0 + h,重新進(jìn)行線性回歸計(jì)算,直至|h| m + 1才能求出上式中的m + 1個(gè)回歸系數(shù)。同樣由多元函數(shù)的極值理論可知,要使Q值最小,a0和aj必須滿足下列條件:式(6-15)經(jīng)整理可得:式(6-16)稱為多元線性回歸模型的正規(guī)方程組。它是一個(gè)m+1元的線性代數(shù)方程組。由于xij和yi已知,故可求得m+1個(gè)待定系數(shù)a0,a1,am

19、。實(shí)際計(jì)算時(shí),一般作如下處理:先將式(6-16)的第一式寫成然后將式(6-17)代入方程組(6-16)的第2至第m+1式,重新組成一個(gè)m元線性方程組,其中有a1,a2,am等m個(gè)待定系數(shù)。通過求解此m元線性方程組,獲得系數(shù)a1,a2,am,再代回式(6-17),求得a0。為簡化計(jì)算,用表示第j個(gè)x的平均值,表示y的平均值,則用Ljk表示第j個(gè)x離差與第k個(gè)x離差乘積之和,則用Lyy表示y離差的平方和,則用Ljy表示第j個(gè)x離差與y離差乘積之和,則將式(6-17)分別代入式(6-16)的第2至m+1式,經(jīng)簡化整理可得如下m元線性方程組:可用主元素消去法求解此式,然后將求得的a1,a2,am代入式

20、(6-17),求出a0,從而完成對多元線性回歸模型的參數(shù)估計(jì)。多元線性回歸的計(jì)算中,常用復(fù)相關(guān)系數(shù)衡量數(shù)據(jù)點(diǎn)之間的線性優(yōu)劣。復(fù)相關(guān)系數(shù)定義如下:式中,U稱為回歸平方和:應(yīng)當(dāng)指出,并非所有曲線都可以按這種方法處理。例如拋物線就不能通過變量變換把它化為直線。但是如果令x1 = x,x2 = x2,則上式就化成一個(gè)包含兩個(gè)自變量的線性方程從而將拋物線按二元線性回歸計(jì)算。對于含多變量的任意多項(xiàng)式也可以通過類似的變換,把它們轉(zhuǎn)化成多元線性回歸計(jì)算。6.2.2程序框圖圖6-6是多元線性回歸的通用計(jì)算程序框圖。圖6-6(a) 多元線性回歸的通用計(jì)算程序框圖(1)圖6-6(b) 多元線性回歸的通用計(jì)算程序框圖

21、(2)程序框圖中的主要變量:N 數(shù)據(jù)點(diǎn)數(shù)M多元線性模型元數(shù)X二維數(shù)組,用于存放原始數(shù)據(jù)的x值Y一維數(shù)組,用于存放原始數(shù)據(jù)的y值YP值YYLLyy值XP一維數(shù)組,用于存放值A(chǔ)二維數(shù)組,用于存放m元線性方程組的系數(shù)LjkB一維數(shù)組,用于存放m元線性方程組的常數(shù)項(xiàng)LjyC一維數(shù)組,用于存放多元線性模型的系數(shù)aj(j = 0,1,M)R復(fù)相關(guān)系數(shù)R0U回歸平方和Q殘差平方和子程序XYF為列主元消去法求解線性方程組的程序,可參見圖5-2和圖5-3。6.2.3計(jì)算實(shí)例例6-4已知某溶液由兩種物質(zhì)組成,cA為物質(zhì)A的濃度(g/L),cB為物質(zhì)B的濃度(g/L), 為溶液的黏度(mPas)。設(shè)數(shù)學(xué)模型為=a0

22、+a1cA+a2cB試根據(jù)下列實(shí)驗(yàn)數(shù)據(jù),確定a0、a1、a2的值。解按圖6-6編寫計(jì)算源程序。源程序:* 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

23、.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 NFor J = 1 To M:X(I, J) = XYA(K):K = K + 1: Next JY(I) = XYA(K): K = K + 1 Next I Call LINEAR2(N, M, X(), Y(), C(), R) Print Tab(4); * Results

24、 * For J = 0 To MPrint 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 AverageFor I = 1 To N: YP = YP + Y(I): Next IYP = YP / NYYL = 0 = = Lyy AverageFor I = 1 To N: YYL = YYL + (Y(I) - YP) * (Y

25、(I) - YP): Next IFor J = 1 To M = = Xj AverageXP(J) = 0For I = 1 To N: XP(J) = XP(J) + X(I, J): Next IXP(J) = XP(J) / NNext JFor J = 1 To M = = LjyXYL = 0For I = 1 To NXYL = XYL + (X(I, J) - XP(J) * (Y(I) - YP)Next IB(J) = XYLNext JFor J = 1 To M = = LjkFor K = 1 To MXXL = 0For I = 1 To NXXL = XXL +

26、 (X(I, J) - XP(J) * (X(I, K) - XP(K)Next IA(J, K) = XXLNext KNext JCall XYF(A(), B(), M, C()C(0) = YP =a0For J = 1 To MC(0) = C(0) - C(J) * XP(J)Next JU = 0: Q = 0 = RFor I = 1 To NYI = C(0)For J = 1 To MYI = YI + C(J) * X(I, J)Next JU = U + (YI - YP) * (YI - YP)Q = Q + (YI - Y(I) * (YI -Y(I)Next IR

27、 = Sqr(U / (U + Q) -or R = SQR(U/YYL) End Sub執(zhí)行結(jié)果:* Results *A(0)=-27.4324958A(1)=0.2327103A(2)=0.4095299R= 0.7472224源程序中將多元線性回歸安排在子程序LINEAR2中,子程序XYF和XLZY見例5-3。例6-5按塞里奧特方法確定例6-2的BOD公式中的參數(shù)。解由前敘述可知,塞里奧特方法計(jì)算機(jī)求解時(shí)要應(yīng)用二元線性回歸,且無常數(shù)項(xiàng),即a00。對a00時(shí)的多元線性回歸,求解回歸方程系數(shù)的m元線性方程組同樣具有方程組(6-18)的形式,不過系數(shù)矩陣和常數(shù)項(xiàng)必須用下列式子進(jìn)行計(jì)算:將例6

28、-4的源程序加以修改即可得到本例的計(jì)算源程序。源程序(子程序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

29、 NT(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2 Next IE = 0.000001: X0 = 0.5: H = 100 * E Do While Abs(H) EFor I = 1 To NX(I, 1) = 1 - Exp(-X0 * T(I)X(I, 2) = T(I) * Exp(-X0 * T(I)Next ICall 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

30、= ; 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 NXYL = XYL + X(I, J) * Y(I) Next I B(J) = XYL Next J For J = 1 To M = Sjk

31、 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

32、+ Q) End Sub執(zhí)行結(jié)果:* Results *L0= 173.428617756483(mg/L)K=0.3306173539389562 (d-1)R= 0.9956168注意按這種方法進(jìn)行計(jì)算時(shí),k的初始值必須在真實(shí)值附近選取,否則將得出錯(cuò)誤的結(jié)果。注意本例源程序中的多元線性回歸子程序LINEAR20,只適用于無常數(shù)項(xiàng)的多元線性方程的回歸。6.3剔除可疑數(shù)據(jù)及其計(jì)算程序6.3.1剔除可疑數(shù)據(jù)的方法在線性回歸計(jì)算中,假定每個(gè)測定數(shù)據(jù)與回歸結(jié)果之間的誤差均在隨機(jī)誤差允許的范圍之內(nèi)。然而,由于測量誤差或過失誤差等多種原因,在一組實(shí)驗(yàn)值中,誤差往往會超出隨機(jī)誤差的允許范圍。這些數(shù)據(jù),稱為

33、可疑數(shù)據(jù)。為保證回歸結(jié)果的可靠性,必須剔除這些可疑的數(shù)據(jù)。剔除可疑數(shù)據(jù),應(yīng)當(dāng)有一個(gè)科學(xué)的標(biāo)準(zhǔn)。這個(gè)標(biāo)準(zhǔn)就是統(tǒng)計(jì)判據(jù),屬于統(tǒng)計(jì)判據(jù)的剔除準(zhǔn)則有多種。以一元線性回歸為例,其代數(shù)模型為y = a + bx。若自變量x無測量誤差,則y的標(biāo)準(zhǔn)偏差為式中,n為原始數(shù)據(jù)點(diǎn)數(shù);m為回歸模型中自變量的個(gè)數(shù),對一元線性回歸m = 1;i為殘差,即 i=yi-a、b是按最小二乘法求出的最佳估計(jì)值。根據(jù)數(shù)理統(tǒng)計(jì)分析,合理的數(shù)據(jù),其殘差不應(yīng)超出 的k倍。若取k = 3,便是常用的3 準(zhǔn)則。據(jù)此,可以把殘差絕對值超過3 的個(gè)別數(shù)據(jù)(xi,yi),判為可疑數(shù)據(jù)而加以剔除。必須指出,3 準(zhǔn)則是以數(shù)據(jù)點(diǎn)數(shù)n為前提的,當(dāng)n為有限

34、值時(shí),3 判據(jù)并不十分可靠。下面介紹一種廣泛采用的判據(jù),即所謂肖維奈特準(zhǔn)則。按肖維奈特準(zhǔn)則,若n次等精度測量中,有某個(gè)測量值yi,其殘差的絕對值超出k ,就可以認(rèn)為是可疑數(shù)據(jù)而予以剔除。表6-1列出了肖維奈特準(zhǔn)則中與n相對應(yīng)的k值。表6-1肖維奈特準(zhǔn)則的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.2635

35、2.455003.29122.04222.28402.5010003.48132.07232.30502.5820003.66142.10242.32602.6450003.89使用這個(gè)準(zhǔn)則時(shí),可根據(jù)回歸結(jié)果,對全部實(shí)驗(yàn)值進(jìn)行逐級檢查,把屬于可疑數(shù)據(jù)的實(shí)驗(yàn)值選出。若發(fā)現(xiàn)不止一個(gè)可疑數(shù)據(jù),則應(yīng)把其中殘差絕對值最大者剔除,然后重新計(jì)算 值。根據(jù)新的 值,再次用肖維奈特準(zhǔn)則進(jìn)行檢查。每次只剔除一個(gè)可疑數(shù)據(jù),其余數(shù)據(jù)重新進(jìn)行回歸,直至回歸所用的數(shù)據(jù)中不再含有可疑數(shù)據(jù)為止。6.3.2剔除可疑數(shù)據(jù)的計(jì)算程序框圖圖6-7是具有剔除可疑數(shù)據(jù)功能的一元線性回歸通用計(jì)算程序框圖,整個(gè)計(jì)算過程分為輸入原始數(shù)據(jù)、一元

36、線性回歸計(jì)算、確定肖維奈特準(zhǔn)則的k值、確定殘差絕對值最大的數(shù)據(jù)點(diǎn)、剔除最可疑數(shù)據(jù)點(diǎn)(即殘差絕對值最大的數(shù)據(jù)點(diǎn))。圖6-7具有剔除可疑數(shù)據(jù)功能的一元線性回歸通用計(jì)算程序框圖程序框圖中的主要變量:N原始數(shù)據(jù)點(diǎn)數(shù)或剔除可疑數(shù)據(jù)后的合格數(shù)據(jù)點(diǎn)數(shù)N1可疑數(shù)據(jù)點(diǎn)數(shù)X一維數(shù)組,用于存放原始數(shù)據(jù)及合格數(shù)據(jù)中的x值Y一維數(shù)組,用于存放原始數(shù)據(jù)或合格數(shù)據(jù)中的y值X1一維數(shù)組,用于存放可疑數(shù)據(jù)點(diǎn)的x值Y1一維數(shù)組,用于存放可疑數(shù)據(jù)點(diǎn)的y值A(chǔ)回歸直線截距B回歸直線斜率R簡單相關(guān)系數(shù)SD標(biāo)準(zhǔn)偏差ER平均相對誤差DALTA絕對值最大的殘差I(lǐng)D殘差絕對值最大的數(shù)據(jù)點(diǎn)序號U肖維奈特準(zhǔn)則的k值子程序LINEAR1A為一元線性回

37、歸計(jì)算子程序,比例6-1中的子程序LINEAR1增加了標(biāo)準(zhǔn)偏差和平均相對誤差的計(jì)算;子程序RULES為肖維奈特準(zhǔn)則中k值的計(jì)算程序。采用類似的方法,可以編寫能剔除可疑數(shù)據(jù)的多元線性回歸計(jì)算程序框圖。6.3.3計(jì)算實(shí)例例6-6在定溫下某反應(yīng)的活化能(E)與壓力(P)呈直線關(guān)系:E = a + bP由實(shí)驗(yàn)測得如下數(shù)據(jù):試用肖維奈特準(zhǔn)則剔除其中的可疑數(shù)據(jù),并確定a和b的值,計(jì)算相關(guān)系數(shù)、標(biāo)準(zhǔn)偏差和平均相對偏差。解按圖6-7的程序框圖編寫源程序。源程序:*Example 6-6-Eg6-6.frm*DefDbl A-H, O-ZPrivate Sub Command1_Click()Dim X(500

38、), Y(500), X1(500), Y1(500)Dim XYA As VariantClsN = 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 = 0For I = 1 To NX(I) = XYA(K): Y(I) = XYA(K + 1): K = K + 2Next I-IC = 1: N1 = 0Do While IC = 1 Call LINEAR1A(N, X(), Y(), A, B, R, SD, ER) Cal

39、l RULES(N, U) DALTA = 0: ID = 0 For I = 1 To N DT = Abs(Y(I) - A - B * X(I) If DT DALTA ThenDALTA = DTID = I End If Next I If DALTA U * SD ThenINN = 0For I = 1 To NIf 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 IfNext IN = N - 1: IC = 1 Else

40、IC = 0 End IfLoop-Print N = ; N: Print N1 = ; N1Print I, X1, Y1For I = 1 To N1Print I, X1(I), Y1(I)Next IPrint A = ; A: Print B = ; B: Print R = ; RPrint 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 = 0For I = 1 To NXT = XT + X(I)

41、: YT = YT + Y(I)XX = XX + X(I) * X(I): YY = YY + Y(I) * Y(I)XY = XY + X(I) * Y(I)Next IXXL = XX - XT * XT / N YYL = YY - YT * YT / NXYL = XY - XT * YT / NB = XYL / XXLA = (YT - B * XT) / NR = XYL / Sqr(XXL * YYL)SD = 0: ER = 0For I = 1 To NSD = SD + (Y(I) - A - B * X(I) * (Y(I) - A - B * X(I)ER = ER

42、 + Abs(Y(I) - A - B * X(I) / (A + B * X(I)Next ISD = Sqr(SD / (N - 2): ER = (ER / N) * 100End SubSub RULES(N, U) Select Case N執(zhí)行結(jié)果: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%從源程序的執(zhí)行結(jié)果看,按肖維奈特準(zhǔn)則,10個(gè)數(shù)據(jù)點(diǎn)中有

43、兩個(gè)屬于可疑數(shù)據(jù),可疑數(shù)據(jù)的活化能測定值分別為80和70。6.4多項(xiàng)式擬合在化學(xué)化工的實(shí)驗(yàn)或科研中,經(jīng)常需要從一組測定數(shù)據(jù),例如從n對(xi,yi)數(shù)據(jù),去求自變量x和因變量y的近似函數(shù)關(guān)系式y(tǒng) = p(x)。從圖形上看,這是由給定的n個(gè)點(diǎn)(xi,yi)(i = 1,2,n)作曲線擬合。在曲線擬合中,多項(xiàng)式擬合問題占特殊的地位。任何函數(shù)在一個(gè)比較小的范圍內(nèi),可以用多項(xiàng)式任意逼近。因此,在比較復(fù)雜的實(shí)際問題中,可以不問y與各因素的確切關(guān)系,而用多項(xiàng)式擬合進(jìn)行分析和計(jì)算。下面以多項(xiàng)式擬合為例,說明曲線擬合的方法和計(jì)算程序。6.4.1方法概述設(shè)用下列m次多項(xiàng)式 :擬合一組數(shù)據(jù)(xi,yi)(i =

44、1,2,n),即曲線y = f (x)上已給定n個(gè)點(diǎn),用多項(xiàng)式求作該曲線的近似圖形。這一問題與前述的插值問題有類似之處。但插值問題要求近似曲線y = p (x)嚴(yán)格地通過所給的n個(gè)點(diǎn),這一要求將會使近似曲線y = p (x)保留數(shù)據(jù)的全部測試點(diǎn)的測量誤差。如果個(gè)別數(shù)據(jù)的誤差很大,那么插值的效果顯然是不夠理想的。鑒于這種情況,考慮放棄嚴(yán)格通過所有結(jié)點(diǎn)(xi,yi)這一要求,而采用別的方法去構(gòu)造近似曲線,以盡可能反映所給數(shù)據(jù)的總趨勢。曲線擬合的常用方法仍然是最小二乘法,即殘差平方和最小法。若以i代表結(jié)點(diǎn)處的殘差,則殘差的平方和為由于xi與yi為已知值,故Q是aj(j = 0,1,2,m)的函數(shù)。由

45、多元函數(shù)的極值理論可知,要使Q最小,則系數(shù)ak必須滿足下式:即(aj)=yi上式變換后得:令 Sl=將上兩式代入式(6-24)可得m+1元線性方程組:式(6-25)稱為正規(guī)方程組。這是m+1元線性方程組,其系數(shù)矩陣是對稱矩陣,可以證明,上述正規(guī)方程組有唯一解。所得的m次多項(xiàng)式Pm(x)確能使殘差平方和Q最小,故Pm(x)即為所求的擬合多項(xiàng)式。它是函數(shù)y = f (x)的近似表達(dá)式,也可用它代替f (x)作微分、積分計(jì)算。此外,必須注意的是上式計(jì)算中1。6.4.2程序框圖圖6-8是多項(xiàng)式擬合的通用計(jì)算程序框圖。圖6-8多項(xiàng)式擬合的通用計(jì)算程序框圖程序框圖中的主要變量:M多項(xiàng)式次數(shù)N原始數(shù)據(jù)組數(shù)X、Y一維數(shù)組,用于存放原始數(shù)據(jù)的x、y值S一維數(shù)組,用于存放正規(guī)方程組的Sl數(shù)值T一維數(shù)組,用于存放正規(guī)方程組各常數(shù)項(xiàng)Tk的值B二維數(shù)組,用于存放正規(guī)方程組的系數(shù)矩陣A一