多元模型回歸與分析ppt課件

1、 多元數(shù)據(jù)模型回歸與分析一、實(shí)驗(yàn)數(shù)據(jù)分析一、實(shí)驗(yàn)數(shù)據(jù)分析q由實(shí)驗(yàn)數(shù)據(jù)回歸模型，得到模型參數(shù)前，對(duì)數(shù)據(jù)自變量間的由實(shí)驗(yàn)數(shù)據(jù)回歸模型，得到模型參數(shù)前，對(duì)數(shù)據(jù)自變量間的線性相關(guān)性進(jìn)展檢驗(yàn)，是發(fā)現(xiàn)回歸模型運(yùn)用的可靠性和準(zhǔn)確線性相關(guān)性進(jìn)展檢驗(yàn)，是發(fā)現(xiàn)回歸模型運(yùn)用的可靠性和準(zhǔn)確性受限制的有效方法。性受限制的有效方法。q因自變量間的線性相關(guān)性，使得無(wú)法區(qū)分它們對(duì)因變量的作因自變量間的線性相關(guān)性，使得無(wú)法區(qū)分它們對(duì)因變量的作用用;q回歸模型參數(shù)時(shí)會(huì)遇到幾乎是奇特的數(shù)據(jù)矩陣，這樣的模型回歸模型參數(shù)時(shí)會(huì)遇到幾乎是奇特的數(shù)據(jù)矩陣，這樣的模型參數(shù)有很大的不確定性參數(shù)有很大的不確定性(95的參數(shù)置信度范圍寬的參數(shù)置信度

2、范圍寬)。 q例：回歸二氧化硫的催化氧化速率方程：例：回歸二氧化硫的催化氧化速率方程： 裝有載鉑氧化鋁催化劑顆粒的微分固定床反響器中，測(cè)定二裝有載鉑氧化鋁催化劑顆粒的微分固定床反響器中，測(cè)定二氧化硫的催化氧化速率?？倝簽檠趸虻拇呋趸俾??？倝簽?90 mmHg時(shí)，記錄流體相的時(shí)，記錄流體相的組成分壓，有下表所示的速率結(jié)果，經(jīng)過這些數(shù)據(jù)求取二氧組成分壓，有下表所示的速率結(jié)果，經(jīng)過這些數(shù)據(jù)求取二氧化硫的催化氧化速率方程?；虻拇呋趸俾史匠?。 二氧化硫的催化氧化速率二氧化硫的催化氧化速率 r分 壓 (atm)mol/g.hSO3SO2O20.020.04280.02550.1860.040.

3、03310.03530.1900.060.02720.04090.1930.080.02360.04430.1950.100.02140.04640.1960.120.02010.04760.197表表8 82 2 二氧化硫的催化氧化速率二氧化硫的催化氧化速率 兩種模型的非線性回歸兩種模型的非線性回歸q1、普通的指數(shù)速率方程方式、普通的指數(shù)速率方程方式 cObSOaSOPPkPr223(8.2.1) k=0.517113.3；a=-1.987.02；b=-0.2164.556；c=6.078124.7 擬合結(jié)果：擬合結(jié)果：參數(shù)的參數(shù)的9595置信度太寬，模型參數(shù)不可靠。置信度太寬，模型參數(shù)不可

