第十二章 回歸分析

1、第十二章 回歸分析前面我們講過曲線擬合問題。曲線擬合問題的特點是，根據(jù)得到的若干有關變量的一組數(shù)據(jù)，尋找因變量與（一個或幾個）自變量之間的一個函數(shù)，使這個函數(shù)對那組數(shù)據(jù)擬合得最好。通常，函數(shù)的形式可以由經(jīng)驗、先驗知識或?qū)?shù)據(jù)的直觀觀察決定，要作的工作是由數(shù)據(jù)用最小二乘法計算函數(shù)中的待定系數(shù)。從計算的角度看，問題似乎已經(jīng)完全解決了，還有進一步研究的必要嗎?從數(shù)理統(tǒng)計的觀點看，這里涉及的都是隨機變量，我們根據(jù)一個樣本計算出的那些系數(shù)，只是它們的一個（點）估計，應該對它們作區(qū)間估計或假設檢驗，如果置信區(qū)間太大，甚至包含了零點，那么系數(shù)的估計值是沒有多大意義的。另外也可以用方差分析方法對模型的誤差進行

2、分析，對擬合的優(yōu)劣給出評價。簡單地說，回歸分析就是對擬合問題作的統(tǒng)計分析。具體地說，回歸分析在一組數(shù)據(jù)的基礎上研究這樣幾個問題：（i）建立因變量與自變量之間的回歸模型（經(jīng)驗公式）；（ii）對回歸模型的可信度進行檢驗；（iii）判斷每個自變量對的影響是否顯著；（iv）診斷回歸模型是否適合這組數(shù)據(jù)；（v）利用回歸模型對進行預報或控制。§1 多元線性回歸回歸分析中最簡單的形式是，均為標量，為回歸系數(shù)，稱一元線性回歸。它的一個自然推廣是為多元變量，形如 (1)，或者更一般地 （2）其中，是已知函數(shù)。這里對回歸系數(shù)是線性的，稱為多元線性回歸。不難看出，對自變量作變量代換，就可將（2）化為（1）

3、的形式，所以下面以（1）為多元線性回歸的標準型。1.1 模型在回歸分析中自變量是影響因變量的主要因素，是人們能控制或能觀察的，而還受到隨機因素的干擾，可以合理地假設這種干擾服從零均值的正態(tài)分布，于是模型記作 （3）其中未知?，F(xiàn)得到個獨立觀測數(shù)據(jù)，由（3）得 （4）記， （5），（4）表為 （6）1.2 參數(shù)估計用最小二乘法估計模型（3）中的參數(shù)。由（4）式這組數(shù)據(jù)的誤差平方和為 （7）求使最小，得到的最小二乘估計，記作，可以推出 （8）將代回原模型得到的估計值 （9）而這組數(shù)據(jù)的擬合值為，擬合誤差稱為殘差，可作為隨機誤差的估計，而 （10）為殘差平方和（或剩余平方和），即。1.3 統(tǒng)計分析不加

4、證明地給出以下結(jié)果：（i）是的線性無偏最小方差估計。指的是是的線性函數(shù)；的期望等于；在的線性無偏估計中，的方差最小。（ii）服從正態(tài)分布 （11）（iii）對殘差平方和，且 （12）由此得到的無偏估計 （13）是剩余方差（殘差的方差），稱為剩余標準差。（iv）對的樣本方差進行分解，有 ， （14）其中是由（10）定義的殘差平方和，反映隨機誤差對的影響，稱為回歸平方和，反映自變量對的影響。1.4 回歸模型的假設檢驗因變量與自變量之間是否存在如模型（1）所示的線性關系是需要檢驗的，顯然，如果所有的 都很小，與的線性關系就不明顯，所以可令原假設為 當成立時由分解式（14）定義的滿足 (15)在顯著性

