第八章--統(tǒng)計回歸模型

1、第八章統(tǒng)計回歸模型回歸分析是研究一個變量 Y與其它假設(shè)干變量 X之間相關(guān)關(guān)系的一種數(shù)學(xué)工具 它是在一組試驗 或觀測數(shù)據(jù)的根底上，尋找被隨機性掩蓋了的變量之間的依存關(guān)系粗略的講，可以理解為用一種確定的函數(shù)關(guān)系去近似代替比擬復(fù)雜的相關(guān)關(guān)系這個函數(shù)稱為回歸函數(shù)回歸分析所研究的主要問題是如何利用變量X、Y的觀察值(樣本)，對回歸函數(shù)進行統(tǒng)計推斷，包括對它進行估計及檢驗與它有關(guān)的假設(shè)等回歸分析包含的內(nèi)容廣泛此處將討論多項式回歸、多元線性回歸、非線性回歸以及逐步回歸、多項式回歸(1) 一元多項式回歸元多項式回歸模型的一般形式為y 0 ixmmX如果從數(shù)據(jù)的散點圖上發(fā)現(xiàn)回歸y與x呈現(xiàn)較明顯的二次(或高次)函

2、數(shù)關(guān)系，那么可以選用一元多項式1. 用函數(shù)polyfit估計模型參數(shù)，其具體調(diào)用格式如下:p=polyfit(x,y,m) p返回多項式系數(shù)的估計值； m設(shè)定多項式的最高次數(shù)； x，y為對應(yīng)數(shù)據(jù)點值.p,S=polyfit(x,y,m) S是一個矩陣，用來估計預(yù)測誤差2. 輸出預(yù)估值與殘差的計算用函數(shù)polyval實現(xiàn)，其具體調(diào)用格式如下：Y=polyval(p,X) 求polyfit所得的回歸多項式在 X處的預(yù)測值 Y.Y,DELTA=polyval(p,X,S) p，S為polyfit的輸出，DELTA為誤差估計.在線性回歸模型中，丫土DELTA以50%的概率包含函數(shù)在 X處的真值.3.

3、模型預(yù)測的置信區(qū)間用polyconf實現(xiàn)，其具體調(diào)用格式如下：Y,DELTA=polyconf(p,X,S,alpha) 求polyfit所得的回歸多項式在 X處的預(yù)測值 Y及預(yù)測值的顯著 性為1-alpha的置信區(qū)間 丫土DELTA，alpha缺省時為0.05.4. 交互式畫圖工具polytool，其具體調(diào)用格式如下：polytool(x,y,m)；polytool(x,y,m,alpha)；用m次多項式擬合x，y的值，默認值為1，alpha為顯著性水平.例1觀測物體降落的距離s與時間t的關(guān)系，得到數(shù)據(jù)如下表，求s.t (s)1/302/303/304/305/306/307/30s (cm

4、)t (s)8/309/3010/3011/3012/3013/3014/30s (cm)解根據(jù)數(shù)據(jù)的散點圖，應(yīng)擬合為一條二次曲線.選用二次模型，具體代碼如下:%輸入數(shù)據(jù)t=1/30:1/30:14/30;s=11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 113.77 129.54 146.48;%多項式系數(shù)擬合p,S=polyfit(t,s,2);那么得回歸模型為：2? 489.2946t65.8896t9.1329.%的擬合值及預(yù)測值 y的置信半徑deltay,dalta=polyc on f(p,t,S

5、);得結(jié)果如下：y=Colu mns 1 through 11Colu mns 12 through 14113.7759 129.5637dalta=Colu mns 1 through 110.0937 0.0865 0.08290.0817 0.08230.0827Colu mns 12 through 140.0829 0.0865%交互式畫圖polytool(t,s,2);polytool所得的交互式圖形如圖8-1所示.圖8-1(2)多元二項式回歸多元二項式回歸模型的一般形式為y 01x1. mXmjkXjXk多元二項式回歸命令：rstoolx,y,model ,alplxa表示 n

6、 m 矩陣；y 表示 n維列向量； alpha 為顯表示由以下4個模型中選擇model1個用字符串輸入，缺省時為線性模型:著性水平缺省時為0.05；1X1purequadratic純二次:1X12jjXj ；in teract ion交叉: y1X1m Xmjkxj xk ；quadratic完全二次：y1X1jk Xj Xk .linear線性：y 。需求量10075807050659010011060收入10006001200500300400130011001300300價格5766875439例2設(shè)某商品的需求量與消費者的平均收入、 均收入為1000價格為6時的商品需求量.j,k m商