4、靠。 q2、根據(jù)原子氧的吸附機(jī)理，得到的速率方程式、根據(jù)原子氧的吸附機(jī)理，得到的速率方程式Smith，Chemical Engineering Kinetics, 3rd Ed., 1981, McGraw-Hill, P.374212/1)(3322SOSOKOSOBPAPPPr(8.2.2) K=73K=73，為反響平衡常數(shù)，為反響平衡常數(shù)A=0.1017A=0.10170.09580.0958；B=16.02B=16.024.33 4.33 擬合結(jié)果：擬合結(jié)果：與方程與方程(8.2.1)(8.2.1)相比，方程參數(shù)的置信度有了顯著改善。相比，方程參數(shù)的置信度有了顯著改善。 對(duì)速率方程的進(jìn)

5、一步分析對(duì)速率方程的進(jìn)一步分析 假設(shè)把方程假設(shè)把方程(8.2.2)改寫為：改寫為： 2/122233)(OSOSOSOPPBPArKP(8.2.3) 將模型參數(shù)代入計(jì)算并以方程左邊為橫坐標(biāo)、右邊為縱坐標(biāo)將模型參數(shù)代入計(jì)算并以方程左邊為橫坐標(biāo)、右邊為縱坐標(biāo)作圖。作圖。0.0200.0300.040PSO30.1850.1900.1950.200PO2-0.10.00.10.2RHS of Eq. (8.2.3)Linear dependence of Pso3 and Po2RHS of Eq. (8.2.3)結(jié)果并不是斜結(jié)果并不是斜率為率為-1-1的直線。的直線。闡明表所給的闡明表所給的速率數(shù)

6、據(jù)沒有速率數(shù)據(jù)沒有足夠的信息來(lái)足夠的信息來(lái)闡明速率方程闡明速率方程中的逆反響奉中的逆反響奉獻(xiàn)。獻(xiàn)。如將如將SO3SO3分壓分壓對(duì)對(duì)O2O2分壓作分壓作圖，這兩分圖，這兩分壓間有近似壓間有近似線性關(guān)系。線性關(guān)系。所以方程所以方程8.2.18.2.1的的置信區(qū)間范置信區(qū)間范圍大。圍大。二、回歸模型的選擇二、回歸模型的選擇1 q例：水飽和蒸汽壓的模型回歸例：水飽和蒸汽壓的模型回歸 q水的蒸汽壓數(shù)據(jù)選用的溫度范圍為水的蒸汽壓數(shù)據(jù)選用的溫度范圍為0120 q三參數(shù)的三參數(shù)的Antoine方程：方程： q四參數(shù)的四參數(shù)的Riedel回歸方程：回歸方程： q五參數(shù)回歸方程參考五參數(shù)回歸方程參考Thek-St

7、iel的蒸汽壓預(yù)測(cè)方程提出：的蒸汽壓預(yù)測(cè)方程提出： )/(lnCTBAP6ln/lnDTTCTBAPTETDCTTBAPln/ln2(8.2.4) (8.2.5) (8.2.6) 水飽和蒸汽壓的模型回歸結(jié)果水飽和蒸汽壓的模型回歸結(jié)果參數(shù)Antoine方程改進(jìn)Thek-Stiel方程A18.5587.5132B-3973.2-10.449C-39.9832.8683D-.064796E-6.8475R20.99999981.0表表8-38-3 水飽和蒸汽壓的方程擬合結(jié)果水飽和蒸汽壓的方程擬合結(jié)果 擬合度非常接近擬合度非常接近1 1，闡明擬合是勝利的，但實(shí)踐上用，闡明擬合是勝利的，但實(shí)踐上用Ant

8、oineAntoine方程來(lái)擬合回歸得到的結(jié)果不理想，闡明僅從擬方程來(lái)擬合回歸得到的結(jié)果不理想，闡明僅從擬合度上來(lái)判別結(jié)果的好壞是不夠的。為什么呢？合度上來(lái)判別結(jié)果的好壞是不夠的。為什么呢？ Model: lnP=A+B/(T+C)y=(18.55832)+(-3973.1923)/(x+(-39.983344) T, K LnP12345678280300320340360380Water vapor fittingwith Antoine Eq.因變量與殘差關(guān)系圖因變量與殘差關(guān)系圖 )ln()ln(,calciobsiiPPq殘差定義：殘差定義： (8.2.7) q調(diào)查模型參數(shù)估計(jì)方法的兩

9、個(gè)根本假設(shè)：調(diào)查模型參數(shù)估計(jì)方法的兩個(gè)根本假設(shè)：q參數(shù)估計(jì)的誤差相互不相關(guān)聯(lián)，是隨機(jī)的。參數(shù)估計(jì)的誤差相互不相關(guān)聯(lián)，是隨機(jī)的。q估計(jì)誤差符合正態(tài)分布。估計(jì)誤差符合正態(tài)分布。檢查模型適宜體系數(shù)據(jù)程度的最有效方法之一檢查模型適宜體系數(shù)據(jù)程度的最有效方法之一是對(duì)因變量與殘差作圖，察看其分布情況。是對(duì)因變量與殘差作圖，察看其分布情況。 Antoine方程擬合的殘差方程擬合的殘差 殘差雖然很小，但其分殘差雖然很小，但其分布不是隨機(jī)的。布不是隨機(jī)的。殘差的分布同正態(tài)分布相殘差的分布同正態(tài)分布相比，有較大的差距。比，有較大的差距。 兩方面的結(jié)果充分闡明了擬合回歸的兩方面的結(jié)果充分闡明了擬合回歸的Antoin

10、eAntoine方程還不能充方程還不能充分反映蒸汽壓與溫度間的關(guān)系，呵斥殘差間存在關(guān)聯(lián)。分反映蒸汽壓與溫度間的關(guān)系，呵斥殘差間存在關(guān)聯(lián)。采用采用RiedelRiedel方程擬合得到的也是類似的結(jié)果。方程擬合得到的也是類似的結(jié)果。lnPResidual Values-0.003-0.002-0.0010.0000.0010.0020.00312345678Residual fromAntoine Eq fitExpectedNormalFrequency Distribution: ResidualsNo of obs02468101214161820222426-0.003-0.003-0.0

11、02-0.002-0.001-5.000e-40.0005.000e-40.0010.0020.0020.0030.0030.004改良改良Thek-Stiel方程方程的擬合結(jié)果方程方程的擬合結(jié)果擬合誤差比擬合誤差比AntoineAntoine方程方程小了近一個(gè)數(shù)量級(jí)，而且小了近一個(gè)數(shù)量級(jí)，而且殘差分布是隨機(jī)分布的。殘差分布是隨機(jī)分布的。誤差分布根本符合正態(tài)誤差分布根本符合正態(tài)分布。分布。改良改良Thek-StielThek-Stiel方程方程描畫水飽和蒸汽壓的適宜模型。方程方程描畫水飽和蒸汽壓的適宜模型。 lnPvPrResidual Values-7e-4-5e-4-3e-4-1e-41e

12、-43e-45e-4-11-10-9-8-7-6-5-4Eq.(8.2.6) fit Frequency Distribution: ResidualsNo of obs01020304050607080-7e-4-6e-4-5e-4-4e-4-3e-4-2e-4-1e-401e-42e-43e-44e-4Eq.(6) fitExpectedNormal二、回歸模型的選擇二、回歸模型的選擇2q前面闡明了模型參數(shù)較少時(shí)會(huì)出現(xiàn)擬合殘差前面闡明了模型參數(shù)較少時(shí)會(huì)出現(xiàn)擬合殘差的分布不是隨機(jī)的，而是呈現(xiàn)某種分布，相的分布不是隨機(jī)的，而是呈現(xiàn)某種分布，相互關(guān)聯(lián)?；リP(guān)聯(lián)。q在模型回歸擬合數(shù)據(jù)的過程中，如模型

13、參數(shù)在模型回歸擬合數(shù)據(jù)的過程中，如模型參數(shù)過多會(huì)出現(xiàn)什么情況？如何判別回歸擬合模過多會(huì)出現(xiàn)什么情況？如何判別回歸擬合模型中有過多的參數(shù)呢？型中有過多的參數(shù)呢？丙烷在氫型絲光沸石上的吸附平衡丙烷在氫型絲光沸石上的吸附平衡q例：選用不同吸附方程擬合丙烷在氫型絲光沸石體系303K的吸附平衡數(shù)據(jù)。q目的：闡明如何對(duì)模型擬合結(jié)果進(jìn)展統(tǒng)計(jì)分析，確定模型擬合的好壞、模型參數(shù)的可靠性和準(zhǔn)確性，從而進(jìn)展擬合模型的選擇。 P, kPaq, mmol/gP, kPaq, mmol/gP, kPaq, mmol/gP, kPaq, mmol/g0.100.09 1.080.4812.670.81115.891.140

14、.140.12 1.470.5116.700.85140.071.170.220.18 1.510.5324.810.90158.901.190.330.24 2.270.5934.280.95176.761.200.410.30 3.220.6443.850.98193.371.220.490.31 4.720.6954.621.02206.811.240.570.36 5.060.7065.791.040.770.41 7.390.7573.191.060.990.4410.260.7994.661.09表表8-4 303K8-4 303K時(shí)丙烷在氫型絲光沸石上的吸附平衡數(shù)據(jù)時(shí)丙烷在氫型絲

15、光沸石上的吸附平衡數(shù)據(jù) 具有代表性的、也是適用性較廣的模型具有代表性的、也是適用性較廣的模型 q1、Lanmuir (L)雙參數(shù)方程：雙參數(shù)方程： q2、Freundlich (F)雙參數(shù)方程：雙參數(shù)方程： nn aPaPm1(8.2.8) q3、BET雙參數(shù)方程：雙參數(shù)方程：q4、Langmuir-Freundlich (LF)三參數(shù)方程：三參數(shù)方程：q5、三參數(shù)方程：、三參數(shù)方程： q6、Toth三參數(shù)方程：三參數(shù)方程： q7、擴(kuò)展的、擴(kuò)展的LF方程五參數(shù)：方程五參數(shù)： q8、(14)式的特殊方式四參數(shù)式的特殊方式四參數(shù) ： naPbnn cxxxcxm()()11xP P/0nn aPa

16、Pmbb1Pbc111exp() n nm/nn PbPmcc()/1nn aPdcPaPcPmbebe1nn aPbc bd PaPcPmbdbd()11(8.2.9) (8.2.10) (8.2.11) (8.2.12) (8.2.13) (8.2.14) (8.2.15) 各模型的計(jì)算結(jié)果各模型的計(jì)算結(jié)果 Eq.(8)Eq.(9)Eq.(10)Eq.(11)Eq.(12)Eq.(13)Eq.(14)Eq.(15)nm1.0840.051/0.9760.0250.4380.0534.62317.221.5350.2570.7580.0880.7690.068a0.5530.1310.446

17、0.034/1.3820.107/1.3290.4961.4270.144b/0.200.018/0.4940.0620.6781.6580.5490.0760.9420.1070.9750.058c/812.1116.3/20.9652.650.3410.0591.9071.5930.0170.003d/0.0240.0480.9120.046e/1.6380.586/s29.09010-27.72110-25.13410-23.79910-22.9933.15410-21.04810-28.68610-3R20.987900.991270.996140.997950.998580.9985

18、90.999860.99990表表8 85 5 吸附等溫線關(guān)聯(lián)的參數(shù)值、方差和回歸系數(shù)吸附等溫線關(guān)聯(lián)的參數(shù)值、方差和回歸系數(shù) 從表中可看出，方程從表中可看出，方程(814)(814)擬合方差逐漸減少，回歸系數(shù)更擬合方差逐漸減少，回歸系數(shù)更接近接近1 1方程方程(12)(12)是經(jīng)過壓力數(shù)據(jù)來(lái)擬合的，故擬合方差和其是經(jīng)過壓力數(shù)據(jù)來(lái)擬合的，故擬合方差和其它方程的結(jié)果不是在同一數(shù)量級(jí)上它方程的結(jié)果不是在同一數(shù)量級(jí)上) )。由方程。由方程(14)(14)的五參數(shù)方的五參數(shù)方式改良的方程式改良的方程(15)(15)式獲得的結(jié)果最好，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)幾乎完全落式獲得的結(jié)果最好，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)幾乎完全落在方程在方程1

