《應(yīng)用回歸分析》

模型選擇問題:對(duì)于模型y=B+Bx+Bx+Bx+e，其中B=1,B=2,B=—1,B=0TOC\o"1-5"\h\z01122330123用隨機(jī)數(shù)的方法產(chǎn)生n=40組數(shù)據(jù)，要求x?U［—10,10］,k=1,2,3,i=1,…，n；e?N(0,1)；iki并且y由iy=B+Bx+Bx+Bx+ei01i12i23i3i得出。對(duì)于這40組隨機(jī)數(shù)據(jù)(y,x,x,x),i=1,…,n，我們建立了以下四ii1i2i3種模型：^①.y=B+Bx+e011y=B+Bx+Bx+eTOC\o"1-5"\h\z01122y—B+Bx+Bx+e01133④.y=B+Bx+Bx+Bx+e0112233運(yùn)用我們所學(xué)的模型選擇的準(zhǔn)則在①?④中選出最佳模型。一、產(chǎn)生隨機(jī)數(shù)對(duì)于這個(gè)問題，我們首先要解決的是根據(jù)原模型及給定的參數(shù)分布產(chǎn)生問題要求的40組隨機(jī)數(shù)(y,x,x,x),i=1,…,n。ii1i2i3我們知道在Matlab中，可以利用R=rand這個(gè)函數(shù)來產(chǎn)生一個(gè)［0,1］上的隨機(jī)數(shù)，并且R是來自［0,1］的均勻分布，即R?U［0,1］；我們利用R=rand(n,k)就可以得到一個(gè)n行k列的來自均勻分布U［0,1］的隨機(jī)數(shù)組成的矩陣。由此我們可以想到，利用R=10—20*rand(40,3)，我們就可以得到x,k=1,2,3,i=1:40,ik我們?cè)谒淖髠?cè)加入全為1的一列，保存在X中。我們要運(yùn)用林德貝格-勒維中心極限定理通過均勻分布U［0,1］的隨機(jī)數(shù)來產(chǎn)生N(0,1)上的隨機(jī)數(shù)。U［0,1］的期望和方差分別為1/2和1/12,所以12個(gè)相互獨(dú)立的U［0,1］和的期望和方差分別為6和1。因此只要產(chǎn)生12個(gè)U［0,1］上的隨機(jī)數(shù)x,x，…,x,計(jì)算x+x++x—6就得到一個(gè)來自N(0,1)的隨機(jī)數(shù)。1212121211118888"18.289689.389186.53223_"0.21794"8.40813_14.75036-4.881492.18124-1.7277613.654431-6.02029-0.00045-6.627590.71993-10.320219.415590.40156-6.067291.7396221.169251-8.57708-8.094448.790580.12487-7.934851-4.60662-2.197332.014840.80141-5.2145010.22782-2.35333-0.537520.553244.362211-1.57050-7.188851.664010.948755.9965915.25433-6.10979-3.13720-1.9320615.6863810.82302-1.53443-2.55947-0.817933.362551-9.261776.341554.160321.32851-22.53651-0.936115.201361.36698-1.70398-7.777561-0.42272-7.730249.69026-0.277337.6074815.368119.42652-9.681271.015223.3249310.222050.201976.656631.034102.276221-2.481206.641467.87567-0.51016-11.11401-3.58271-9.573612.551810.801244.2094312.08970-4.253896.037630.432239.8655112.65127-0.009430.20625-0.442505.869471-9.759640.578233.210130.58418v-18.513319.245228.80762-9.03261e—0.44931Y—11.132131-7.70336-3.63944-8.40664-0.06311-10.83031-8.265749.151388.94646-0.48189-25.16471-5.923688.57109-4.75716-2.44547-21.863918.02575-0.433004.61761-1.7166015.7679C14.762588.065401.54329-1.279771.1799813.29286-6.36297-0.957421.0754615.024161-3.59456-6.35094-8.854740.713570.8753917.26894-4.448791.64512-0.1250319.861641-4.424557.00269-9.66105-0.39874-15.250517.86476-3.192113.970901.7040021.625631-3.07515-0.37190-4.02198-0.33701-5.1154010.11652-9.45949-3.326781.3918812.084411-5.58103-2.97983-0.782530.02059-7.161651-4.30074-6.00661-3.96211-1.63332-3.228191-8.074410.92405-3.33056-1.70961-17.78241-7.818451.352176.43735-1.05981-17.048813.31674-6.506287.43971-1.9170612.222701-3.974928.33060-9.981610.11329-15.167116.043807.336586.57758__-1.9456_1_3.80542_83X=413因此我們得到了40組數(shù)據(jù)(x,x,x,e),i=1,…,40，將其代入模型i1i2i3iy=B+Bx+Bx+Bx+ei01i12i23i3i就得到了上頁中以矩陣形式表示的40組隨機(jī)數(shù)(y,x,x,x),i=1,…,40。ii1i2i3二、模型選擇準(zhǔn)則這里我們有五種模型選取準(zhǔn)則：1、平均平方和準(zhǔn)則對(duì)于一個(gè)選模型，假設(shè)模型中含有p個(gè)回歸變量，記：MS=1SSEpn—pp其中SSE是在此選模型下的殘差平方和。計(jì)算多個(gè)選模型的MS，我們認(rèn)為ppMS越小的模型效果越好。p2、Gp準(zhǔn)則同樣的，我們對(duì)選模型計(jì)算：SSEG=——n+2pp￡2其中￡2是全模型下的￡2的最小二乘估計(jì)。G越小，模型效果越好。p3、AIC準(zhǔn)則Y,Y，…，Y是一個(gè)樣本，記含有k個(gè)參數(shù)的模型的似然函數(shù)為TOC\o"1-5"\h\z12nL(9IY，…，Y),9的MLE為J，則AIC準(zhǔn)則要求1kAIC=ln厶(0八IY，…，Y)—k1k的值越大，選模型的效果越好。進(jìn)一步地，在線性模型場(chǎng)合，我們有nAIC=—lnSSE+p2p的值越小越好。4、CV準(zhǔn)則將40組原始數(shù)據(jù)的第i組數(shù)據(jù)刪去，利用剩下的39組數(shù)據(jù)對(duì)選模型進(jìn)行最