7、品價格的統(tǒng)計數(shù)據(jù)如下，建立回歸模型，預(yù)測平解 選擇純二次模型，即y 01x12X2211 X1222 X2 .%輸入數(shù)據(jù)x1=1000 600 1200 500 300 400 1300 1100 1300 300; x2=5 7 6 6 8 7 5 4 3 9;x=x1 x2;y=100 75 80 70 50 65 90 100 110 60;%多元二項式回歸rstoolx,y,purequadratic;得如下結(jié)果：圖8-2得到一個如下列圖的交互式畫面，左邊是 x1=1000固定時的曲線y x1及其置信區(qū)間，右邊是x2 =6固定時的曲線 yx2及其置信區(qū)間.用鼠標移動圖中的十字線，或在圖

8、下方窗口內(nèi)輸入，可改變x1,x2.在左邊圖形下方的方框中輸入1000，右邊圖形下方的方框中輸入6,那么畫面左邊的“PredictedY1 下方的數(shù)據(jù)變?yōu)?,.在畫面左下方單擊Export在出現(xiàn)的窗體中單擊 ok按鈕，那么beta、rmse和residuals都傳送到Matlab工作區(qū)中.在Matlab工作區(qū)中輸入命令：beta,rmse,得結(jié)果：beta故回歸模型為：y 110.5313 0.1464為 26.5709x2 0.0001x2 1.8475x2,，說明此回歸模型的顯著性較好.二、多元線性回歸多元線性回歸模型的一般形式為y 01x1Xmm在Matlab統(tǒng)計工具箱中使用函數(shù)regre

9、ss實現(xiàn)多元線性回歸.具體調(diào)用格式為:b=regress(Y,X)b,b in t,r,ri nt,stats=regress(Y,X,alpha)Y1X11X12X1 m其中 Y丫2 , X1X21X22x2m.對于 兀線性回歸，取m 1即可.b為輸出冋量;Yn1Xn1xn2xnmb, bi nt表示回歸系數(shù)估計值和它們的置信區(qū)間;檢驗回歸模型的統(tǒng)計量，有四個數(shù)值：相關(guān)系數(shù)r表示殘差；rint表示殘差的置信區(qū)間；stats表示用于2 2R、F值、與F值對應(yīng)的概率P、s的值.相關(guān)系數(shù)R2越接近1,說明回歸方程越顯著;FF1 m, n m 1時拒絕H0, F越大，說明回歸方程越顯著；與F對應(yīng)的概

10、率P 時拒絕H0，回歸模型成立；alpha表示顯著性水平缺省時為0.05.殘差及其置信區(qū)間可以用命令rcoplotr,rint畫出.例3某湖泊八年來湖水中COD濃度實測值y與影響因素，如湖區(qū)工業(yè)產(chǎn)值X、總?cè)丝跀?shù)X2、捕魚量X3、降水量X4的資料，建立y的水質(zhì)分析模型.湖水濃度與影響因素數(shù)據(jù)表X1X2X3X4y解 作出因變量y與各自變量的樣本散點圖作散點圖的目的主要是觀察因變量y與各自變量間是否有比擬好的線性關(guān)系，以便選擇恰當?shù)臄?shù)學(xué)模型形式圖8-3、圖8-4、圖8-5、圖8-6分別為y與Xi、x2、x3、x4的散點圖.從圖中可以看出這些點大致分布在一條直線旁邊，因此有較好的線性關(guān)系，可以采用線性

11、回歸圖8-4 y與X2的散點圖圖8-3 y與xi的散點圖圖8-5 y與X3的散點圖在Matlab中實現(xiàn)回歸的具體代碼如下:%輸入數(shù)據(jù)圖8-6 y與X4的散點圖x1=1.376 1.375 1.387 1.401 1.412 1.428 1.445 1.477;x2=0.450 0.475 0.485 0.500 0.535 0.545 0.550 0.575;x3=2.170 2.554 2.676 2.713 2.823 3.088 3.122 3.262; x4=0.8922 1.1610 0.5346 0.9589 1.0239 1.0499 1.1065 1.1387; x=o nes