19、515式的曲線上式的曲線上( (見以下圖見以下圖) )。方程方程(13)和方程和方程(15)的擬合結(jié)果的擬合結(jié)果 0.11.010.0100.0p, kPa0.00.51.0n, mmol/gExp. dataToth Eq.(13)Eq.(15)方程方程(15)(15)式獲得的結(jié)果最好，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)幾乎完式獲得的結(jié)果最好，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)幾乎完全落在方程全落在方程1515式的曲線上。式的曲線上。判別模型參數(shù)能否過少的根據(jù)判別模型參數(shù)能否過少的根據(jù)q經(jīng)過對(duì)方程經(jīng)過對(duì)方程(13)和五參數(shù)方程和五參數(shù)方程(15)的殘差進(jìn)展分析，的殘差進(jìn)展分析，q方程方程(13)因參數(shù)過少，吸附量的計(jì)算誤差與實(shí)驗(yàn)吸附量之間

20、存在著某種分因參數(shù)過少，吸附量的計(jì)算誤差與實(shí)驗(yàn)吸附量之間存在著某種分布。布。q方程方程(15) 計(jì)算誤差在零的兩邊是隨機(jī)分布的，看不出規(guī)律性。因此，擬合計(jì)算誤差在零的兩邊是隨機(jī)分布的，看不出規(guī)律性。因此，擬合計(jì)算誤差有無(wú)規(guī)律性的分布是判別模型參數(shù)能否過少的根據(jù)。計(jì)算誤差有無(wú)規(guī)律性的分布是判別模型參數(shù)能否過少的根據(jù)。因此，擬合計(jì)算誤差有無(wú)規(guī)律性的分布是判別模型參數(shù)能因此，擬合計(jì)算誤差有無(wú)規(guī)律性的分布是判別模型參數(shù)能否過少的根據(jù)。否過少的根據(jù)。 Predicted versus Residual Valuesn, mmol/gResidual Values-0.08-0.06-0.04-0.020

