8建模與仿真統(tǒng)計回歸模型

1、1/332021-11-5 牙膏的銷售量牙膏的銷售量 軟件開發(fā)人員的薪金軟件開發(fā)人員的薪金酶促反應(yīng)酶促反應(yīng)第十章第十章 統(tǒng)計回歸模型統(tǒng)計回歸模型2/332021-11-5回歸模型是用回歸模型是用統(tǒng)計分析方法建立的最常用的一類模型統(tǒng)計分析方法建立的最常用的一類模型 數(shù)學(xué)建模的基本方法數(shù)學(xué)建模的基本方法機(jī)理分析機(jī)理分析測試分析測試分析通過對數(shù)據(jù)的統(tǒng)計分析，找出與數(shù)據(jù)擬合最好的模型通過對數(shù)據(jù)的統(tǒng)計分析，找出與數(shù)據(jù)擬合最好的模型 不涉及回歸分析的數(shù)學(xué)原理和方法不涉及回歸分析的數(shù)學(xué)原理和方法 通過實例討論如何選擇不同類型的模型通過實例討論如何選擇不同類型的模型 對軟件得到的結(jié)果進(jìn)行分析，對模型進(jìn)行改進(jìn)對

2、軟件得到的結(jié)果進(jìn)行分析，對模型進(jìn)行改進(jìn) 由于客觀事物內(nèi)部規(guī)律的復(fù)雜及人們認(rèn)識程度的限制由于客觀事物內(nèi)部規(guī)律的復(fù)雜及人們認(rèn)識程度的限制,無法分析實際對象內(nèi)在的因果關(guān)系，建立合乎機(jī)理規(guī)無法分析實際對象內(nèi)在的因果關(guān)系，建立合乎機(jī)理規(guī)律的數(shù)學(xué)模型。律的數(shù)學(xué)模型。 3/332021-11-510.1 牙膏的銷售量牙膏的銷售量 問問題題建立牙膏銷售量與價格、廣告投入之間的模型建立牙膏銷售量與價格、廣告投入之間的模型 預(yù)測在不同價格和廣告費(fèi)用下的牙膏銷售量預(yù)測在不同價格和廣告費(fèi)用下的牙膏銷售量 收集了收集了30個銷售周期本公司牙膏銷售量、價格、個銷售周期本公司牙膏銷售量、價格、廣告費(fèi)用，及同期其它廠家同類牙

3、膏的平均售價廣告費(fèi)用，及同期其它廠家同類牙膏的平均售價 9.260.556.804.253.70307.930.055.803.853.80298.510.256.754.003.7527.38-0.055.503.803.851銷售量銷售量(百萬支百萬支)價格差價格差（元）（元）廣告費(fèi)用廣告費(fèi)用(百萬元百萬元)其它廠家其它廠家價格價格(元元)本公司價本公司價格格(元元)銷售銷售周期周期4/332021-11-5基本模型基本模型y 公司牙膏銷售量公司牙膏銷售量x1其它廠家與本公司其它廠家與本公司價格差價格差x2公司廣告費(fèi)用公司廣告費(fèi)用110 xy222210 xxy55.566.577.577

4、.588.599.510 x2y-0.200.20.40.677.588.599.510 x1y22322110 xxxyx1, x2解釋變量解釋變量(回歸變量回歸變量, 自變量自變量) y被解釋變量（因變量）被解釋變量（因變量） 0, 1 , 2 , 3 回歸系數(shù)回歸系數(shù) 隨機(jī)隨機(jī)誤差（誤差（均值為零的均值為零的正態(tài)分布隨機(jī)變量）正態(tài)分布隨機(jī)變量）5/332021-11-5matlab 統(tǒng)計工具箱統(tǒng)計工具箱 模型求解模型求解b,bint,r,rint,stats=regress(y,x,alpha) 輸入輸入 x= n 4數(shù)數(shù)據(jù)矩陣據(jù)矩陣, 第第1列為全列為全1向量向量1 2221xxxal

