短期氣候?qū)嵙?xí)報告六

1、南京信息工程大學(xué)短期氣候?qū)嶒灒▽嵙?xí)）報告實驗名稱 夏季區(qū)域降水的定量預(yù)測 日期2017.6.6得分指導(dǎo)教師 系大氣科學(xué)專業(yè)大氣科學(xué)年級2014級班次2姓名車楚玉學(xué)號20141301043一、目的要求：目的：掌握短期氣候預(yù)測中物理統(tǒng)計預(yù)測的基本步驟。要求：能運(yùn)用提供的資料和方法子程序，編寫或補(bǔ)充完成程序當(dāng)中的部分片斷，了解區(qū)域降水的預(yù)測方法及其建立過程，輸出實驗要求的相應(yīng)結(jié)果，并就方法對區(qū)域降水的擬合及試驗預(yù)測效果進(jìn)行分析。二、資料和方法：資料：1、前期1月的nino3.4指數(shù)（來自cpc）;2、西太平洋副高脊線、西太平洋副高西伸脊點(diǎn)、亞洲極渦面積、南方濤動指數(shù)（來自中國氣象局整編的 74個環(huán)流

2、指數(shù)）；3、夏季華北區(qū)域10站的降水量，距平百分率。方法：回歸分析（mregrssion.for）是用來尋找若干變量之間統(tǒng)計關(guān)系的一種方法，利用所找 到的統(tǒng)計關(guān)系對某一變量作出未來時刻的估計，稱為回歸預(yù)報值。三、實習(xí)步驟：說明所用資料方法； 計算方法的簡單介紹； 輸出反映回歸效果的參數(shù)及回歸系數(shù)，并就相關(guān)參數(shù)分析回歸效果；預(yù)測量與回歸方程計算的估計值和觀測值的歷年曲線變化圖（19522001年），并附簡單的說明；輸出獨(dú)立預(yù)測試驗的觀測與預(yù)測值。四、結(jié)果（圖表并解釋說明）程序如下：programainintegerparameter=50integerparameter=5 realdimens

3、io(k,n):x realdimensionj):y realdimensionk+1):a realdimensionk+1,k+1):b realdimensionk):vrealq,s,r,u,ind(6,60),year(50),y1(50),pre(7) integer i,jc x(5,50) y(50) a(6) b(6,6),v(5)c open the input data fileopen0,file= 'd:duanqishixi6shixi.txt')c read the data and give data to x and y do i=1,50