21、.000.020.040.060.00.20.40.60.81.01.21.4Predicted versus Residual Valuesn, mmol/gResidual Values-0.020-0.015-0.010-0.0050.0000.0050.0100.0150.0200.00.20.40.60.81.01.21.4方程方程(13)(13)的擬合誤差的擬合誤差 方程方程(15)(15)的擬合誤差的擬合誤差 判別模型參數(shù)能否過多的根據(jù)判別模型參數(shù)能否過多的根據(jù)Eq.(14)Eq.(15)nm0.7580.0880.7690.068a1.3290.4961.4270.144b0.

22、9420.1070.9750.058c1.9071.5930.0170.003d0.0240.0480.9120.046e1.6380.586/s21.04810-28.68610-3R20.999860.99990q在方程在方程(14)的計(jì)算結(jié)果中，有的計(jì)算結(jié)果中，有些參數(shù)些參數(shù)95%的置信度較大，闡的置信度較大，闡明這些參數(shù)之間有聯(lián)絡(luò)，不是明這些參數(shù)之間有聯(lián)絡(luò)，不是獨(dú)立的。獨(dú)立的。q而對(duì)于方程而對(duì)于方程(14)的五參數(shù)方式，的五參數(shù)方式，即方程即方程(15)，其一切參數(shù)的，其一切參數(shù)的95%置信度都較小?，F(xiàn)實(shí)上，置信度都較小?，F(xiàn)實(shí)上，方程方程(15)就是據(jù)此分析對(duì)吸附就是據(jù)此分析對(duì)吸附平衡