5、水平下有分位數(shù)，若，接受；否則，拒絕。注意 拒絕只說明與的線性關系不明顯，可能存在非線性關系，如平方關系。還有一些衡量與相關程度的指標，如用回歸平方和在樣本方差中的比值定義 （16）稱為相關系數(shù)，越大，與相關關系越密切，通常，大于0.8（或0.9）才認為相關關系成立。1.5 回歸系數(shù)的假設檢驗和區(qū)間估計當上面的被拒絕時，不全為零，但是不排除其中若干個等于零。所以應進一步作如下個檢驗： 由（11）式，是對角線上的元素，用代替，由（11）（13）式，當成立時 （17）對給定的，若，接受；否則，拒絕。(17)式也可用于對作區(qū)間估計（），在置信水平下，的置信區(qū)間為 （18）其中。1.6 利用回歸模型進

6、行預測當回歸模型和系數(shù)通過檢驗后，可由給定的預測，是隨機的，顯然其預測值（點估計）為 （19）給定可以算出的預測區(qū)間（區(qū)間估計），結(jié)果較復雜，但當較大且接近平均值時，的預測區(qū)間可簡化為 （20）其中是標準正態(tài)分布的分位數(shù)。對的區(qū)間估計方法可用于給出已知數(shù)據(jù)殘差的置信區(qū)間，服從均值為零的正態(tài)分布，所以若某個的置信區(qū)間不包含零點，則認為這個數(shù)據(jù)是異常的，可予以剔除。1.7 Matlab實現(xiàn)Matlab統(tǒng)計工具箱用命令regress實現(xiàn)多元線性回歸，用的方法是最小二乘法，用法是：b=regress(Y,X)其中Y,X為按（5）式排列的數(shù)據(jù)，b為回歸系數(shù)估計值。b,bint,r,rint,stats=

7、regress(Y,X,alpha)這里Y,X同上，alpha為顯著性水平（缺省時設定為0.05），b,bint為回歸系數(shù)估計值和它們的置信區(qū)間，r,rint為殘差（向量）及其置信區(qū)間，stats是用于檢驗回歸模型的統(tǒng)計量，有三個數(shù)值，第一個是（見（16）式），第二個是（見（15）式），第3個是與對應的概率，拒絕，回歸模型成立。殘差及其置信區(qū)間可以用rcoplot(r,rint)畫圖。例1 合金的強度與其中的碳含量有比較密切的關系，今從生產(chǎn)中收集了一批數(shù)據(jù)如下表： 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 42.0 41.5 45.0 45.5

8、45.0 47.5 49.0 55.0 50.0試先擬合一個函數(shù)，再用回歸分析對它進行檢驗。解 先畫出散點圖：x=0.1:0.01:0.18;y=42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0;plot(x,y,'+')可知與大致上為線性關系。設回歸模型為 （21）用regress和rcoplot編程如下：clc,clearx1=0.1:0.01:0.18'y=42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0'x=ones(9,1),x1;b,bint,r,rint,stats=regre

9、ss(y,x);b,bint,stats,rcoplot(r,rint)得到b =27.4722 137.5000bint =18.6851 36.2594 75.7755 199.2245stats =0.7985 27.7469 0.0012即，的置信區(qū)間是18.6851,36.2594，的置信區(qū)間是75.7755,199.2245；，?？芍Ｐ停?1）成立。觀察命令rcoplot(r,rint)所畫的殘差分布，除第8個數(shù)據(jù)外其余殘差的置信區(qū)間均包含零點，第8個點應視為異常點，將其剔除后重新計算，可得b =30.7820 109.3985bint =26.2805 35.2834 76.9

10、014 141.8955stats =0.9188 67.8534 0.0002應該用修改后的這個結(jié)果。例2 某廠生產(chǎn)的一種電器的銷售量與競爭對手的價格和本廠的價格有關。下表是該商品在10個城市的銷售記錄。元 120 140 190 130 155 175 125 145 180 150元 100 110 90 150 210 150 250 270 300 250個 102 100 120 77 46 93 26 69 65 85試根據(jù)這些數(shù)據(jù)建立與和的關系式，對得到的模型和系數(shù)進行檢驗。若某市本廠產(chǎn)品售價160（元），競爭對手售價170（元），預測商品在該市的銷售量。解 分別畫出關于和關于

