課件-5-Matlab軟件包與多元回歸(共17頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上Matlab軟件包與多元回歸 回歸分析的方法以及“回歸（Regress）”這個(gè)名稱的起源，統(tǒng)計(jì)學(xué)史上一般歸功于英國(guó)生物學(xué)家、統(tǒng)計(jì)學(xué)家Francis Galton（著名的Darwin的表哥，cousin of ）。 Francis Galton（18221911） K. Pearson（18571936） 問(wèn)題：考察人體的某項(xiàng)指標(biāo)，比如說(shuō)：身高。假設(shè)父輩的身高為，孩子的身高為，身高具有遺傳性，但是，父輩的身高又不能完全決定孩子的身高。局部說(shuō)來(lái)，父輩是高個(gè)子，子女也是高個(gè)子，父輩是矮個(gè)子，子女也是矮個(gè)子。總體觀察，這個(gè)結(jié)論不對(duì)（否則，長(zhǎng)期進(jìn)化之后，人類的身高應(yīng)該兩極分化：

2、很高的人與很矮的人）。F Galton的學(xué)生K Pearson觀察了1078對(duì)夫婦，以每對(duì)夫婦的平均身高為，取他們的一個(gè)成年兒子的身高為，全體父母身高的平均值吋。K Pearson觀察得到：父母的身高中等，即，（吋）時(shí)，（吋），看起來(lái)遺傳起作用；父母的身高較高時(shí)，（吋）時(shí)，（吋），子女的身高變小了；父母的身高較矮時(shí)，（吋）時(shí)，（吋），子女的身高變大了。觀察說(shuō)明，父母的身高到了高和矮的邊界時(shí)，總有某種力量將子女的身高拉向中心平均值，這就叫回歸（Regress）。正是由于這種現(xiàn)象，Galton用“Regress”一詞來(lái)描述與的關(guān)系，盡管，這種“向中心回歸”的現(xiàn)象只在特殊的領(lǐng)域觀察得到，不具有普遍性

3、，但是，統(tǒng)計(jì)學(xué)界使用習(xí)慣了，所以，沿用至今?；貧w（Regress）的這種性質(zhì)，可以得到理論上的說(shuō)明。設(shè)，的相關(guān)系數(shù)，則可以證明，線性回歸方程中的系數(shù)滿足，當(dāng)時(shí)，這一點(diǎn)解釋了Galton指出的“Regress”現(xiàn)象：父輩身高的方差，子女身高的方差，在一代之間變化不大，故可以假設(shè)。于是，令，從線性回歸方程得到：，以及，這樣，就可以得到對(duì)于具體的，有，從可以看出：當(dāng)時(shí)，盡管，也增加，但是，增加一個(gè)單位，并不相應(yīng)增加一個(gè)單位，相當(dāng)于打了折，這就解釋了“向中心回歸”的現(xiàn)象。注：在許多實(shí)際問(wèn)題中有，此時(shí)，“向中心回歸”的現(xiàn)象就不存在了。（一）一般多元回歸一般在生產(chǎn)實(shí)踐和科學(xué)研究中，人們得到了參數(shù)和因變量的

4、數(shù)據(jù)，需要求出關(guān)系式，這時(shí)就可以用到回歸分析的方法。如果只考慮是線性函數(shù)的情形，當(dāng)自變量只有一個(gè)時(shí)，即，中時(shí)，稱為一元線性回歸，當(dāng)自變量有多個(gè)時(shí)，即，中時(shí)，稱為多元線性回歸。進(jìn)行線性回歸時(shí)，有4個(gè)基本假定： 待定參數(shù)（系數(shù)）是線性關(guān)系； 殘差是獨(dú)立的； 殘差滿足正態(tài)分布。 殘差滿足方差奇性（所謂方差齊性指的就是我們要比較的幾組數(shù)據(jù)是獨(dú)立的、且服從同方差的正態(tài)分布）；在Matlab軟件包中有一個(gè)做一般多元回歸分析的命令regress，調(diào)用格式如下：b, bint, r, rint, stats = regress(y,X,alpha) 或者b, bint, r, rint, stats = re

5、gress(y,X) 此時(shí)，默認(rèn)置信度alpha 0.05。這里，y是一個(gè)的列向量，X是一個(gè)的矩陣，其中第一列是全1向量（這一點(diǎn)對(duì)于回歸來(lái)說(shuō)很重要，這一個(gè)全1列向量對(duì)應(yīng)回歸方程的常數(shù)項(xiàng)），一般情況下，需要人工造一個(gè)全1列向量。回歸方程具有如下形式（待定參數(shù)具有線性關(guān)系）：其中，是殘差。在返回項(xiàng)b,bint,r,rint,stats中，是回歸方程的系數(shù)；是一個(gè)矩陣，它的第行表示的(1-alpha)可信區(qū)間；是的殘差列向量；是矩陣，它的第行表示第個(gè)殘差的(1-alpha)可信區(qū)間；注釋：殘差與殘差區(qū)間杠桿圖，最好在0點(diǎn)線附近比較均勻的分布，而不呈現(xiàn)一定的規(guī)律性，如果是這樣，就說(shuō)明回歸分析做得比較理

6、想。 一般的，返回4個(gè)值：值、F_檢驗(yàn)值、閾值，與顯著性概率相關(guān)的值（如果這個(gè)值不存在，則，只輸出前3項(xiàng)）。注釋：（1）一般說(shuō)來(lái)，值越大越好。（2）人們一般用以下統(tǒng)計(jì)量對(duì)回歸方程做顯著性檢驗(yàn)：F_檢驗(yàn)、t_檢驗(yàn)、以及相關(guān)系數(shù)檢驗(yàn)法。Matlab軟件包輸出F_檢驗(yàn)值和閾值。一般說(shuō)來(lái)，F(xiàn)_檢驗(yàn)值越大越好，特別的，應(yīng)該有F_檢驗(yàn)值。 我國(guó)著名統(tǒng)計(jì)學(xué)家許寶祿（19101970）教授證明：F檢驗(yàn)有多方面的優(yōu)良性。（3）與顯著性概率相關(guān)的值應(yīng)該滿足。如果，則說(shuō)明回歸方程中有多余的自變量，可以將這些多余的自變量從回歸方程中剔除（見下面逐步回歸的內(nèi)容）。這幾個(gè)技術(shù)指標(biāo)說(shuō)明擬合程度的好壞。這幾個(gè)指標(biāo)都好，就說(shuō)明

7、回歸方程是有意義的。例1（Hamilton，1987）數(shù)據(jù)如下：序號(hào)YX1X2112.372.239.66212.662.578.94312.003.874.40411.933.106.64511.063.394.91613.032.838.52713.133.028.04811.442.149.05912.863.047.711010.843.265.111111.203.395.051211.562.358.511310.832.766.591412.633.904.901512.463.166.96第一步 分析數(shù)據(jù)在Matlab軟件包中分析是否具有線性關(guān)系，并作圖觀察，M文件opt_ha

8、nmilton_1987：x1=2.23,2.57,3.87,3.10,3.39,2.83,3.02,2.14,3.04,3.26,3.39,2.35,2.76,3.90,3.16;x2=9.66,8.94,4.40,6.64,4.91,8.52,8.04,9.05,7.71,5.11,5.05,8.51,6.59,4.90,6.96;y=12.37,12.66,12.00,11.93,11.06,13.03,13.13,11.44,12.86,10.84,11.20,11.56,10.83,12.63,12.46;corrcoef(x1,y)corrcoef(x2,y)plot3(x1,x2

9、,y,*)得到結(jié)果：ans = 1.0000 0.0025 0.0025 1.0000ans = 1.0000 0.4341 0.4341 1.0000即，corrcoef(x1,y)0.0025，corrcoef(x2,y)0.4341，說(shuō)明沒(méi)有非常明顯的單變量線性關(guān)系。圖形如下：也看不出有線性關(guān)系，但是，旋轉(zhuǎn)圖形，可以看出所有點(diǎn)幾乎在一個(gè)平面上。這說(shuō)明，在一個(gè)平面上，滿足線性關(guān)系：或者，換成一個(gè)常見的形式 其中，是殘差。于是，在Matlab軟件包中做線性多元回歸，寫一個(gè)M文件opt_regress_hamilton：x1=2.23,2.57,3.87,3.10,3.39,2.83,3.02

10、,2.14,3.04,3.26,3.39,2.35,2.76,3.90,3.16;x2=9.66,8.94,4.40,6.64,4.91,8.52,8.04,9.05,7.71,5.11,5.05,8.51,6.59,4.90,6.96;y=12.37,12.66,12.00,11.93,11.06,13.03,13.13,11.44,12.86,10.84,11.20,11.56,10.83,12.63,12.46;e=ones(15,1);x=e,x1,x2;b,bint,r,rint,stats=regress(y,x,0.05)rcoplot(r,rint)其中，rcoplot（Res

11、idual case order plot）表示畫出殘差與殘差區(qū)間的杠桿圖。執(zhí)行后得到：b = -4.5154 3.0970 1.0319bint = -4.6486 -4.3822 3.0703 3.1238 1.0238 1.0399r = 0.0113 -0.0087 -0.0102 -0.0069 0.0101 -0.0106 -0.0037 -0.0105 0.0049 -0.0136 0.0057 0.0163 -0.0023 0.0110 0.0071rint = -0.0087 0.0314 -0.0303 0.0128 -0.0301 0.0098 -0.0299 0.016

12、2 -0.0106 0.0308 -0.0313 0.0102 -0.0252 0.0178 -0.0299 0.0089 -0.0174 0.0272 -0.0331 0.0058 -0.0161 0.0275 -0.0027 0.0354 -0.0236 0.0190 -0.0079 0.0299 -0.0156 0.0298stats = 1.0e+004 * 0.0001 3.9222 0 0.0000即，?？尚哦?5，且，與顯著性概率相關(guān)的，這說(shuō)明，回歸方程中的每個(gè)自變量的選取，都是有意義的。殘差杠桿圖：從杠桿圖看出，所有的殘差都在0點(diǎn)附近均勻分布，區(qū)間幾乎都位于之間，即，沒(méi)有發(fā)現(xiàn)高

13、杠桿點(diǎn)，也就是說(shuō)，數(shù)據(jù)中沒(méi)有強(qiáng)影響點(diǎn)、異常觀測(cè)點(diǎn)。 綜合起來(lái)看，以上回歸結(jié)果（回歸函數(shù)、擬合曲線或曲面）近乎完美。（二）回歸診斷的必要性例2（Anscombe，1973）數(shù)據(jù)如下：序號(hào)108.049.147.468.06.5828.06.958.146.7785.76313.07.588.7412.7487.71498.818.777.1188.845118.339.267.8188.476149.968.18.8487.04767.246.136.0885.25844.263.105.391912.591210.849.138.1585.561074.827.266.4487.911155.

14、684.745.7386.89解：A 在Matlab7.1中，Import Data：A=上表中的數(shù)據(jù)。（1）畫圖觀察，并做與的線性回歸畫圖觀察x1=A(:,2);y1=A(:,3);plot(x1,y1,-*)旋轉(zhuǎn)圖形，看出數(shù)據(jù)點(diǎn)都在一個(gè)平面上，但不一定在同一條直線上。做線性回歸：x1=A(:,2);y1=A(:,3);plot(x1,y1,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y1,x)rcoplot(r,rint)var(r)程序執(zhí)行后，得到結(jié)果：b = 5.2818 0.2743即，系數(shù)的95%可信區(qū)間為：bint = 2

15、.6131 7.9504 -0.0206 0.5691回歸誤差的方差為：Var= 2.7663擬合指標(biāo)為：stats = 0.3297 4.4277 0.0647 3.0737即，與顯著性概率相關(guān)的，這說(shuō)明，擬合效果不好。殘差杠桿圖為：殘差杠桿圖顯示，第一個(gè)數(shù)據(jù)有異常。（2）畫圖觀察，并做與的線性回歸畫圖觀察x1=A(:,2);y2=A(:,4);plot(x1,y1,-*) 有一個(gè)數(shù)據(jù)點(diǎn)不在直線上。做線性回歸x1=A(:,2);y2=A(:,4);plot(x1,y2,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y2,x)rcoplo

16、t(r,rint)程序執(zhí)行后得到：b = 5.7742 0.2134即，回歸誤差的方差為：var = 3.3036擬合指標(biāo)不太好：stats = 0.1996 2.2447 0.1683 3.6707殘差杠桿圖為：殘差杠桿圖顯示，有2個(gè)數(shù)據(jù)異常。（3）畫圖觀察，并做與的線性回歸畫圖觀察 做線性回歸x1=A(:,2);y3=A(:,5);plot(x1,y1,-*)e=ones(11,1);x=e,x1;b,bint,r,rint,stats=regress(y3,x)rcoplot(r,rint)var(r)結(jié)果如下：b = 5.0264 0.3059即，回歸誤差的方差為：Var=2.4250

17、擬合指標(biāo)不太好：stats = 0.4112 6.2846 0.0335 2.6944（4）畫圖觀察，并做與的線性回歸畫圖觀察 做線性回歸b = 3.0017 0.4999即，回歸誤差的方差為：Var=1.3742擬合指標(biāo)不太好：stats = 0.6667 18.0033 0.0022 1.5269殘差杠桿圖為：B 在Mathematica5.0中，做最小二乘擬合：A=0.0000,8.0400,9.1400,7.4600,8.0000,6.5800,8.0000,6.9500,8.1400,6.7700,8.0000,5.7600,13.0000,7.5800,8.7400,12.7400

18、,8.0000,7.7100,9.0000,8.8100,8.7700,7.1100,8.0000,8.8400,11.0000,8.3300,9.2600,7.8100,8.0000,8.4700,14.0000,9.9600,8.1000,8.8400,8.0000,7.0400,6.0000,7.2400,6.1300,6.0800,8.0000,5.2500, 4.0000,4.2600,3.1000,5.3900,19.0000,12.5000,12.0000,10.8400,9.1300,8.1500,8.0000,5.5600,7.0000,4.8200,7.2600,6.440

19、0,8.0000,7.9100,5.0000,5.6800,4.7400,5.7300,8.0000,6.8900;n=LengthA;x1=TableAj,1,j,n;y1=TableAj,2,j,n;d1=Threadx1,y1;f1=Fitd1,1,x,x;Printf1 = ,f1y2=TableAj,3,j,n;d2=Threadx1,y2;f2=Fitd2,1,x,x;Printf2 = ,f2y3=TableAj,4,j,n;d3=Threadx1,y3;f3=Fitd3,1,x,x;Printf3 = ,f3x2=TableAj,5,j,n;y4=TableAj,6,j,n;d

20、4=Threadx2,y4;f4=Fitd4,1,x,x;Printf4 = ,f4程序執(zhí)行后，得到結(jié)果：f1 = 5.28176 +0. xf2 = 5.77417 +0. xf3 = 5.02645 +0. xf4 = 3.00173 +0. x注：這個(gè)例子說(shuō)明，應(yīng)該對(duì)數(shù)據(jù)擬合、回歸結(jié)果作出評(píng)價(jià)，即，要做回歸診斷。（三）逐步回歸假設(shè)已有數(shù)據(jù)X 和Y，在Matlab軟件包中，使用stepwise命令進(jìn)行逐步回歸，得到回歸方程，其中是隨機(jī)誤差。stepwise命令的使用格式如下：stepwise(X,Y)注意：應(yīng)用stepwise命令做逐步回歸，數(shù)據(jù)矩陣X的第一列不需要人工加一個(gè)全1向量，程序

21、會(huì)自動(dòng)求出回歸方程的常數(shù)項(xiàng)（intercept）。在應(yīng)用stepwise命令進(jìn)行運(yùn)算時(shí)，程序不斷提醒將某個(gè)變量加入（Move in）回歸方程，或者提醒將某個(gè)變量從回歸方程中剔除（Move out）。注釋：使用stepwise命令進(jìn)行逐步回歸，既有剔除變量的運(yùn)算，也有引入變量的運(yùn)算，它是目前應(yīng)用較為廣泛的一種多元回歸方法。在運(yùn)行stepwise(X,Y)命令時(shí)，默認(rèn)顯著性概率。 例2（Hald,1960）Hald數(shù)據(jù)是關(guān)于水泥生產(chǎn)的數(shù)據(jù)。某種水泥在凝固時(shí)放出的熱量（單位：卡克）與水泥中4種化學(xué)成分所占的百分比有關(guān)：在生產(chǎn)中測(cè)得13組數(shù)據(jù)：序號(hào)X1X2X3X4Y172666078.52129155

22、274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.4求出關(guān)系式。解：（1）本問(wèn)題涉及的數(shù)據(jù)是5維的，不能畫圖觀察。先做異常值分析。X=7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23

23、,34;11,66,9,12;10,68,8,12;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;A=X,Y;mahal(A,A)程序執(zhí)行后得到結(jié)果：ans = 5.6803 3.6484 6.7002 3.3676 3.3839 4.4300 4.0080 6.5067 3.0849 7.5016 5.1768 2.4701可以認(rèn)為數(shù)據(jù)都是正常的。（2）一般多元回歸。在Matlab軟件包中寫一個(gè)M文件opt_cement_1:X=7,26,6,60;1,29,15,52;11,56,8,20;

24、11,31,8,47; 7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44; 2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12; 10,68,8,12;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;a1=ones(13,1);A=a1,X;b,bint,r,rint,stat=regress(Y,A)rcoplot(r,rint)程序執(zhí)行后得到：b = 62.4054 1.5511 0.5102 0.1019 -0.1441

25、bint = -99.1786 223.9893 -0.1663 3.2685 -1.1589 2.1792 -1.6385 1.8423 -1.7791 1.4910r = 0.0048 1.5112 -1.6709 -1.7271 0.2508 3.9254 -1.4487 -3.1750 1.3783 0.2815 1.9910 0.9730 -2.2943rint = -4.0390 4.0485 -3.2331 6.2555 -5.3126 1.9707 -6.5603 3.1061 -4.5773 5.0788 -0.5623 8.4132 -6.0767 3.1794 -6.89

26、63 0.5463 -3.5426 6.2993 -3.0098 3.5729 -2.2372 6.2191 -4.1338 6.0797 -6.9115 2.3228stat = 0.9824 111.4792 0.0000 5.9830以及殘差杠桿圖：于是，我們得到：并且，殘差杠桿圖顯示，殘差均勻分布在0點(diǎn)線附近，在stat返回的4個(gè)值中，0.9824，說(shuō)明模型擬合的很好。F_檢驗(yàn)值111.47920.000，符合要求。但是，與顯著性概率相關(guān)的值5.98300.05，這說(shuō)明，回歸方程中有些變量可以剔除。（3）逐步回歸在Matlab軟件包中寫一個(gè)M文件opt_cement_2:X=7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47; 7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44; 2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12; 10,68,8,12;Y=78.5,74.3,104.3,87.6,95.9,10