23、實(shí)際作進(jìn)一步研討而獲得平衡實(shí)際作進(jìn)一步研討而獲得的。的。因此，擬合參數(shù)因此，擬合參數(shù)95%95%的置信度能否較大是判別模型參數(shù)能否的置信度能否較大是判別模型參數(shù)能否過多的根據(jù)。過多的根據(jù)?；貧w模型的選擇總結(jié)回歸模型的選擇總結(jié)q模型參數(shù)較少時(shí)會(huì)出現(xiàn)擬合殘差的分布不是模型參數(shù)較少時(shí)會(huì)出現(xiàn)擬合殘差的分布不是隨機(jī)的，而是呈現(xiàn)某種分布，相互關(guān)聯(lián)。殘隨機(jī)的，而是呈現(xiàn)某種分布，相互關(guān)聯(lián)。殘差的分布偏離正態(tài)分布較遠(yuǎn)。差的分布偏離正態(tài)分布較遠(yuǎn)。 q模型參數(shù)過多會(huì)出現(xiàn)某些參數(shù)模型參數(shù)過多會(huì)出現(xiàn)某些參數(shù)95%的置信度的置信度較大，闡明這些參數(shù)之間有聯(lián)絡(luò)，不是獨(dú)立較大，闡明這些參數(shù)之間有聯(lián)絡(luò)，不是獨(dú)立的。的。習(xí)題習(xí)題

24、q研討二硫化碳飽和蒸汽壓的模型回歸問題。P266，Ex8.3二硫化碳的根本性質(zhì)：二硫化碳的根本性質(zhì)： 臨界溫度為臨界溫度為273.05273.05 臨界壓力為臨界壓力為72.87 atm72.87 atm。 溫度，蒸汽壓，mmHg溫度，蒸汽壓，mmHg-701.610198.0-603.520297.5-507.130432.7-4014.040616.7-3026.250995.6-2046.5601170.4-1078.8701558.00127.3Statistica的非線性估計(jì)的非線性估計(jì)q非線性估計(jì)方法非線性估計(jì)方法qUser-Specified Regression, least

25、square 可以計(jì)算95%置信區(qū)間“l(fā)east square與與“Custom Loss 比較比較q“l(fā)east square計(jì)算結(jié)果采用計(jì)算結(jié)果采用Levenberg-Marquardt方法方法Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qLevel of confidence: 95.0% ( alpha=0.050)EstimateStandarderrort-valuedf = 28p-levelLo. ConfLimitUp. ConfLimitnmabc

26、d0.7637600.032019 23.85354 0.000000 0.698173 0.8293471.4359130.069004 20.80900 0.000000 1.294564 1.5772620.9822270.028027 35.04555 0.000000 0.924816 1.0396380.0163330.0016779.73910 0.000000 0.012898 0.0197680.9168660.022152 41.38886 0.000000 0.871489 0.962243Model: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p

27、*b*(1+c*p*d). (propane-mor.sta)Dep. var: q Loss: (OBS-PRED)*2Final loss: .002032166 R=.99975 Variance explained: 99.950%N=33nmabcdEstimate0.763768 1.435910 0.982209 0.016336 0.916840q“Custom Loss 計(jì)算結(jié)果采用計(jì)算結(jié)果采用Quasi-Newton法法qMatrix ill conditioned; cannot compute standard errors.EConf“Custom Loss 的方差分

28、析與迭代步的方差分析與迭代步驟驟 “Custom Custom Loss Loss 無(wú)迭無(wú)迭代歷史紀(jì)錄，代歷史紀(jì)錄，協(xié)方差分析結(jié)協(xié)方差分析結(jié)果已出現(xiàn)病態(tài)。果已出現(xiàn)病態(tài)。Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qEffect1Sum of Sqares2DF3Mean Squares4F-value5p-valueRegressionResidualTotalCorrected TotalRegression vs.Corrected Total21.505475.

29、000004.30109459262.190.0000000.0020328.000000.00007321.5075033.000004.0675032.0000021.505475.000004.30109433.840.000000Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qLossFunctionnmabcd123456789104.629683 0.100000 0.100000 0.100000 0.100000 0.1000004.294027 0.1