12、(8,1) x1 x2 x3 x4;y=5.19 5.30 5.60 5.82 6.00 6.06 6.45 6.95;%多元線性回歸b,bi nt,r,ri nt,stats=regress(y,x);得如下結(jié)果：bint =rint =stats =故回歸模型為：y 13.9849 13.1920X, 2.4228x2 0.0754x3 0.1897x4 ,此外，由 stats 的值可知 R20.9846 , F 47.9654 , P 0.0047。%殘差分析，作殘差圖rcoplot(r,ri nt)圖8-7從殘差圖可以看出， 除第一和第六個數(shù)據(jù)外， 其余數(shù)據(jù)的殘差離零點均較近， 且殘差

13、的置信區(qū)間均 包含零點 . 第一和第六個數(shù)據(jù)可視為異常點，將其去掉后重新擬合可得新的回歸模型 .三、非線性回歸非線性回歸模型的一般形式為y f(X, ), x(Xi，X2，, Xm) ,( 1, 2 , k)，其中f 對回歸系數(shù) 是非線性的 . 非線性回歸包括如下幾個用于計算回歸參數(shù)、預(yù)測輸出、置信區(qū)間以及輸出交互圖像的函數(shù).1. 非線性最小二乘參數(shù)估計對于非線性方程的的系數(shù)估計通常采用最小二乘估計，又叫做非線性最小二乘回歸在Matlab中采用 nlinfit 實現(xiàn)，其調(diào)用格式如下： beta=nlinfit(X,y,fun,beta0) 返回非線性回歸方程系數(shù)的最小二乘估計值. 非線性方程由

14、 fun 給定， fun 為用戶提供形如y f (beta,x)的函數(shù)，beta為待估系數(shù)；betaO為回歸系數(shù) 的初值.beta,r,J=nlinfit(X,y,fun,beta0)返回回歸系數(shù) beta、殘差r、Jacobi矩陣J.2. 最小二乘估計參數(shù)的置信區(qū)間求非線性最小二乘估計的系數(shù)的置信區(qū)間用 nlparci 計算。其輸入為 nlinfit 函數(shù)的輸出 beta， r， J。函數(shù)具體調(diào)用格式如下：ci=nlparci(beta,r,J) 返回系數(shù) beta 的 95%置信區(qū)間 .ci=nlparci(beta,r,J,alpha) 返回系數(shù)beta的100(1-alpha)%置信區(qū)

15、間.3. 最小二乘估計模型的預(yù)測輸出及其置信區(qū)間 非線性最小二乘估計模型的預(yù)測輸出及其置信區(qū)間用 nlpredci 計算 . 函數(shù)具體調(diào)用格式如下：ypred,delta=nlpredci(fun,inputs,beta,r,J)返回模型在對應(yīng)in puts處的輸出預(yù)測值 ypred，給出95%的置信區(qū)間ypred-delta, ypred+delta； ypred,delta=nlpredci(fun,inputs,beta,r,J,alpha)返回模型在對應(yīng) inputs 處的輸出預(yù)測值 ypred ，給出 100(1-alpha)% 的置信區(qū)間 ypred-delta, ypred+de

16、lta .4. 非線性擬合和預(yù)測的交互圖形工具nlintool 是非線性擬合和預(yù)測的交互圖形工具 . 函數(shù)具體調(diào)用格式如下： nlintool(X,y,fun,beta0)返回 X， y 的非線性最小二乘法的曲線擬合圖，并畫出 95%的置信區(qū)間；nlintool(x,y,fun,beta0,alpha)給出曲線擬合圖及 100(1-alpha)% 的置信區(qū)間；nlintool(x,y,fun, beta0,alpha, xnamey,name)給出曲線擬合圖及 100(1-alpha)%的置信區(qū)間，標出 x, y變量名稱.例4在化工生產(chǎn)中獲得的氯氣的級分y隨生產(chǎn)時間x下降，假定在x 8時，y與

17、x之間有如下形式的非線性模型：a (0.49 a)eb(x 8)現(xiàn)收集 43 組數(shù)據(jù)， 如下所示 . 要求利用該數(shù)據(jù)求 a， b 的值以及它們 95%的置信區(qū)間 . 并畫出擬合曲線 此處設(shè) a， b .x:8.00,8.00,10.00,10.00,10.00,10.00,12.00,12.00,12.00,14.00,14.00,14.00,16.00,16.00,16.00,18.00,18.00, 20.00,20.00,20.00,20.00,22.00,22.00,24.00,24.00,24.00,26.00,26.00,26.00,28.00,28.00,30.00,30.00,