11、的散點圖，可以看出與有較明顯的線性關系，而與之間的關系則難以確定，我們將作幾種嘗試，用統(tǒng)計分析決定優(yōu)劣。設回歸模型為 （22）編寫如下程序：x1=120 140 190 130 155 175 125 145 180 150'x2=100 110 90 150 210 150 250 270 300 250'y=102 100 120 77 46 93 26 69 65 85'x=ones(10,1),x1,x2;b,bint,r,rint,stats=regress(y,x);b,bint,stats得到b =66.5176 0.4139 -0.2698bint =-

12、32.5060 165.5411 -0.2018 1.0296 -0.4611 -0.0785stats =0.6527 6.5786 0.0247可以看出結(jié)果不是太好：，取時回歸模型（22）可用，但取則模型不能用；較??；的置信區(qū)間包含了零點。下面將試圖用的二次函數(shù)改進它。1.8 多項式回歸如果從數(shù)據(jù)的散點圖上發(fā)現(xiàn)與呈較明顯的二次（或高次）函數(shù)關系，或者用線性模型（1）的效果不太好，就可以選用多項式回歸。1.8.1 一元多項式回歸一元多項式回歸可用命令polyfit實現(xiàn)。例3 將17至29歲的運動員每兩歲一組分為7組，每組兩人測量其旋轉(zhuǎn)定向能力，以考察年齡對這種運動能力的影響?，F(xiàn)得到一組數(shù)據(jù)如

13、下表：年 齡 17 19 21 23 25 27 29第一人 20.48 25.13 26.15 30.0 26.1 20.3 19.35第二人 24.35 28.11 26.3 31.4 26.92 25.7 21.3試建立二者之間的關系。解 數(shù)據(jù)的散點圖明顯地呈現(xiàn)兩端低中間高的形狀，所以應擬合一條二次曲線。選用二次模型 （23）編寫如下程序：x0=17:2:29;x0=x0,x0;y0=20.48 25.13 26.15 30.0 26.1 20.3 19.35. 24.35 28.11 26.3 31.4 26.92 25.7 21.3;p,s=polyfit(x0,y0,2); p得到

14、p =-0.2003 8.9782 -72.2150即，。上面的s是一個數(shù)據(jù)結(jié)構(gòu)，用于計算其它函數(shù)的計算，如y,delta=polyconf(p,x0,s);y得到的擬合值，及預測值的置信區(qū)間半徑delta。用polytool(x0,y0,2)，可以得到一個如上圖的交互式畫面，在畫面中綠色曲線為擬合曲線，它兩側(cè)的紅線是的置信區(qū)間。你可以用鼠標移動圖中的十字線來改變圖下方的值，也可以在窗口內(nèi)輸入，左邊就給出的預測值及其置信區(qū)間。通過左下方的Export下拉式菜單，可以輸出回歸系數(shù)等。這個命令的用法與下面將介紹的rstool相似。1.8.2 多元二項式回歸統(tǒng)計工具箱提供了一個作多元二項式回歸的命令

15、rstool，它也產(chǎn)生一個交互式畫面，并輸出有關信息，用法是rstool(x,y,model,alpha)其中輸入數(shù)據(jù)x,y分別為矩陣和維向量，alpha為顯著性水平（缺省時設定為0.05），model由下列4個模型中選擇1個（用字符串輸入，缺省時設定為線性模型）： linear(線性)： purequadratic(純二次)： interaction（交叉）： quadratic(完全二次)：我們再作一遍例2 商品銷售量與價格問題，選擇純二次模型，即 （24）編程如下：x1=120 140 190 130 155 175 125 145 180 150'x2=100 110 90 1

16、50 210 150 250 270 300 250'y=102 100 120 77 46 93 26 69 65 85'x=x1 x2;rstool(x,y,'purequadratic') 得到一個如圖所示的交互式畫面，左邊是（=151）固定時的曲線及其置信區(qū)間，右邊是（=188）固定時的曲線及其置信區(qū)間。用鼠標移動圖中的十字線，或在圖下方窗口內(nèi)輸入，可改變。圖左邊給出的預測值及其置信區(qū)間，就用這種畫面可以回答例2提出的“若某市本廠產(chǎn)品售價160（元），競爭對手售價170（元），預測該市的銷售量”問題。圖的左下方有兩個下拉式菜單，一個菜單Export用以向