30、84348 0.196978 0.173970 0.528173 0.5570214.227658 0.191476 0.221943 0.188990 0.587318 0.5716294.102290 0.204036 0.266635 0.217677 0.697367 0.5994433.874137 0.224709 0.339318 0.270956 0.885848 0.6511153.429648 0.259087 0.453795 0.373779 1.198671 0.7516122.556049 0.312350 0.616018 0.561417 1.656081 0.

31、9378231.036472 0.392841 0.756563 0.881424 1.932017 1.2618210.899994 0.404754 0.764595 0.929852 1.734310 1.3103860.755383 0.425557 0.992230 0.962118 0.684109 1.321683Covariance matrix cannot be computed. Statistica非線性估計(jì)的殘差分析非線性估計(jì)的殘差分析殘差分布情況殘差分布情況殘差對(duì)預(yù)測(cè)值作圖殘差對(duì)預(yù)測(cè)值作圖對(duì)方程對(duì)方程(15)的殘差分析的殘差分析Predicted versus Re

32、sidual Values123456789101112131415161718192021222324252627282930313233-0.20.00.20.40.60.81.01.21.4Predicted Values-0.020-0.015-0.010-0.0050.0000.0050.0100.0150.020Residual ValuesFrequency Distribution:-0.020-0.015-0.010-0.0050.0000.0050.0100.0150.020Residual Values01234567891011No of obsHistogram of

33、 Histogram of residualsresidualsResidual vs. Residual vs. PredictedPredicted誤差正態(tài)分布圖誤差正態(tài)分布圖Frequency Distribution: Residuals Expected Normal-0.0030-0.0025-0.0020-0.0015-0.0010-0.00050.00000.00050.00100.00150.00200.00250.00300.003502468101214161820222426No of obsFrequency Distribution:-0.0030-0.0025-0

34、.0020-0.0015-0.0010-0.00050.00000.00050.00100.00150.00200.00250.00300.00350.0040Residual Values0510152025303540No of obs“l(fā)east square計(jì)算的誤差正態(tài)分布圖計(jì)算的誤差正態(tài)分布圖“Custom Loss 計(jì)算的誤差正態(tài)分布圖計(jì)算的誤差正態(tài)分布圖qAntoine 方程擬合結(jié)果方程擬合結(jié)果非線性函數(shù)的管理非線性函數(shù)的管理三、三、MATLAB的擬合函數(shù)的擬合函數(shù)q多項(xiàng)式擬合函數(shù)多項(xiàng)式擬合函數(shù)polyfitq非線性最小二乘法非線性最小二乘法qlsqnonlin()非線性最小二

35、乘優(yōu)化問題非線性最小二乘優(yōu)化問題qlsqcurvefit()非線性最小二乘曲線擬合非線性最小二乘曲線擬合qnlinfit()前兩種的簡(jiǎn)化版本前兩種的簡(jiǎn)化版本qnlparci()計(jì)算參數(shù)的置信區(qū)間計(jì)算參數(shù)的置信區(qū)間qnlpredci()計(jì)算預(yù)測(cè)值的置信區(qū)間計(jì)算預(yù)測(cè)值的置信區(qū)間qnlintool()nlinfit()、nlparci()、nlpredci()的集成的集成圖形用戶界面擬合函數(shù)圖形用戶界面擬合函數(shù)q留意：不同的擬合函數(shù)命令，其優(yōu)化目的函數(shù)定義留意：不同的擬合函數(shù)命令，其優(yōu)化目的函數(shù)定義以及調(diào)用方式不同，留意區(qū)分！以及調(diào)用方式不同，留意區(qū)分！一一Polyfit的運(yùn)用的運(yùn)用qp = pol

36、yfit(x0, y0, n)q其中其中x0和和y0分別為察看節(jié)點(diǎn)和察看值向量；分別為察看節(jié)點(diǎn)和察看值向量；n表示插值多項(xiàng)式表示插值多項(xiàng)式的次數(shù)；輸出值的次數(shù)；輸出值p表示插值多項(xiàng)式的系數(shù)。表示插值多項(xiàng)式的系數(shù)。q例：某實(shí)驗(yàn)中測(cè)得一組數(shù)據(jù)，其值如下：例：某實(shí)驗(yàn)中測(cè)得一組數(shù)據(jù)，其值如下：q x l 2 3 4 5q y 1.3 1.8 2.2 2.9 3.5q 知知x和和y成線性關(guān)系，即成線性關(guān)系，即y=kx+b，求系數(shù)，求系數(shù)k和和b x=1 2 3 4 5;y=1.3 1.8 2.2 2.9 3.5;p=polyfit(x,y,1)y1=polyval(p,x);plot(x,y1)hold

