成信工统计天气预报实验指导01回归分析

统计天气预报实验指导实验一回归分析一、目的和要求：回归分析是统计天气预报中最基础和最重要的方法，是必须掌握的方法之一。该方法既能作为天气分析模型，又能作为天气预报模型，在实际工作中有极强的实用意义。通过该实验，掌握运用物理统计方法建立区域降水预测的基本步骤，运用提供的资料和方法子程序，编写或补充程序当中的部分片断，了解区域降水的预报方法及其建立过程，输出实验要求的相应结果，并就方法对区域降水的拟合及实验预报效果进行分析，从而深刻理解回归统计模型的意义，理解和掌握统计分析中的核心概念如平均值、方差、协方差、相关系数等的分析计算，掌握多元回归模型和逐步回归模型的计算步骤、原理和计算差异、显著性检验的意义，为统计结果的分析、统计模型的实际应用打下一定的基础。二、实验的主要内容：我国属于季风气候，降水量的季节差异很大，年降水的大部分集中在夏半年。影响我国夏季降水的物理因子很多，根据多年的研究，在业务预测中常考虑的有副热带高压和季风活动、中高纬阻塞与冷空气的活动、海温特别是与ENSO相联系的赤道太平洋海温等。这里选用河北省10个地区500hPa高度距平和（1月、2～3月、4月），石家庄地面WSW-WNW风速合计（12月下旬一次年1月下旬）作为预报因子，采用多元回归分析和逐步回归分析，建立区域降水的预报方程。具体如下：2.1利用河北省10个地区500hPa高度距平和（1月、2～3月、4月），石家庄地面WSW-WNW风速合计（12月下旬一次年1月下旬）作为预报因子；2.2运用多元回归方法，对1951-1968年河北省6～8月降水量建立预报方程；2.3对1969-1971年进行多元回归预测实验；2.4运用逐步回归方法，对1951-1968年河北省6～8月降水量建立预报方程；2.5对1969-1971年进行逐步回归预测实验；2.6完成实验报告。三、步骤：3.1熟悉资料方法3.1.1资料资料预报因子和降水量，依次按列排放，时间从1951年-1971年，其中因子序列包括1月、2～3月、4月500hPa高度距平和，地面WSW-WNW风速合计共4项；降水量为河北省逐年6～8月。3.1.2多元回归方法回归分析是用来寻找若干变量之间统计关系的一种方法，利用所找到的统计关系对某一变量作出未来时刻的估计，称为回归预报值。回归作为一种统计模型，包括线性和非线性回归。由于线性回归比较简单，理论上也比较严谨，在气象中不少气象要素之间可以近似地存在这种线性的关系，因此，这种线性回归应用得比较多。预报量常常与多个因子之间存在线性相关关系，因而大部分气象统计预报中的回归分析都是多元回归。设预报量与m个因子的关系是线性的，则多元线性回归方程可表示为：其中Y是预报量，是预报因子，e为残差误差，，，，…为估计的回归系数。3.1.3逐步回归方法在气象预报中，对预报量的预报常常需要从可能影响预报y的许多因素中挑选一批关系较好的作为预报因子，应用多元线性回归的方法建立回归方程来作预报。逐步回归分析是在已选定的一批因子中得到“最优”的回归方程的一种常用方法。以下是几种实际计算方案来选择“最优”回归方程。逐步剔除方案：从包含全部变量的回归方程中逐步剔除不显著的因子。但此方案计算量大。逐步引进方案：在一批待选的因子中，考察它们对预报量y的方差贡献，挑选所有因子中方差贡献最大者，经统计检验是显著的，进入回归方程。但此方案不一定保证最后的方程是“最优”的。双重检验的逐步回归方案：将因子一个个引入，引入因子的条件是该因子的方差贡献是显著的。同时，每引入一个新因子后，要对老因子逐个检验，将方差贡献变为不显著的因子剔除。这种方案是利用求解线性方程中求逆同时并行的方法。此方案综合前两种方案的优点克服其缺点，所以本实验选择此方案。3.2效果分析3.2.1回归拟合效果的参数分析(1)残差平方和(SSR)(2)标准差(3)复相关系数其中当R近似等于1，则相对误差将近似0，说明回归效果好。(4)偏相关系数其中当偏相关系数越大，回归中因子的重要性越大。(5)回归方差(6)多元回归的F检验(7)方程中已引入l个因子时,其中第k个因子的方差贡献(8)对第k个因子作F检验3.2.2预测实验表1采用多元回归分析河北省6～8月降水量预测实验结果年份196919701971观测值预测值表2采用逐步回归分析河北省6～8月降水量预测实验结果年份196919701971观测值预测值3.3编写程序根据本实验的具体情况，参考主程序读入资料和预报量文件，建立回归方程并输出相关的参数与回归系数；然后根据拟合的系数，编写程序计算1951-1968年拟合估计值和对1969-1971年进行预测实验的预测值。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样本长度：N；因子个数：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调用回归子程序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输出回归系数及相关参数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=