17、Matlab工作區(qū)傳送數(shù)據(jù)，包括beta(回歸系數(shù))，rmse（剩余標準差），residuals(殘差)。模型（24）的回歸系數(shù)和剩余標準差為beta =-312.5871 7.2701 -1.7337 -0.0228 0.0037rmse =16.6436另一個菜單model用以在上述4個模型中選擇，你可以比較以下它們的剩余標準差，會發(fā)現(xiàn)以模型（24）的rmse=16.6436最小。§2 非線性回歸和逐步回歸本節(jié)介紹怎樣用Matlab統(tǒng)計工具箱實現(xiàn)非線性回歸和逐步回歸。2.1 非線性回歸非線性回歸是指因變量對回歸系數(shù)（而不是自變量）是非線性的。Matlab統(tǒng)計工具箱中的nlinfi

18、t，nlparci，nlpredci，nlintool，不僅給出擬合的回歸系數(shù)，而且可以給出它的置信區(qū)間，及預測值和置信區(qū)間等。下面通過例題說明這些命令的用法。例4 在研究化學動力學反應過程中，建立了一個反應速度和反應物含量的數(shù)學模型，形式為 其中是未知的參數(shù)，是三種反應物（氫，戊烷，異構(gòu)戊烷）的含量，是反應速度。今測得一組數(shù)據(jù)如下表，試由此確定參數(shù)，并給出其置信區(qū)間。的參考值為（0.1，0.05，0.02，1，2）。序號 反應速度 氫 戊烷 異構(gòu)戊烷 1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120

19、 5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.32 285 300 10 13 3.13 285 190 120解 首先，以回歸系數(shù)和自變量為輸入變量，將要擬合的模型寫成函數(shù)文件huaxue.m：function yhat=huaxue(beta,x);yhat=(beta(4)*x(2)-x(3)/beta(5)./(1+beta(1)*x(1)+.beta(2)*x(

20、2)+beta(3)*x(3);然后，用nlinfit計算回歸系數(shù)，用nlparci計算回歸系數(shù)的置信區(qū)間，用nlpredci計算預測值及其置信區(qū)間，編程如下：clc,clearx0= 1 8.55 470 300 10 2 3.79 285 80 10 3 4.82 470 300 120 4 0.02 470 80 120 5 2.75 470 80 10 6 14.39 100 190 10 7 2.54 100 80 65 8 4.35 470 190 65 9 13.00 100 300 54 10 8.50 100 300 120 11 0.05 100 80 120 12 11.

21、32 285 300 10 13 3.13 285 190 120;x=x0(:,3:5);y=x0(:,2);beta=0.1,0.05,0.02,1,2; %回歸系數(shù)的初值betahat,f,j=nlinfit(x,y,'huaxue',beta); %f,j是下面命令用的信息betaci=nlparci(betahat,f,j);betaa=betahat,betaci %回歸系數(shù)及其置信區(qū)間yhat,delta=nlpredci('huaxue',x,betahat,f,j) %y的預測值及其置信區(qū)間的半徑，置信區(qū)間為yhat±delta。用n

22、lintool得到一個交互式畫面，左下方的Export可向工作區(qū)傳送數(shù)據(jù)，如剩余標準差等。使用命令nlintool(x,y,'huaxue',beta)可看到畫面，并傳出剩余標準差rmse= 0.1933。2.2 逐步回歸實際問題中影響因變量的因素可能很多，我們希望從中挑選出影響顯著的自變量來建立回歸模型，這就涉及到變量選擇的問題，逐步回歸是一種從眾多變量中有效地選擇重要變量的方法。以下只討論線性回歸模型（1）式的情況。變量選擇的標準，簡單地說就是所有對因變量影響顯著的變量都應選入模型，而影響不顯著的變量都不應選入模型，從便于應用的角度應使模型中變量個數(shù)盡可能少。若候選的自變量

23、集合為，從中選出一個子集，設中有個自變量，由和因變量構(gòu)造的回歸模型的誤差平方和為，則模型的剩余標準差的平方，為數(shù)據(jù)樣本容量。所選子集應使盡量小，通?；貧w模型中包含的自變量越多，誤差平方和越小，但若模型中包含有對影響很小的變量，那么不會由于包含這些變量在內(nèi)而減少多少， 卻因的增加可能使反而增大，同時這些對影響不顯著的變量也會影響模型的穩(wěn)定性，因此可將剩余標準差最小作為衡量變量選擇的一個數(shù)量標準。逐步回歸是實現(xiàn)變量選擇的一種方法，基本思路為，先確定一初始子集，然后每次從子集外影響顯著的變量中引入一個對影響最大的，再對原來子集中的變量進行檢驗，從變得不顯著的變量中剔除一個影響最小的，直到不能引入和剔

24、除為止。使用逐步回歸有兩點值得注意，一是要適當?shù)剡x定引入變量的顯著性水平和剔除變量的顯著性水平，顯然，越大，引入的變量越多；越大，剔除的變量越少。二是由于各個變量之間的相關性，一個新的變量引入后，會使原來認為顯著的某個變量變得不顯著，從而被剔除，所以在最初選擇變量時應盡量選擇相互獨立性強的那些。在Matlab統(tǒng)計工具箱中用作逐步回歸的是命令stepwise，它提供了一個交互式畫面，通過這個工具你可以自由地選擇變量，進行統(tǒng)計分析，其通常用法是：stepwise(x,y,inmodel,alpha)其中x是自變量數(shù)據(jù)，y是因變量數(shù)據(jù)，分別為和矩陣，inmodel是矩陣x的列數(shù)的指標，給出初始模型中

25、包括的子集（缺省時設定為全部自變量），alpha為顯著性水平。stepwise命令產(chǎn)生三個圖形窗口：Stepwise Table，Stepwise History，Stepwise Plot。Stepwise Table窗口中列出了一個統(tǒng)計表，包括回歸系數(shù)及其置信區(qū)間，模型的統(tǒng)計量（RMSE R-square,F,p等，其含義與regress,rstool相同）。你可以通過這些統(tǒng)計量的變化來確定模型。Stepwise History窗口顯示RMSE的值及其置信區(qū)間。Stepwise Plot窗口，顯示回歸系數(shù)及其置信區(qū)間，綠色表明在模型中的變量，紅色表明從模型中移去的變量，兩邊有虛線或?qū)嵕€，虛