小二乘估計(jì)，將第i組數(shù)據(jù)(x,x,x)代入模型中得出y。對(duì)i=l,2,…，40ili2i3i重復(fù)進(jìn)行上述操作40次，最后計(jì)算CV=-￡(y-y)2nzzi=1CV越小，選模型效果越好。5、BIC準(zhǔn)則SSEBIC=—+plogn其中&2是全模型下的b2的最小二乘估計(jì)，BIC越小，選模型效果越好。三、模型選擇在以上幾種準(zhǔn)則中需要用到全模型下的一些數(shù)據(jù)，所以我們先就全模型即第④種模型進(jìn)行分析。1、全模型y=卩+卩x+卩x+卩x+e0112233將所有數(shù)據(jù)導(dǎo)入到Minitab軟件中，可以得到:0.83392.00381-1.035982.00381-1.03598，SSE=49.5，b2=1.37569-0.02939由此,y=0.8339+2.00381x-1.03598x-0.02939x1231SSE=MS1SSE1SSE=MSpn-ppn-pSSE丄-n+2p=1.98183nAIC=lnSSE+p=81.039452p在Matlab中利用循環(huán)可以求得CV,定義一個(gè)nx1階的Y1用以保存每次得到的y，并且輸入如下循環(huán)語句：i>>fori=1:40A=X;B=Y;A1=A(i,:);B1=B(i,:);A(i,:)=[];B(i,:)=[];R=regress(B,A);YO=A1*R;Yl(i,l)=YO;A=X;B=Y;end于是得到：1“CV=￡(y-y)2=1.52538niii=1SSEBIC=p+plogn=40.78801b22、選模型①y=卩+卩x+e011將X的第3、4列刪去，然后和上面一樣我們可以得到:入「0.961]0=，SSE=1566.3丄9630」p由此，y=0.961+1.9630x丿1MS=—^SSE=40.16154pn—ppSSEG=p—n+2p=1100.552pb2nAIC=尹SSEp+p=148-12941“CV=￡(y—y)2=43.27734(只需將上述循環(huán)中的第二行改為niii=1A=X(:,[12]);B=Y;即可)SSEBIC=p+plogn=1140.154b23、選模型②y=0+0x+0x+e01122刪去X中的第4列，進(jìn)行回歸，得到：ii=10.8281P二2.00337,SSE=50.7-1.03221」卩所以y二0.8281+2.00337x—1.03221x12MS=1SSE=1.33421pn—ppSSEG=p—n+2p=0.85412Pb2nAIC=—InSSE+p=80.518522p1“CV=—^(y—y)2=1.50043niii=1SSEBIC=p+plogn=40.05823b24、選模型③y=卩+卩x+卩x+eTOC\o"1-5"\h\z01133刪去X中的第3列，用同樣的方法回歸，得：_0.937_P=1.9619，SSE=1549.90.1101J"所以y=0.937+1.9619x+0.1101x13MS=—SSE=40.78684pn—ppSSEGp=￡—n+2p=1090.6310nAIC=—lnSSE+p=148.91892pCV=1￡(y—y)2=45.7901niiSSEBIC=p+plogn=1129.835b2四、結(jié)論將上述四種模型計(jì)算所得的MS,G,AIC,CV,BIC數(shù)據(jù)統(tǒng)計(jì)到同一表pp格中進(jìn)行直觀比較。MSpGpAICCVBIC模型140.161531100.552148.129443.277341140.154模型21.334210.8541180.518511.5004340.