18、30.00,y:0.49,0.49,0.48,0.47,0.48,0.47,0.46,0.46,0.45,0.43,0.45,0.43,0.43,0.44,0.43,0.43,0.46,0.42,0.42,0.43,0.41,解首先編寫函數(shù)fun，表示出需要擬合模型的函數(shù)形式 %擬%合模型的函數(shù)形式function y=fun(beta0,x)a=beta0(1);b=beta0(2);y=a+(0.49-a)*exp(-b*(x-8);編寫程序?qū)?shù)據(jù)進行最小二乘擬合和處理%輸%入數(shù)據(jù)x=8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 1

19、4.00 14.00 14.00 16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00 22.00 22.00 24.00 24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00 34.00 36.00 36.00 38.00 38.00 40.00 42.00;y=0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43 0.43 0.46 0.42 0.42

20、0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41 0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39;beta0=0.30 0.02;beta,r,J=nlinfit(x,y,fun,beta0);beta%系%數(shù)初值%系%數(shù)最小二乘估計值為betaci=nlparci(beta,r,J);nlintool(x,y,fun,betaO,生產(chǎn)時間 x,氯氣級分 y);%求% 95%置信區(qū)間%顯%示擬合結(jié)果和置信區(qū)間得如下結(jié)果：beta=ci=可知a與b的最小二

21、乘估計值分別為和，它們95%的置信區(qū)間分別為0.3805,0.3988和 0.0771,0.1251.所以用最小二乘法對該模型系數(shù)進行估計，最終的模型表達式為0.3896 0.1004e0.1011 (x 8)最終得到的擬合圖如下所示圖8-8在圖8-8中，中間的實線為擬合曲線，兩側(cè)的虛線標出其95%置信區(qū)間.四、逐步回歸實際問題中影響因變量的因素可能很多，希望從中挑選出影響顯著的自變量來建立回歸模型.這就涉及到了變量選擇的問題.逐步回歸就是一種從眾多變量中有效地選擇重要變量的方法逐步回歸的命令是 stepwise，它提供了一個交互式畫面，通過此工具可以自由地選擇變量，進行統(tǒng)計分析，其調(diào)用格式為

22、stepwisex,y,inmodel,alpha.其中x表示自變量數(shù)據(jù)， n m階矩陣；y表示因變量數(shù)據(jù)，n 1階矩陣；inmodel表示矩陣的列 數(shù)指標缺省時設(shè)定為全部自變量；alpha表示顯著性水平缺省時為0.5.運行 stepwise命令時產(chǎn)生三個圖形窗口：Stepwise Plot， Stepwise Table， Stepwise History.在Stepwise Plot窗口，虛線表示該變量的擬合系數(shù)與0無顯著差異，實線表示有顯著差異，紅色線表示從模型中移去的變量，藍色線表示存在模型中的變量，單擊一條會改變其狀態(tài)Stepwise Table窗口中列出了一個統(tǒng)計表，包括回歸系數(shù)及

23、其置信區(qū)間，以及模型的統(tǒng)計量剩余標準差RMSE、相關(guān)系數(shù)R-square、F值、與F對應(yīng)的概率P .例5水泥凝固時放出的熱量 y與水泥中4種化學(xué)成分x1、x2、x3、x4有關(guān)，今測得一組數(shù)據(jù)如下, 試用逐步回歸法確定一個線性模型序號12345678910111213x17111117113122111110x226295631525571315447406668x3615886917221842398x46052204733226442226341212y解 在Matlab中實現(xiàn)逐步回歸的具體代碼如下：%數(shù)據(jù)輸入：x1=7 1 11 11 7 11 3 1 2 21 1 11 10;x2=26

24、 29 56 31 52 55 71 31 54 47 40 66 68;x3=6 15 8 8 6 9 17 22 18 4 23 9 8;x4=60 52 20 47 33 22 6 44 22 26 34 12 12;x=x1 x2 x3 x4;y=78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4;%逐步回歸%先在初始模型中取全部自變量：stepwise(x,y)得圖8-9所示的結(jié)果.圖8-9從圖8-9中可以看出，x3、x4的P值較大，其對模型影響不顯著假設(shè)要移去x3、x4兩個變量，在 Matla

25、b中輸入：stepwise(x,y,1,2)得圖8-10所示的結(jié)果.圖 8-10移去變量x3和x4后模型具有顯著性，雖然剩余標準差RMSE沒有太大的變化，但是統(tǒng)計量F的值明顯增大，因此新的回歸模型更好對變量y和x1、x2作線性回歸，在 Matlab中輸入：X=o nes13,1 x1 x2;b=regressy,X得結(jié)果：b=討論題1.根據(jù)經(jīng)驗，在人的身高相等的情況下，血壓的收縮壓y與體重x1,年齡x2有關(guān),現(xiàn)在收集了 13個男子的有關(guān)數(shù)據(jù)，如下表所示，試建立y與x1、x2的回歸模型x1x250202030305060504055404020y1201411241261171251231251