26、線表示該變量的擬合系數(shù)與零無顯著差異，實線則表明有顯著差異。在這個窗口中還有Scale Inputs和Export按鈕。按下Scale Inputs表明對于輸入數(shù)據(jù)的每列進行正態(tài)化處理，使其標準差為1。點擊Export產(chǎn)生一個菜單，表明了要傳送給Matlab工作區(qū)的參數(shù)，它們給出了統(tǒng)計計算的一些結(jié)果。下面通過一個例子說明stepwise的用法。例5 水泥凝固時放出的熱量與水泥中4種化學成分有關，今測得一組數(shù)據(jù)如下，試用逐步回歸來確定一個線性模型序號 1 7 26 6 60 78.5 2 1 29 15 52 74.3 3 11 56 8 20 104.3 4 11 31 8 47 87.6 5

27、 7 52 6 33 95.9 6 11 55 9 22 109.2 7 3 71 17 6 102.7 8 1 31 22 44 72.5 9 2 54 18 22 93.1 10 21 47 4 26 115.9 11 1 40 23 34 83.8 12 11 66 9 12 113.3 13 10 68 8 12 109.4編寫程序如下：clc,clearx0=1 7 26 6 60 78.5 2 1 29 15 52 74.3 3 11 56 8 20 104.3 4 11 31 8 47 87.6 5 7 52 6 33 95.9 6 11 55 9 22 109.2 7 3 71 17 6 102.7 8 1 31 22 44 72.5 9 2 54 18 22 93.1 10 21 47 4 26 115.9 11 1 40 23 34 83.8 12 11 66 9 12 113.3 13 10 68 8 12 109.4;x=x0(:,2:5);y=x0(:,6);stepwise(x,y)得到Stepwise Table如下：可以看出，不顯著