5、pha(置信置信水平水平,0.05) 22322110 xxxyb 的的估計值估計值 bintb的置信區(qū)間的置信區(qū)間 r 殘差向量殘差向量y-xb rintr的置信區(qū)間的置信區(qū)間 stats檢驗統(tǒng)計量檢驗統(tǒng)計量 r2,f, p yn維數(shù)據(jù)向量維數(shù)據(jù)向量輸出輸出 由數(shù)據(jù)由數(shù)據(jù) y,x1,x2估計估計 參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間17.32445.7282 28.92061.30700.6829 1.9311 -3.6956-7.4989 0.1077 0.34860.0379 0.6594 r2=0.9054 f=82.9409 p=0.0000 0 1 2 36/332021-

6、11-5結(jié)果分析結(jié)果分析y的的90.54%可由模型確定可由模型確定 參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間17.32445.7282 28.92061.30700.6829 1.9311 -3.6956-7.4989 0.1077 0.34860.0379 0.6594 r2=0.9054 f=82.9409 p=0.0000 0 1 2 322322110 xxxyf遠(yuǎn)超過遠(yuǎn)超過f檢驗的臨界值檢驗的臨界值 p遠(yuǎn)小于遠(yuǎn)小于 =0.05 2的置信區(qū)間包含零點(diǎn)的置信區(qū)間包含零點(diǎn)(右端點(diǎn)距零點(diǎn)很近右端點(diǎn)距零點(diǎn)很近) x2對因變量對因變量y 的的影響不太顯著影響不太顯著x22項顯著項顯著 可將可

7、將x2保留在模型中保留在模型中 模型從整體上看成立模型從整體上看成立7/332021-11-522322110 xxxy銷售量預(yù)測銷售量預(yù)測 價格差價格差x1=其它廠家其它廠家價格價格x3-本公司本公司價格價格x4估計估計x3調(diào)整調(diào)整x4控制價格差控制價格差x1=0.2元，投入廣告費(fèi)元，投入廣告費(fèi)x2=650萬元萬元銷售量預(yù)測區(qū)間為銷售量預(yù)測區(qū)間為 7.8230，8.7636（置信度（置信度95%）上限用作庫存管理的目標(biāo)值上限用作庫存管理的目標(biāo)值 下限用來把握公司的現(xiàn)金流下限用來把握公司的現(xiàn)金流 若估計若估計x3=3.9，設(shè)定，設(shè)定x4=3.7，則可以，則可以95%的把握的把握知道銷售額在知道

8、銷售額在 7.8320 3.7 29（百萬元）以上（百萬元）以上控制控制x1通過通過x1, x2預(yù)測預(yù)測y2933.822322110 xxxy(百萬支百萬支)8/332021-11-5模型改進(jìn)模型改進(jìn)x1和和x2對對y的的影響?yīng)毩⒂绊應(yīng)毩?22322110 xxxy21422322110 xxxxxy參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間17.32445.7282 28.92061.30700.6829 1.9311 -3.6956-7.4989 0.1077 0.34860.0379 0.6594 r2=0.9054 f=82.9409 p=0.0000 0 1 2 3參數(shù)參數(shù)參數(shù)估

9、計值參數(shù)估計值置信區(qū)間置信區(qū)間29.113313.7013 44.525211.13421.9778 20.2906 -7.6080-12.6932 -2.5228 0.67120.2538 1.0887 -1.4777-2.8518 -0.1037 r2=0.9209 f=72.7771 p=0.0000 3 0 1 2 4x1和和x2對對y的影響有的影響有交互作用交互作用9/332021-11-5兩模型銷售量預(yù)測兩模型銷售量預(yù)測比較比較21422322110 xxxxxy22322110 xxxy2933. 8 y(百萬支百萬支)區(qū)間區(qū)間 7.8230，8.7636區(qū)間區(qū)間 7.8953，

10、8.7592 3272. 8 y(百萬支百萬支)控制價格差控制價格差x1=0.2元，投入廣告費(fèi)元，投入廣告費(fèi)x2=6.5百萬元百萬元預(yù)測區(qū)間長度更短預(yù)測區(qū)間長度更短 略有增加略有增加 y 10/332021-11-5x2=6.5x1=0.2 -0.200.20.40.67.588.59x1y -0.200.20.40.67.588.59x1y 56787.588.599.510 x2y 567888.599.51010.5x2y 22322110 xxxy21422322110 xxxxxy兩模型兩模型 與與x1, ,x2關(guān)系的關(guān)系的比較比較y 11/332021-11-5交互作用影響的討論交