26、321231321551472.下表為1980-1991年間以1987年不變價計算的美國個人消費支出y與美國國內(nèi)生產(chǎn)總值 x的數(shù)據(jù)單位：10億美元.1在直角坐標系下，作x與y的散點圖，并判斷y與x是否存在線性相關(guān)關(guān)系？2試求y與x的一元線性回歸方程.3對所得的回歸方程作顯著性檢驗0.05.4假設(shè)x0=4500，試求對應(yīng)的y0的點預(yù)測和包含概率為95%的區(qū)間預(yù)測.年份yx年份yx1980198619811987198219881983198919841990198519913.某地人事部門研究中學(xué)教師的薪金與他們的資歷、性別、教育程度、及培訓(xùn)情況等因素之間的關(guān)系，要建立一個數(shù)學(xué)模型，分析人事策略

27、的合理性，特別是考察女教師是否受到不公正的待遇，以及她們的婚姻狀況是否會影響收入 為此，從當?shù)亟處熤须S機選中 3414位進行觀察，然后從中保存了 90 個觀察對象，得到關(guān)鍵數(shù)據(jù)如下表 Zx1x2x3x4x5x6x79987000000101514110000102818110100125019110000102819010100102819000000101827000001107230000000129030110000120430010000135231012010120431000100110438000000111841110000112742000000125942110100112

28、742110000112742000100109547000001111352000001146252012010118254110000140454000100118254000000159455112110145966000100123767110100Zx1x2x3x4x5x6x7123767010100149675010000142478110100142479010000134791110100134392000001131094000100181410300211015341030000001430103110000143911111010019461141131102216114

29、114110183411411411114161170000012052139110100208714000211122641540021112201158114011299215911511116951620100001792167110100169017300000118271740000012604175112110172019901000017202090000002159209014100185221001000021042131101001852220000001185222200000022102221100002266223010000202722311000018522270

30、00100185223200000119952350000012616245113110232425311010018522570100012054260000000Zx1x2x3x4x5x6x72617284113110194828711000017202900100012604308112110185239110101194231900010020273251100001942326110100172032911010020483370000002334346112111172035500000119423571100002117380110001274238711211127404031

31、1211119424061101002266437010000243645301000020674580100002000464112110表中符號介紹如下：Z-月薪元;x1-工作時間月;x2=1-男性，x2=0-女性；x3=1-男性或單身女性，x3=0-已婚女性； x4-學(xué)歷取值0-6，值越大表示學(xué)歷越高；x5=1-受雇于重點中學(xué)，x5=0-其它；x6=1-受過培訓(xùn)的畢業(yè)生， x6=0-未受過培訓(xùn)的畢業(yè)生或受過培訓(xùn)的肄業(yè)生；x7=1-已兩年以上未從事教學(xué)工作，x7=0-其它。注意組合x2,x3=1,1，0,1，0,0的含義.進行變量選擇，建立變量x1-x7與Z的回歸模型不一定包括每個自變量，

32、說明教師的薪金與哪些變量的關(guān)系密切，是否存在性別和婚姻狀況上的差異為了數(shù)據(jù)處理上的方便，建議對薪金取對數(shù)后作為因變量4.熱敏電阻器的電阻y與溫度x之間的關(guān)系模型為2yiexp -x 3試根據(jù)下表中的數(shù)據(jù)給出參數(shù)1， 2，3的估計xyxyxyxy50347807016370908266110442755286107513720957030115382060236508011540100600512033076519630859744105514712528725.財政收入預(yù)測問題：財政收入與國民收入、工業(yè)總產(chǎn)值、農(nóng)業(yè)總產(chǎn)值、總?cè)丝?、就業(yè)人口、固 定資產(chǎn)投資等因素有關(guān)下表列出了 1952-1981年的原始數(shù)據(jù)，試構(gòu)造預(yù)測模型.年份國民收入 億元工業(yè)總產(chǎn) 值億元農(nóng)業(yè)總產(chǎn)值 億元總?cè)丝?萬人就業(yè)人口 萬人固定資產(chǎn)投資 億元財政收入億元1952598349461574822072944