37、 on;plot(x,y,b*);p = 0.55 0.6911.522.533.544.5511.522.533.5也就是，也就是，k為為0.55，b為為0.69二二lsqcurvefit的運(yùn)用的運(yùn)用q方程目的函數(shù)qFind coefficients beta that best fit the equationq調(diào)用方式調(diào)用方式qbeta= lsqcurvefit(fun, beta0,xdata,ydata) qbeta= lsqcurvefit(fun, beta0,xdata,ydata,lb,ub) qbeta = lsqcurvefit(fun, beta0,xdata,ydat

38、a,lb,ub,options) xdata和和ydata分別為察看節(jié)點(diǎn)和察看值向量；分別為察看節(jié)點(diǎn)和察看值向量；fun自定義的非線性擬合模型；自定義的非線性擬合模型；beta0擬合參數(shù)的初始值；擬合參數(shù)的初始值；beta擬合模型中的參數(shù)；擬合模型中的參數(shù)；lb,ub擬合參數(shù)初值的邊境值擬合參數(shù)初值的邊境值, lb = beta 1 % two output arguments J = . % Jacobian of the function evaluated at betaend options = optimset(Jacobian,on)qFor examplefunction F =

39、 myfun(beta,xdata)F = beta(1)*exp(beta(2)*xdata);lsqcurvefit運(yùn)用例如運(yùn)用例如% Assume you determined xdata and ydata experimentallyxdata = 0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3;ydata = 455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5;beta0 = 100; -1 % Starting guessbeta,resnorm = lsqcurvefit(my

40、fun, beta0,xdata,ydata)function F = myfun(beta,xdata)F = beta (1)*exp(beta(2)*xdata);beta = 498.8309 -0.1013 resnorm = 9.504921b xybe三三lsqnonlin的運(yùn)用的運(yùn)用q方程目的函數(shù)qFind coefficients beta that best fit the equationq調(diào)用方式調(diào)用方式qbeta = lsqnonlin(fun, beta 0) qbeta = lsqnonlin(fun, beta 0,lb,ub) qbeta = lsqnonli

41、n(fun, beta 0,lb,ub,options) fun自定義的優(yōu)化函數(shù)；自定義的優(yōu)化函數(shù)；beta 0優(yōu)化參數(shù)的初始值；優(yōu)化參數(shù)的初始值；beta擬合模型中的參數(shù)；擬合模型中的參數(shù)；lb,ub擬合參數(shù)初值的邊境值擬合參數(shù)初值的邊境值, lb = beta s=dsolve(Dy1=-k1*y1, Dy2=k1*y1-k2*y2, y1(0)=1, y2(0)=0) s = y1: 1x1 sym y2: 1x1 sym s.y1 s.y2 simplify(s.y2) ans = exp(-k1*t) ans = -(-k1+k2)/(k1-k2)2*k1*exp(-k2*t)-1/

42、(k1-k2)*k1*exp(-k1*t)ans = k1*(exp(-k2*t)-exp(-k1*t)/(k1-k2)tkAeC1)(21121tktkBeekkkC對(duì)積分式進(jìn)展參數(shù)估計(jì)的源程序?qū)Ψe分式進(jìn)展參數(shù)估計(jì)的源程序 function seqcurvefit11clear all;load seqdata;beta0=0.005 0.001; lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(seqfun,beta0,lb,ub,t,c); ci = nlparci