11、互作用影響的討論2221 . 06712. 07558. 72267.301xxyx價格差價格差 x1=0.1 價格差價格差 x1=0.32223 . 06712. 00513. 84535.321xxyx21422322110 xxxxxy5357. 72x加大廣告投入使銷售量增加加大廣告投入使銷售量增加 （ x2大于大于6百萬元）百萬元）價格差較小時增加價格差較小時增加的速率更大的速率更大 56787.588.599.51010.5x1=0.1x1=0.3x2y 1 . 03 . 011xxyy價格優(yōu)勢會使銷售量增加價格優(yōu)勢會使銷售量增加 價格差較小時更需要靠廣告價格差較小時更需要靠廣告來

12、吸引顧客的眼球來吸引顧客的眼球 12/332021-11-5完全二次多項式模型完全二次多項式模型 22521421322110 xxxxxxymatlab中有命令中有命令rstool直接求解直接求解00.20.47.588.599.5105.566.57x1x2y ),(543210從輸出從輸出 export 可得可得13/332021-11-5function y=ill(t,x)a=1;b=0.3;y=a*x(1)*x(2)-b*x(1),-a*x(1)*x(2);ts=0:50;x0=0.02,0.98;t,x=ode45(ill,ts,x0);t,xplot(t,x(:,1),t,x(

13、:,2);gridfigureplot(x(:,2),x(:,1),grid,運(yùn)行結(jié)果如下：附附matlab部分程序部分程序 14/332021-11-515/332021-11-516/332021-11-510.2 軟件開發(fā)人員的薪金軟件開發(fā)人員的薪金資歷資歷 從事專業(yè)工作的年數(shù)；管理從事專業(yè)工作的年數(shù)；管理 1= =管理人員，管理人員，0= =非管理人非管理人員；教育員；教育 1= =中學(xué)，中學(xué)，2= =大學(xué)，大學(xué)，3= =更高程度更高程度建立模型研究薪金與資歷、管理責(zé)任、教育程度的關(guān)系建立模型研究薪金與資歷、管理責(zé)任、教育程度的關(guān)系分析人事策略的合理性，作為新聘用人員薪金的參考分析人事

14、策略的合理性，作為新聘用人員薪金的參考 編編號號薪金薪金資資歷歷管管理理教教育育0113876111021160810303187011130411283102編編號號薪金薪金資資歷歷管管理理教教育育422783716124318838160244174831601451920717024619346200146名軟件開發(fā)人員的檔案資料名軟件開發(fā)人員的檔案資料 17/332021-11-5分析與假設(shè)分析與假設(shè) y 薪金，薪金，x1 資歷（年）資歷（年）x2 = = 1 管理人員，管理人員，x2 = = 0 非管理人員非管理人員1= =中學(xué)中學(xué)2= =大學(xué)大學(xué)3= =更高更高其它中學(xué),x013其

15、它大學(xué),x014資歷每加一年薪金的增長是常數(shù)；資歷每加一年薪金的增長是常數(shù)；管理、教育、資歷之間無交互作用管理、教育、資歷之間無交互作用 教教育育443322110 xaxaxaxaay線性回歸模型線性回歸模型 a0, a1, , a4是待估計的回歸系數(shù)，是待估計的回歸系數(shù)， 是隨機(jī)誤差是隨機(jī)誤差 中學(xué)：中學(xué)：x3=1, x4=0 ；大大學(xué)：學(xué)：x3=0, x4=1； 更高：更高：x3=0, x4=0 18/332021-11-5模型求解模型求解443322110 xaxaxaxaay參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間a011032 10258 11807 a1546 484 608

16、 a26883 6248 7517 a3-2994 -3826 -2162 a4148 -636 931 r2=0.957 f=226 p=0.000r2,f, p 模型整體上可用模型整體上可用資歷增加資歷增加1年薪年薪金增長金增長546 管理人員薪金多管理人員薪金多6883 中學(xué)程度薪金比更中學(xué)程度薪金比更高的少高的少2994 大學(xué)程度薪金比更大學(xué)程度薪金比更高的多高的多148 a4置信區(qū)間包含零點(diǎn)，置信區(qū)間包含零點(diǎn)，解釋不可靠解釋不可靠! !中學(xué)：中學(xué)：x3=1, x4=0;大大學(xué)：學(xué)：x3=0, x4=1; 更高：更高：x3=0, x4=0. x2 = = 1 管理，管理，x2 = =

17、0 非管理非管理x1資歷資歷( (年年) )19/332021-11-5殘差分析方法殘差分析方法 結(jié)果分析結(jié)果分析443322110 xaxaxaxaay殘差殘差yyee 與資歷與資歷x1的關(guān)系的關(guān)系 05101520-2000-1000010002000e與管理與管理教育組合的關(guān)系教育組合的關(guān)系 123456-2000-1000010002000殘差全為正，或全為負(fù)，管理殘差全為正，或全為負(fù)，管理教教育組合處理不當(dāng)育組合處理不當(dāng) 殘差大概分成殘差大概分成3個水平，個水平， 6種管理種管理教育組合混在教育組合混在一起，未正確反映一起，未正確反映 。應(yīng)在模型中增加管理應(yīng)在模型中增加管理x2與教育

18、與教育x3, x4的交互項的交互項 組合組合123456管理管理010101教育教育112233管理與教育的組合管理與教育的組合20/332021-11-5426325443322110 xxaxxaxaxaxaxaay進(jìn)一步的模型進(jìn)一步的模型增加管理增加管理x2與教育與教育x3, x4的交互項的交互項參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間a01120411044 11363a1497486 508a270486841 7255a3-1727-1939 -1514a4-348-545 152a5-3071-3372 -2769a618361571 2101r2=0.999 f=554 p

19、=0.000r2, ,f有改進(jìn)，所有回歸系數(shù)置信有改進(jìn)，所有回歸系數(shù)置信區(qū)間都不含零點(diǎn)，模型完全可用區(qū)間都不含零點(diǎn)，模型完全可用 消除了不正?，F(xiàn)象消除了不正常現(xiàn)象 異常數(shù)據(jù)異常數(shù)據(jù)( (33號號) )應(yīng)去掉應(yīng)去掉 05101520-1000-5000500e x1 123456-1000-5000500e 組合組合21/332021-11-5去掉異常數(shù)據(jù)后去掉異常數(shù)據(jù)后的結(jié)果的結(jié)果參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間a01120011139 11261a1498494 503a270416962 7120a3-1737-1818 -1656a4-356-431 281a5-3056-3

20、171 2942a619971894 2100r2= 0.9998 f=36701 p=0.000005101520-200-1000100200e x1 123456-200-1000100200e 組合組合r2： 0.957 0.999 0.9998f： 226 554 36701 置信區(qū)間長度更短置信區(qū)間長度更短殘差殘差圖十分正常圖十分正常最終模型的結(jié)果可以應(yīng)最終模型的結(jié)果可以應(yīng)用用22/332021-11-5模型應(yīng)用模型應(yīng)用 制訂制訂6種管理種管理教育組合人員的教育組合人員的“基礎(chǔ)基礎(chǔ)”薪金薪金( (資歷為資歷為0）組合組合管理管理教育教育系數(shù)系數(shù)“基礎(chǔ)基礎(chǔ)”薪金薪金101a0+a39

21、463211a0+a2+a3+a513448302a0+a410844412a0+a2+a4+a619882503a011200613a0+a218241426325443322110 xxaxxaxaxaxaxaay中學(xué)：中學(xué)：x3=1, x4=0 ；大；大學(xué)：學(xué)：x3=0, x4=1； 更高：更高：x3=0, x4=0 x1= = 0； x2 = = 1 管理，管理，x2 = = 0 非管理非管理大學(xué)程度管理人員比大學(xué)程度管理人員比更高更高程度管理人員的薪金高程度管理人員的薪金高 大學(xué)程度非管理人員比大學(xué)程度非管理人員比更高更高程度非管理人員的薪金略低程度非管理人員的薪金略低 23/332

22、021-11-5對定性因素對定性因素( (如管理、教育如管理、教育) )，可以，可以引入引入0- -1變量變量處理，處理，0- -1變量的個數(shù)應(yīng)比定性因素的水平少變量的個數(shù)應(yīng)比定性因素的水平少1 軟件開發(fā)人員的薪金軟件開發(fā)人員的薪金殘差分析方法殘差分析方法可以發(fā)現(xiàn)模型的缺陷，可以發(fā)現(xiàn)模型的缺陷，引入交互作用項引入交互作用項常常能夠改善模型常常能夠改善模型 剔除異常數(shù)據(jù)剔除異常數(shù)據(jù)，有助于得到更好的結(jié)果，有助于得到更好的結(jié)果注：可以直接對注：可以直接對6種管理種管理教育組合引入教育組合引入5個個0- -1變量變量 24/332021-11-510.3 酶促反應(yīng)酶促反應(yīng) 問問題題研究酶促反應(yīng)（研究

23、酶促反應(yīng)（酶催化反應(yīng)）酶催化反應(yīng)）中嘌呤霉素對反中嘌呤霉素對反應(yīng)速度與底物應(yīng)速度與底物（反應(yīng)物）（反應(yīng)物）濃度之間關(guān)系的影響濃度之間關(guān)系的影響 建立數(shù)學(xué)模型，反映該酶促反應(yīng)的速度與底建立數(shù)學(xué)模型，反映該酶促反應(yīng)的速度與底物濃度以及經(jīng)嘌呤霉素處理與否之間的關(guān)系物濃度以及經(jīng)嘌呤霉素處理與否之間的關(guān)系 設(shè)計了兩個實驗設(shè)計了兩個實驗 ：酶經(jīng)過嘌呤霉素處理；酶未：酶經(jīng)過嘌呤霉素處理；酶未經(jīng)嘌呤霉素處理。實驗數(shù)據(jù)見下表經(jīng)嘌呤霉素處理。實驗數(shù)據(jù)見下表: 方方案案底物濃度底物濃度(ppm)0.020.060.110.220.561.10反應(yīng)反應(yīng)速度速度處理處理7647971071231391591521912

24、01207200未處理未處理6751848698115131124144158160/25/332021-11-5線性化模型線性化模型 經(jīng)嘌呤霉素處理后實驗數(shù)據(jù)的估計結(jié)果經(jīng)嘌呤霉素處理后實驗數(shù)據(jù)的估計結(jié)果 參數(shù)參數(shù)參數(shù)估計值（參數(shù)估計值（10-3）置信區(qū)間（置信區(qū)間（10-3） 15.1073.539 6.676 20.2470.176 0.319r2=0.8557 f=59.2975 p=0.00008027.195/11104841. 0/122xxy21xy111121對對 1 , 2非線性非線性 對對 1, 2線性線性 x12126/332021-11-5線性化模型結(jié)果分析線性化模型結(jié)

25、果分析 x較大時，較大時，y有較大偏差有較大偏差 1/x較小時有很好的較小時有很好的線性趨勢，線性趨勢，1/x較大較大時出現(xiàn)很大的起落時出現(xiàn)很大的起落 參數(shù)估計時，參數(shù)估計時，x較小較?。?/x很大）的數(shù)據(jù)控很大）的數(shù)據(jù)控制了回歸參數(shù)的確定制了回歸參數(shù)的確定 0102030405000.0050.010.0150.020.0251/y1/xxy112100.511.5050100150200250 xxy21xy27/332021-11-5beta,r,j = nlinfit (x,y,model,beta0) beta的置信區(qū)間的置信區(qū)間matlab 統(tǒng)計工具箱統(tǒng)計工具箱 輸入輸入 x自變量

26、自變量數(shù)據(jù)矩陣數(shù)據(jù)矩陣y 因變量數(shù)據(jù)向量因變量數(shù)據(jù)向量beta 參數(shù)的估計值參數(shù)的估計值r 殘差，殘差，j 估計預(yù)估計預(yù)測誤差的測誤差的jacobi矩陣矩陣 model 模型的函數(shù)模型的函數(shù)m文件名文件名beta0 給定的參數(shù)初值給定的參數(shù)初值 輸出輸出 betaci =nlparci(beta,r,j) 非線性模型參數(shù)估計非線性模型參數(shù)估計function y=f1(beta, x)y=beta(1)*x./(beta(2)+x);xxy21x= ; y= ;beta0=195.8027 0.04841;beta,r,j=nlinfit(x,y,f1,beta0);betaci=nlparc

27、i(beta,r,j);beta, betaci beta0線性化線性化模型估計結(jié)果模型估計結(jié)果 28/332021-11-5非線性模型結(jié)果分析非線性模型結(jié)果分析參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間 1212.6819197.2029 228.1609 20.06410.0457 0.0826 畫面左下方的畫面左下方的export 輸出其它統(tǒng)計結(jié)果。輸出其它統(tǒng)計結(jié)果。拖動畫面的十字線，得拖動畫面的十字線，得y的預(yù)測值和預(yù)測區(qū)間的預(yù)測值和預(yù)測區(qū)間剩余標(biāo)準(zhǔn)差剩余標(biāo)準(zhǔn)差s= 10.9337xxy21最終反應(yīng)速度為最終反應(yīng)速度為半速度點(diǎn)半速度點(diǎn)(達(dá)到最終速度一半達(dá)到最終速度一半時的時的x值值

28、)為為6831.21210641. 02其它輸出其它輸出命令命令nlintool 給出交互畫面給出交互畫面00.511.5050100150200250o 原始數(shù)據(jù)原始數(shù)據(jù)+ 擬合結(jié)果擬合結(jié)果 00.20.40.60.81-5005010015020025029/332021-11-5混合反應(yīng)混合反應(yīng)模型模型 x1為底物濃度，為底物濃度， x2為一示性變量為一示性變量 x2=1表示經(jīng)過處理，表示經(jīng)過處理，x2=0表示未經(jīng)處理表示未經(jīng)處理 1是未經(jīng)處理的最終反應(yīng)速度是未經(jīng)處理的最終反應(yīng)速度 1是經(jīng)處理后最終反應(yīng)速度的增長值是經(jīng)處理后最終反應(yīng)速度的增長值 2是未經(jīng)處理的反應(yīng)的半速度點(diǎn)是未經(jīng)處理的反

29、應(yīng)的半速度點(diǎn) 2是經(jīng)處理后反應(yīng)的半速度點(diǎn)的增長值是經(jīng)處理后反應(yīng)的半速度點(diǎn)的增長值 在同一模型中考慮嘌呤霉素處理的影響在同一模型中考慮嘌呤霉素處理的影響xxy2112221211)(xxxxy）（30/332021-11-5o 原始數(shù)據(jù)原始數(shù)據(jù)+ 擬合結(jié)果擬合結(jié)果 混合模型求解混合模型求解用用nlinfit 和和 nlintool命令命令,17001,6001,05. 00201. 002估計結(jié)果和預(yù)測估計結(jié)果和預(yù)測剩余標(biāo)準(zhǔn)差剩余標(biāo)準(zhǔn)差s= 10.4000 參數(shù)參數(shù)參數(shù)估計值參數(shù)估計值置信區(qū)間置信區(qū)間 1160.2802145.8466 174.7137 20.04770.0304 0.0650 152.403532.4130 72.3941 20.0164-0.0075 0.0403 2置信區(qū)間包含零點(diǎn)，置信區(qū)間包含零點(diǎn)，表明表明 2對因變量對因變量y的影響不顯著的影響不顯著12221211)(xxxxy）（參數(shù)初值參數(shù)初值(基于對數(shù)據(jù)的分析基于對數(shù)據(jù)的分析)經(jīng)嘌呤霉素處理的作用不影響半速度點(diǎn)參數(shù)經(jīng)嘌呤霉素處理的作用不影響半速度點(diǎn)參數(shù)未經(jīng)未