成信工統(tǒng)計(jì)天氣預(yù)報(bào)實(shí)驗(yàn)指導(dǎo)01回歸分析

統(tǒng)計(jì)天氣預(yù)報(bào)實(shí)驗(yàn)指導(dǎo)實(shí)驗(yàn)一回歸分析一、目的和要求：回歸分析是統(tǒng)計(jì)天氣預(yù)報(bào)中最基礎(chǔ)和最重要的方法，是必須掌握的方法之一。該方法既能作為天氣分析模型，又能作為天氣預(yù)報(bào)模型，在實(shí)際工作中有極強(qiáng)的實(shí)用意義。通過(guò)該實(shí)驗(yàn)，掌握運(yùn)用物理統(tǒng)計(jì)方法建立區(qū)域降水預(yù)測(cè)的基本步驟，運(yùn)用提供的資料和方法子程序，編寫或補(bǔ)充程序當(dāng)中的部分片斷，了解區(qū)域降水的預(yù)報(bào)方法及其建立過(guò)程，輸出實(shí)驗(yàn)要求的相應(yīng)結(jié)果，并就方法對(duì)區(qū)域降水的擬合及實(shí)驗(yàn)預(yù)報(bào)效果進(jìn)行分析，從而深刻理解回歸統(tǒng)計(jì)模型的意義，理解和掌握統(tǒng)計(jì)分析中的核心概念如平均值、方差、協(xié)方差、相關(guān)系數(shù)等的分析計(jì)算，掌握多元回歸模型和逐步回歸模型的計(jì)算步驟、原理和計(jì)算差異、顯著性檢驗(yàn)的意義，為統(tǒng)計(jì)結(jié)果的分析、統(tǒng)計(jì)模型的實(shí)際應(yīng)用打下一定的基礎(chǔ)。二、實(shí)驗(yàn)的主要內(nèi)容：我國(guó)屬于季風(fēng)氣候，降水量的季節(jié)差異很大，年降水的大部分集中在夏半年。影響我國(guó)夏季降水的物理因子很多，根據(jù)多年的研究，在業(yè)務(wù)預(yù)測(cè)中?？紤]的有副熱帶高壓和季風(fēng)活動(dòng)、中高緯阻塞與冷空氣的活動(dòng)、海溫特別是與ENSO相聯(lián)系的赤道太平洋海溫等。這里選用河北省10個(gè)地區(qū)500hPa高度距平和（1月、2～3月、4月），石家莊地面WSW-WNW風(fēng)速合計(jì)（12月下旬一次年1月下旬）作為預(yù)報(bào)因子，采用多元回歸分析和逐步回歸分析，建立區(qū)域降水的預(yù)報(bào)方程。具體如下：2.1利用河北省10個(gè)地區(qū)500hPa高度距平和（1月、2～3月、4月），石家莊地面WSW-WNW風(fēng)速合計(jì)（12月下旬一次年1月下旬）作為預(yù)報(bào)因子；2.2運(yùn)用多元回歸方法，對(duì)1951-1968年河北省6～8月降水量建立預(yù)報(bào)方程；2.3對(duì)1969-1971年進(jìn)行多元回歸預(yù)測(cè)實(shí)驗(yàn)；2.4運(yùn)用逐步回歸方法，對(duì)1951-1968年河北省6～8月降水量建立預(yù)報(bào)方程；2.5對(duì)1969-1971年進(jìn)行逐步回歸預(yù)測(cè)實(shí)驗(yàn)；2.6完成實(shí)驗(yàn)報(bào)告。三、步驟：3.1熟悉資料方法3.1.1資料資料預(yù)報(bào)因子和降水量，依次按列排放，時(shí)間從1951年-1971年，其中因子序列包括1月、2～3月、4月500hPa高度距平和，地面WSW-WNW風(fēng)速合計(jì)共4項(xiàng)；降水量為河北省逐年6～8月。3.1.2多元回歸方法回歸分析是用來(lái)尋找若干變量之間統(tǒng)計(jì)關(guān)系的一種方法，利用所找到的統(tǒng)計(jì)關(guān)系對(duì)某一變量作出未來(lái)時(shí)刻的估計(jì)，稱為回歸預(yù)報(bào)值?；貧w作為一種統(tǒng)計(jì)模型，包括線性和非線性回歸。由于線性回歸比較簡(jiǎn)單，理論上也比較嚴(yán)謹(jǐn)，在氣象中不少氣象要素之間可以近似地存在這種線性的關(guān)系，因此，這種線性回歸應(yīng)用得比較多。預(yù)報(bào)量常常與多個(gè)因子之間存在線性相關(guān)關(guān)系，因而大部分氣象統(tǒng)計(jì)預(yù)報(bào)中的回歸分析都是多元回歸。設(shè)預(yù)報(bào)量與m個(gè)因子的關(guān)系是線性的，則多元線性回歸方程可表示為：其中Y是預(yù)報(bào)量，是預(yù)報(bào)因子，e為殘差誤差，，，，…為估計(jì)的回歸系數(shù)。3.1.3逐步回歸方法在氣象預(yù)報(bào)中，對(duì)預(yù)報(bào)量的預(yù)報(bào)常常需要從可能影響預(yù)報(bào)y的許多因素中挑選一批關(guān)系較好的作為預(yù)報(bào)因子，應(yīng)用多元線性回歸的方法建立回歸方程來(lái)作預(yù)報(bào)。逐步回歸分析是在已選定的一批因子中得到“最優(yōu)”的回歸方程的一種常用方法。以下是幾種實(shí)際計(jì)算方案來(lái)選擇“最優(yōu)”回歸方程。逐步剔除方案：從包含全部變量的回歸方程中逐步剔除不顯著的因子。但此方案計(jì)算量大。逐步引進(jìn)方案：在一批待選的因子中，考察它們對(duì)預(yù)報(bào)量y的方差貢獻(xiàn)，挑選所有因子中方差貢獻(xiàn)最大者，經(jīng)統(tǒng)計(jì)檢驗(yàn)是顯著的，進(jìn)入回歸方程。但此方案不一定保證最后的方程是“最優(yōu)”的。雙重檢驗(yàn)的逐步回歸方案：將因子一個(gè)個(gè)引入，引入因子的條件是該因子的方差貢獻(xiàn)是顯著的。同時(shí)，每引入一個(gè)新因子后，要對(duì)老因子逐個(gè)檢驗(yàn)，將方差貢獻(xiàn)變?yōu)椴伙@著的因子剔除。這種方案是利用求解線性方程中求逆同時(shí)并行的方法。此方案綜合前兩種方案的優(yōu)點(diǎn)克服其缺點(diǎn)，所以本實(shí)驗(yàn)選擇此方案。3.2效果分析3.2.1回歸擬合效果的參數(shù)分析(1)殘差平方和(SSR)(2)標(biāo)準(zhǔn)差(3)復(fù)相關(guān)系數(shù)其中當(dāng)R近似等于1，則相對(duì)誤差將近似0，說(shuō)明回歸效果好。(4)偏相關(guān)系數(shù)其中當(dāng)偏相關(guān)系數(shù)越大，回歸中因子的重要性越大。(5)回歸方差(6)多元回歸的F檢驗(yàn)(7)方程中已引入l個(gè)因子時(shí),其中第k個(gè)因子的方差貢獻(xiàn)(8)對(duì)第k個(gè)因子作F檢驗(yàn)3.2.2預(yù)測(cè)實(shí)驗(yàn)表1采用多元回歸分析河北省6～8月降水量預(yù)測(cè)實(shí)驗(yàn)結(jié)果年份196919701971觀測(cè)值預(yù)測(cè)值表2采用逐步回歸分析河北省6～8月降水量預(yù)測(cè)實(shí)驗(yàn)結(jié)果年份196919701971觀測(cè)值預(yù)測(cè)值3.3編寫程序根據(jù)本實(shí)驗(yàn)的具體情況，參考主程序讀入資料和預(yù)報(bào)量文件，建立回歸方程并輸出相關(guān)的參數(shù)與回歸系數(shù)；然后根據(jù)擬合的系數(shù)，編寫程序計(jì)算1951-1968年擬合估計(jì)值和對(duì)1969-1971年進(jìn)行預(yù)測(cè)實(shí)驗(yàn)的預(yù)測(cè)值。3.3.1多元回歸程序ccccccccccccccccccccccccccccccc多元回歸子程序cccccccccccccccccccccccccccccccccsubroutineDYHG(x,y,m,mm,n,a,q,s,r,v,u,b,dyy) dimensionx(m,n),y(n),a(mm),b(mm,mm),v(m) b(1,1)=n do20j=2,mm b(1,j)=0.0do10i=1,n10 b(1,j)=b(1,j)+x(j-1,i)b(j,1)=b(1,j)20continuedo50i=2,mm do40j=i,mm b(i,j)=0.0 do30k=1,n30b(i,j)=b(i,j)+x(i-1,k)*x(j-1,k)b(j,i)=b(i,j)40continue50continuea(1)=0.0do60i=1,n60a(1)=a(1)+y(i)do80i=2,mm a(i)=0.0 do70j=1,n70a(I)=a(i)+x(i-1,j)*y(j)80continuecallCHOLESKY(b,mm,1,a,l) yy=0.0 do90i=1,n90yy=yy+y(i)/nq=0.0 dyy=0.0 u=0.0ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccdo110i=1,n p=a(1) do100j=1,m100p=p+a(j+1)*x(j,i)q=q+(y(i)-p)*(y(i)-p) dyy=dyy+(y(i)-yy)*(y(i)-yy) u=u+(yy-p)*(yy-p)110continuecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc s=sqrt(q/n) r=sqrt(1.0-q/dyy) do150j=1,m p=0.0 do140i=1,n pp=a(1) do130k=1,m if(k.ne.j)pp=pp+a(k+1)*x(k,i)130continuep=p+(y(i)-pp)*(y(i)-pp)140continuev(j)=sqrt(1.0-q/p)150continuereturn endC!PerformtheCHOLESKYDecompositionsubroutineCHOLESKY(a,n,m,d,l) dimensiona(n,n),d(n,m) l=1 if(a(1,1)+1.0.eq.1.0)then l=0 write(*,30) return endif a(1,1)=sqrt(a(1,1)) do10j=2,n10 a(1,j)=a(1,j)/a(1,1)do100i=2,n do20j=2,i20a(i,i)=a(i,i)-a(j-1,i)*a(j-1,i)if(a(i,i)+1.0.eq.1.0)then l=0 write(*,30) return endif30format(1x,'fail')a(i,i)=sqrt(a(i,i)) if(i.ne.n)then do50j=i+1,n do40k=2,i40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)50a(i,j)=a(i,j)/a(i,i)endif100continuedo130j=1,m d(1,j)=d(1,j)/a(1,1) do120i=2,n do110k=2,i110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)d(i,j)=d(i,j)/a(i,i)120continue130continuedo160j=1,m d(n,j)=d(n,j)/a(n,n)do150k=n,2,-1 do140i=k,n140d(k-1,j)=d(k-1,j)-a(k-1,i)*d(i,j)d(k-1,j)=d(k-1,j)/a(k-1,k-1)150continue160continuereturn end 參考主程序(范例)C樣本長(zhǎng)度：N；因子個(gè)數(shù)：K；PROGRAMMAININTEGER,PARAMETER::N=18 INTEGER,PARAMETER::K=4REAL,DIMENSION(K,N)::X REAL,DIMENSION(N)::Y REAL,DIMENSION(K+1)::A REAL,DIMENSION(K+1,K+1)::B REAL,DIMENSION(K)::V REALQ,S,R,UCOPENTHEINPUTDATAFILEOPEN(10,FILE='x-1.TXT')OPEN(20,FILE='y.TXT')CREADTHEDATADOI=1,N READ(10,*)X(1,I),X(2,I),X(3,I),X(4,I)READ(20,*)Y(I) ENDDO CLOSE(10) MM=K+1C調(diào)用回歸子程序callDYHG(X,Y,K,MM,N,A,Q,S,R,V,U,B,DYY) write(*,88)A(1)88format(/1x,'b0='f9.5)do89j=2,MM89write(*,100)j-1,A(j)100format(1x,'b',i2,'=',f9.5)ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccwrite(*,20)Q,S,R20format(1x,'Q=',f13.6,3x,'S=',f13.6,3x,'R=',f13.6)write(*,22)U,DYY22format(1x,'U=',f13.6,3x,'DYY=',f13.6)write(*,30)(i,V(i),i=1,K)30format(1x,'V(',i2,')=',f13.6)write(*,40)U40format(1x,'U=',f13.6)C輸出回歸系數(shù)及相關(guān)參數(shù)open(6,file='table.txt') write(6,180)180format(/2x,'regressioncoefficients:')write(6,88)A(1) do189j=2,MM189write(6,100)j-1,A(j)write(6,200)200format(1x,'GenericAnalysisofVarianceTablefortheMultiple&LinearRegression')write(6,202)202format(1x,'-----------------------------------------------------------&-------------------') write(6,204)204format(3x,'SourcedfSSMS')write(6,202) write(6,206)N-1,DYY206format(1x,'Totaln-1=',i2,'SST=',f13.4)u2=U/real(K) write(6,280)K,U,U2280format(1x,'RegressionK=