43、(beta,residual,jacobian); function y = seqfun(beta,t,c)ca=exp(-beta(1)*t);cb=beta(1)/(beta(2)-beta(1)*(exp(-beta(1)*t)-exp(-beta(2)*t);y=ca-c(:,1) cb-c(:,2); 輸出計(jì)算結(jié)果的源程序輸出計(jì)算結(jié)果的源程序 % print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1 = %.4f %.4fn,beta(1),ci(1,2)-beta(1)fprintf(t

44、k2 = %.4f %.4fn,beta(2),ci(2,2)-beta(2)fprintf( The sum of the squares is: %.1enn,sum(residual.2)% plot of fit resultstc=linspace(0,max(t),200);y_row,y_col=size(c);zeroc=zeros(200,y_col);yc = seqfun(beta,tc,zeroc);plot(t,c(:,1),ro,tc,yc(:,1),r-);hold onplot(t,c(:,2),b+,tc,yc(:,2),b-);xlabel(Time);yl

45、abel(Concentration);hold off用函數(shù)求擬合值。用函數(shù)求擬合值。留意轉(zhuǎn)置，與函數(shù)定義留意轉(zhuǎn)置，與函數(shù)定義的矩陣維數(shù)一致！的矩陣維數(shù)一致！參數(shù)擬合結(jié)果參數(shù)擬合結(jié)果 seqcurvefit11Optimization terminated: relative function value changing by less than OPTIONS.TolFun. Estimated Parameters by Lsqnonlin(): k1 = 0.0412 0.0018 k2 = 0.0121 0.0008 The sum of the residual squares

46、is: 2.7e-00302040608010012014016018020000.10.20.30.40.50.60.70.80.91TimeConcentration方法方法2：對(duì)微分式進(jìn)展參數(shù)估計(jì)：對(duì)微分式進(jìn)展參數(shù)估計(jì)function seqcurvefit21clear all;load seqdata;y_row,y_col=size(c); beta0=0.005 0.001;c0=1 0;lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(seqfun,beta

47、0,lb,ub,t,c,y_col,c0);ci = nlparci(beta,residual,jacobian); 擬合模型和目的函數(shù)的定義擬合模型和目的函數(shù)的定義function y = seqfun(beta,t,c,y_col,c0) % Objective functiontspan = 0 max(t);tt yy = ode45(modeleqs,tspan,c0,beta); for col = 1:y_col yc(:,col) = spline(tt,yy(:,col),t);endy=c(:,1)-yc(:,1); c(:,2)-yc(:,2); function dy

48、dt = modeleqs(t,y,beta) % Model equationdydt = -beta(1)*y(1); beta(1)*y(1)-beta(2)*y(2); beta,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(seqfun,beta0,lb,ub,t,c,y_col,c0);輸出計(jì)算結(jié)果的源程序輸出計(jì)算結(jié)果的源程序 % print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1 = %.4f %.4fn,b

49、eta(1),ci(1,2)-beta(1)fprintf(t k2 = %.4f %.4fn,beta(2),ci(2,2)-beta(2)fprintf( The sum of the squares is: %.1enn,sum(residual.2)% plot of fit resultstspan = 0 max(t);tt yc = ode45(modeleqs,tspan,c0,beta);tc=linspace(0,max(t),200);yca = spline(tt,yc(:,1),tc);plot(t,c(:,1),ro,tc,yca,r-);hold onycb =

50、spline(tt,yc(:,2),tc);plot(t,c(:,2),b+,tc,ycb,b-);xlabel(Time);ylabel(Concentration);hold off 首先解微分方程求擬合首先解微分方程求擬合值。然后用樣條插值逼值。然后用樣條插值逼近。近。參數(shù)擬合結(jié)果參數(shù)擬合結(jié)果 seqcurvefit21Optimization terminated: relative function value changing by less than OPTIONS.TolFun. Estimated Parameters by Lsqnonlin(): k1 = 0.0412

51、0.0018 k2 = 0.0121 0.0008 The sum of the residual squares is: 5.0e-003 02040608010012014016018020000.10.20.30.40.50.60.70.80.91TimeConcentration運(yùn)用運(yùn)用MATLAB進(jìn)展參數(shù)估計(jì)例如進(jìn)展參數(shù)估計(jì)例如q例2：青霉素發(fā)酵過程動(dòng)力學(xué)參數(shù)估計(jì)： q 在間歇發(fā)酵罐中研討青霉素發(fā)酵過程動(dòng)力學(xué)，微生物Penicillium chrysogenum在一定的控制條件下生長(zhǎng)，細(xì)胞生長(zhǎng)速率可以用邏輯模型描畫：q 青霉素的消費(fèi)速率模型為： q初始條件為：t=0， y1=0.29

52、，y2=0 。211111kyykdtdy24132ykykdtdyTimehoursCellconcentrationdry weightPenicillinconcentrationunits/ml00.180100.120220.480.0089341.460.0642461.560.2266581.730.4373701.990.6943822.621.2459942.881.43151063.432.04021183.371.92781303.922.18481423.962.42041543.582.46151663.582.2831783.342.70781903.472.654

53、2運(yùn)用運(yùn)用MATLAB的的dsolve求積分求積分 s=dsolve(Dy1=k1*y1*(1-y1/k1), Dy2=k3*y1-k4*y2,y1(0)=0.29,y2(0)=0) s = y1: 1x1 sym y2: 1x1 sym s.y1 s.y2 simplify(s.y2) ans = k1/(1+1/29*exp(-k1*t)*(100*k1-29) ans = Int(k3*k1/(1+1/29*exp(-k1*_z1)*(100*k1-29)*exp(k4*_z1),_z1 = 0 . t)*exp(-k4*t) ans = 29*k3*k1*Int(exp(k4*_z1)/(29+100*exp(-k1*_z1)*k1-29*exp(-k1*_z1),_z1 = 0 . t)*exp(-k4*t) 沒有確定的積分表達(dá)式?jīng)]有確定的積分表達(dá)式參數(shù)估計(jì)的源程序參數(shù)估計(jì)的源程序function PenicilliumEstclear all;load PenicilliumData % Load experimental data y_row,y_col=size(y); % Nonlinear least square estimate using lsqnonlin()beta0 = 0.1 4.0