4、|read (10,pre(1)=a(1)+a(2)*26.50+a(3)*13.00+a(4)*100.00+a(5)*196.00+a(6)* 180 format (/2x, 'regression coefficients:' )write (6,88) a(1)do 189 j=2,mm189 write (6,100) j-1,a(j)write (6,200)200 format (/1x, 'generic analysis of variance table for the multiple) year(i),x(1:5,i),y(i)enddomm=

5、k+1call dyhg(x,y,k,mm,n,a,q,s,r,v,u,b,dyy) write (*,88) a(1)88 format (/1x, 'b 0=' ,f19.5) do 89 j=2,mm89write (*,100) j-1,a(j)100format (1x, 'b' ,i2, '=' ,f9.5)open(11,file= 'd:duanqishixi6guji.txt' ) open(12,file= 'd:duanqishixi6guji.grd' ,form= 'binary&

6、#39;) open(13,file= 'd:duanqishixi6prediction.txt')do j=1,50y1(j)=a(1)+a(2)*x(1,j)+a(3)*x(2,j)+a(4)*x(3,j)+a(5)*x(4,j)+a(6)* *x(5,j)enddopre(4)=a(1)+a(2)do 140 i=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)150continue160continuereturnend27.10+a(3)*13.00+a(4)*110.00+a(5

7、)*220.00+a(6)*3.00 pre(5)=a(1)+a(2)*25.64+a(3)*14.00+a(4)*105.00+a(5)*215.00+a(6)*13.00 pre(6)=a(1)+a(2)*27.26+a(3)*15.00+a(4)*90.00+a(5)*185.00+a(6)*(-8.00)pre(7)=a(1)+a(2)*24.71+a(3)*17.00+a(4)*125.00+a(5)*222.00+a(6)*13.00 print *,prewrite (11,*)(y(i),y1(i),i=1,50)write (12)(y(i),y1(i),i=1,50)wri

8、te (13,*)(pre(i),i=1,7)ccccccccccccccccccccccccccccccccccccccccccccccccccccccwrite (*,20)q,s,r20 format (1x, 'q=' ,f13.6,3x,'s=' ,f13.6,3x,'r=' ,f13.6)write (*,22)u,dyy22 format (1x, 'u=' ,f13.6,3x, 'dyy=' ,f13.6) write (*,30)(i,v(i),i=1,k)30 format (1x, '

9、v(' ,i2, ')=' ,f13.6) write (*,40)u40 format (1x, 'u=' ,f13.6)open(6,file= 'table' )! output datawrite (6,180)f=(u/ real (k)/(q/ real (n-k-1) write (6,220) f220 format (/1x, 'f=msr/mse=',f13.4)write (6,202) close (6) stop endsubroutine dyhg(x,y,m,mm,n,a,q,s,r,v,u,

10、b,dyy) dimension x(m,n),y(n),a(mm),b(mm,mm),v(m) b(1,1)=ndo 20 j=2,mm b(1,j)=0.0 do 10 i=1,n 10b(1,j)=b(1,j)+x(j-1,i)b(j,1)=b(1,j)20 continue do 50 i=2,mm do 40 j=i,mm b(i,j)=0.0do 100 j=1,m |100p=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)110 continuecccccccc

11、cccccccccccccccccccccccccc s= sqrt (q/n) r= sqrt (1.0-q/dyy) | do 150 j=1,m |p=0.0do 140 i=1,n pp=a(1)do 130 k=1,mj if (k.ne.j)pp=pp+a(k+1)*x(k,i)130 continuep=p+(y(i)-pp)*(y(i)-pp)140continuev(j)= sqrt (1.0-q/p)endif30 format (1x, 'fail') a(i,i sqrt (a(i,i) if (i.ne.n) then do 50 j=i+1,n do

12、 40 k=2,i40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)50a(i,j)=a(i,j)/a(i,i)endif100 continue do 130 j=1,m. d(1,j尸d(1,j)/a(1,1) do 120 i=2,n do 110 k=2,i110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)d(i,j)=d(i,j)/a(i,i)120continue130continuedo 160 j=1,m d(n,j)=d(n,j)/a(n,n)do 150 k=n,2,-1表一 1952-2001年回歸方程計算的估計值和觀測值的歷年曲線變化圖

13、2002200320042005200620072008觀測值-29.75-17.450.64-10.63-22.96-19.403.30預(yù)測值30.8053339.2322031.9930128.4085024.2936114.32182-2.514375表二 2002-2008年預(yù)測試驗的觀測與預(yù)測值五、結(jié)果討論:分析表一，白色實線為 1952-2001年的觀測真實值，綠色實線為1952-2001年的回歸方程計算估計值，由表可以看出，某些年份，預(yù)測值與觀測值的趨勢大致相同，而在某些年份趨勢相反，如1959年為極大值，但是估計值卻為極小值，且各個年份的兩個值的差別也很大，可以看出，通過回歸方

14、程預(yù)測的值不是很理想。通過表二可以看出，預(yù)報觀測值與回歸方程的預(yù)測值差別十分顯著，進(jìn)一步表明了多元線形回歸方程的預(yù)測值并不理想。id* 2.00pre(2)=a(1)+a(2)*27.76+a(3)*12.00+a(4)*90.0+a(5)*197.00+a(6)*(-1.00) pre(3)=a(1)+a(2)*26.74+a(3)*12.0+a(4)*120.00+a(5)*199.00+a(6)* *(-11.00)* linear regression' )write (6,202)202 format (/1x,'*')write (6,204)204 for

15、mat (/3x, 'source df ss ms')write (6,202)write (6,206) n-1,dyy206 format (/1x, 'total n-1=',i2, ' sst=' ,f13.4)u2=u/ real (k)write (6,208) k,u,u2208 format (/1x, 'regression k=',i2, ' ssr=' ,f13.4, ' msr=ssr/k='*,f13.4)q2=q/ real (n-k-1)write (6,209)n-

16、k-1,q,q2209 format (/1x, 'residual n-k-1=',i2, ' sse=' ,f13.4, ' mse=sse/(n-k-1)*=' ,f13.4)do 30 k=1,n30b(i,j尸b(i,j)+x(i-1,k)*x(j-1,k)b(j,i)=b(i,j)40continue50 continue a(1)=0.0do 60 i=1,n60a(1)=a(1)+y(i)do 80 i=2,mm a(i)=0.0do 70 j=1,n70a(i)=a(i)+x(i-1,j)*y(j)80continuecall cholesky(b,mm,1,a,l) yy=0.0do 90 i=1,n90 yy=yy+y(i)/n q=0.0 dyy=0.0 u=0.0cccccccccccccccccccccccccccccccccc do 110 i=1,n p=a(1)150continuereturnendsubroutine cholesky(a,n,m,d,l) ! perform the chole