气象统计方法实习报告

1、实习一求500hPa高度场气候场、距平场和均方差场 31、资料介绍 32要求 33、实习结果 31）、 FORTRAN 源程序 3（2）、 grads 文件 5（3） 、实习结果 6实习二 计算给定数据资料的简单相关系数和自相关系数 71、资料介绍 72、要求 73、实习结果 7（ 1）、 Fortran 源程序 7（ 2 ）、程序运行结果： 10实习三 分析中国夏季降水线性趋势的分布特征 101. 资料介绍及要求： 102. 实习结果 10（1） .matlab 程序 10（ 2） .程序运行结果 11实习四 求给定数据的一元线性回归方程 11111 、资料介绍及要求2、实习结果 12（1）

2、、 MATLAB 程序 12（ 2 ）、程序运行结果 12（ 3）、结果分析 13实习五 对给定的海温数据进行 EOF 分析131、资料介绍 132、要求 143、实习结果： 14（1）、 FORTRAN 源程序 14（ 2）空间场和时间序列的 ctl 文件 15（ 3）运行结果 16（ 4）分析 16实习三（附加） 计算给定数据的 1 1年滑动平均和累积距平 171、资料介绍 172、要求 173、实习结果 18实习四（附加） 求给定数据的多元线性回归方程 191、说明 192、要求 193、实习结果： 20（1） Matlab 源程序 20（ 2）运行结果 23实习一 求 500hPa 高

3、度场气候场、距平场和均方差场1、资料介绍有一 500hPa高度场资料，文件名，范围： 60150E, 040N.时段： 共 48个月。水平分辨率： *，格点数： 37*17。2要求编fortran程序，求500hPa高度场的( 1 )气候场；(2)距平场；( 3)均方差场。并能用Grads做出图形，实习报告中气候场、距平场、均方差场任意给出 两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。注：给出了如何用 fortran 读取 ASCII 码资料.3、实习结果1) 、FORTRA源程序program ex_gradsimplicit noneinteger,parameter:nx

4、=37,ny=17,nz=4,nt=12integer i,j,iz,itreal var(nx,ny,nz,nt),cl(nx,ny,nt),sum,jp(nx,ny,nz,nt),jfc(nx,ny,nt)! Opening fileopen(10,file='g:gradsdata')do iz=1,nzdo it=1,ntread(10,1000)read(10,3000) (var(i,j,iz,it),i=1,nx),j=1,ny)enddoenddo1000 format(2i7)2000 format3000 format4000 formatclose(10)

5、!Outputopen(16,file='g:gradsdata',form='binary') do iz=1,nzdo it=1,ntwrite(16) (var(i,j,iz,it),i=1,nx),j=1,ny)enddoenddo!Calculating the Climatological Fielddo it=1,ntdo i=1,nxdo j=1,nysum=0do iz=1,nz sum=sum+var(i,j,iz,it) enddo cl(i,j,it)=sum/4 enddo enddo enddo! Output climate-fil

6、e open(12,file='g:gradsdata',form='binary') do it=1,ntwrite(12) (cl(i,j,it),i=1,nx),j=1,ny) enddo!Calculating the Anomalydo iz=1,nzdo it=1,ntdo i=1,nxdo j=1,ny jp(i,j,iz,it)=var(i,j,iz,it)-cl(i,j,it) enddo enddo enddo enddo open(13,file='g:gradsdata',form='binary') !O

7、utput anomaly-filedo iz=1,nzdo it=1,ntwrite(13) (jp(i,j,iz,it),i=1,nx),j=1,ny)enddoenddo !Calculating the Mean-square Deviation do it=1,ntdo i=1,nxdo j=1,nysum=0do iz=1,nz sum=sum+(jp(i,j,iz,it)*2 enddojfc(i,j,it)=sqrt(sum/4)enddoenddoenddo!Output mean-square deviation-file open(14,file='g:grads

8、data',form='binary') do it=1,ntwrite(14) (jfc(i,j,it),i=1,nx),j=1,ny)enddoend(2) 、 grads 文件'open g:gradsdata*.ctl' (* 为所求场对应的 ctl 文件名 ) 'set lat 0 40''set lon 60 150''set lev 500''enable print g:gradsdata*.gmf'( * 为所求场名称)i=1while(i<=48( 或 12)

9、9;set t 'i'd h''print''c'i=i+1 endwhile 'disable print'3）、实习结果 、原始场1982年1月 1982年7月 结果分析：冬季（此处以 1月为代表）等高线分布整体平缓，表明高度场分布相对均匀，且北部接近极地位势 高度低，赤道地区位势高度高，这与太阳直射点在1月在南半球，使北半球整体的辐射吸收随纬度增加而减小有关。北半球的气块受热随纬度递增而递减，因而膨胀率递减，故南方接近赤道地区的气体膨胀 大，位势高，而北方近极地地区气体膨胀小，位势低。夏季（此处以 7 月为代表），大

10、洋上空出现副热带高压（588线位置），在东亚地区存在较为明显的位势高度槽，即东亚大槽。夏季在30° N以北的地区位于西风带中，从图中可看岀明显的西风带长波特征。除东亚大槽外，在中亚地区也存在一长波槽，这些槽线发生长波调整时，会在部分地区发生较剧烈 的天气变化。此外，从图中可以看岀，1982年7月副热带高压脊线的平均位置位于25° N ,125 ° E附近，我国华东地区位于副高北侧西南气流控制下，西南气流为水汽输送的主要通道，该地区发生降水较为频 繁。从图中还可以看岀冬夏季的位势高度分布存在明显差异，这与太阳直射点的年纪变化密切相关。 、气候场2月8月结果分析：气候

11、场是多年数据中同时段的平均值序列，表征了区域内多年平均的位势高度变化。从气候场图形可以看岀多年平均的500hPa高度场中，冬季（此处以多年平均的2月气候场为代表）等高线较平直，大陆上等高线较稀疏，而海洋上等高线较密集，这表明大陆上空的位势高度变化率较海 洋小。此外，冬季等位势高度线分布平直，还说明冬季的天气变化显着。夏季（此处以多年平均的 8月气候场为代表）在太平洋上有副热带高压，副高是深厚的系统，所以 可以一直延伸到500hPa高度。东亚地区存在明显的西风带长波槽，即东亚大槽，东亚大槽的移动和变化 配合副热带高压和夏季风的影响，会使我国大部分地区天气造成强烈变化，如形成大范围降水或强对流 天

12、气。同时，孟加拉湾处存在一低涡。由于高原的阻挡作用，这一系统对我国的影响并不显着。图中副 高脊线8月的多年平均位置位于 30 ° N以北，且东部长波槽位于 110 ° E附近，故由气候场分析，华北地区 位于长波槽前，又槽线受到副高的阻挡作用，因而华北地区容易形成降水。 、距平场1982年1月1982年7月结果分析：距平场指示了位势高度的震荡趋势，因距平的平均值为0，则大于 0的值表明位势高度偏高，小于 0的值表示位势高度偏低。从图中可看岀冬季（此处以 1982年1月距平场做代表）在大陆位势高度为正距平，而在大洋则存在 明显的负距平。则由距平场的性质得，冬季在大洋上位势高度偏

13、高，在大陆上位势高度偏低。其原因 是，海水的比热大于陆地，则冬季海洋温度比陆地高，所以海洋上气块膨胀更多，位势高度更高。夏季（此处以 1982年7月距平场做代表）相反，大洋上位势高度偏低，而陆地偏高。在70° E90° E附近位势偏高的原因在于夏季青藏高原接受太阳辐射使之在对流层中层形成热源，位势高度因比大陆其他位置高。由此可见，位势高度的变化不仅与海陆差异有关，同时也与地形有关。在海洋上副热带高压 所在的位置存在证据平值。、均方差场6月12月结果分析 ：均方差场反应同一时段内的位势高度变化幅度的大小。 由图可以看出，整体位势高度在大陆上的变化幅度比海洋小， 且海洋上冬季的

14、变化幅度比夏季 大，而陆地上相反，冬季的变化幅度比夏季小。因为陆地的比热小，所以陆地在夏季白天与夜间的温差 大于冬季，对应的高度场震荡就比冬季剧烈。海洋上的位势高度变化幅度的影响方面温差为次要因素， 其主要受到副热带高压，西风带长波槽脊影响，位势高度根据天气系统的移动而变化，所以震荡幅度较 大。除受到天气系统影响外，海洋上的位势高度场还受到洋流等因素的影响。实习二 计算给定数据资料的简单相关系数和自相关系数1、资料介绍根据下表中年平均气温和冬季平均气温的等级数据进行下列计算：1）计算两个气温之间的简单相关系数。2）分别找出两个气温数据自相关系数绝对值最大的滞后时间长度。 （滞后长度T最大取10

15、）2、要求实习报告中附出简单相关系数或自相关系数程序。答案：r=年平均气温在滞后长度j=3、冬季序列在j=4最大。3、实习结果（ 1）、 Fortran 源程序PROGRAM EXAMIMPLICIT NONEINTEGER,PARAMETER:N=20INTEGER i,j,k,ty,tw,tywREAL:avr_y=0,avr_w=0,sy=0,sw=0,rxy=0,max_y=0,max_w=0,max_yw=0REAL y(N),w(N)DATA y/,DATA w/,REAL syy(N),sww(N),r(N),rty(N),rtw(N),rtyw(N),rxy_ty(N),rxy

16、_tw(N),rxy_tyw(N) !求两数组平均值DO i=1,Navr_y=avr_y+y(i) avr_w=avr_w+w(i)END DOavr_y=avr_y/N avr_w=avr_w/N!简单相关系数DO j=1,N syy(j)=(y(j)-avr_y)*2 sy=sy+syy(j)sww(j)=(w(j)-avr_w)*2 sw=sw+sww(j)END DOsy=sqrt(sy/N) sw=sqrt(sw/N)DO j=1,Nr(j)=(y(j)-avr_y)/sy)*(w(j)-avr_w)/sw)rxy=rxy+r(j)END DOrxy=rxy/NPRINT &quo

17、t;(/'1970-1989 年全年平均气温与冬季平均气温的简单相关系数 rxy=',",rxy k=0!自相关系数DO ty=1,N/2DO i=1,N-tyrty(i)=(y(i)-avr_y)/sy)*(y(i+ty)-avr_y)/sy) rxy_ty(ty)=rxy_ty(ty)+rty(i)END DOrxy_ty(ty)=rxy_ty(ty)/(N-ty)rxy_ty(ty)=ABS(rxy_ty(ty)IF(rxy_ty(ty)>max_y) THEN max_y=rxy_ty(ty) k=tyEND IFEND DOPRINT "(&

18、#39; 全年平均气温绝对值最大自相关系数 rxy_ty=',/,' 滞后时间长度 k=0DO tw=1,N/2DO i=1,N-tw rtw(i)=(w(i)-avr_w)/sw)*(w(i+tw)-avr_w)/sw) rxy_tw(tw)=rxy_tw(tw)+rtw(i)END DO rxy_tw(tw)=rxy_tw(tw)/(N-tw) rxy_tw(tw)=ABS(rxy_tw(tw)IF(rxy_tw(tw)>max_w) THEN max_w=rxy_tw(tw) k=twEND IFEND DOPRINT "(' 冬季平均气温绝对值最

19、大自相关系数 rxy_tw=',/,' 滞后时间长度 k=0!落后交叉相关系数DO tyw=1,N/2DO i=1,N-tyw rtyw(i)=(y(i)-avr_y)/sy)*(w(i+tyw)-avr_w)/sw) rxy_tyw(tyw)=rxy_tyw(tyw)+rtyw(i)END DO rxy_tyw(tyw)=rxy_tyw(tyw)/(N-tyw) rxy_tyw(tyw)=ABS(rxy_tyw(tyw)IF(rxy_tyw(tyw)>max_yw) THEN max_yw=rxy_tyw(tyw) k=tywEND IFEND DOrxy_tyw=&#

20、39;,/,' 滞后时间长度PRINT "(' 全年平均温度与冬季平均气温之间的落后交叉相关系数 k=',I2)",rxy_tyw(k),kk=',I2)",rxy_ty(k),kk=',I2)",rxy_tw(k),kEND2）、程序运行结果：实习三 分析中国夏季降水线性趋势的分布特征1. 资料介绍及要求：利用数据， 编写求 1982-2006 年中国 160 站各站夏季降水线性倾 向率，给出分布图，并进行简单分析。给出了阅读资料的 fortran 程序。数据在文件夹中单独给出。2. 实习结果1) .matlab

21、 程序%编写求 1982-2006 年中国 160 站各站夏季降水线性倾向率clear allclcfid=fopen( 'E:/' , 'rt' );tline=fgets(fid);data1=fscanf(fid,'%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f%f %f %f %f %f %f %f %f %f',28,160);data2=data1'fclose(fid);for i=1:160;j(i,1:25)=data2(i,4:28);n1=1982:1

22、:2006;pp(i,:)=polyfit(n1,j(i,1:25),1);endb=pp(:,1);jd=data2(:,3);wd=data2(:,2);jdc=75:135;wdc=18:.5:55;bz=griddata(jd,wd,b,jdc,wdc', 'cubic' );c=contour(jdc,wdc,bz)xlabel( ' 精度 ' );ylabel(' 纬度 ' );title('1982-2006 年中国 160 站各站夏季降水线性倾向率分布图)(2) 程序运行结果实习四 求给定数据的一元线性回归方程1、

23、资料介绍及要求利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流 指标之间的一元线性回归方程，并对回归方程进行检验。年份气温T环流指标195132195225:1953201954261955271956241957281958024195915196016196124196230196322196430196524196633196726196820196932197035答案：F=>F尸，回归方程显着2、实习结果(1)、MATLAB程序%实习四 求给定数据的一元线性回归方程ClimateData=xlsread('F: 气象统计方法实验数据 气象统计 实验四 数

24、据 .xls' )%从Excel 文件读取数据x=ClimateData(:,3);%提取 ClimateData的第三列，即环流指标y=ClimateData(:,2);%提取 ClimateData的第三列，即气温 Txdata=ones(size(x,1),1),x;%在原始数据x的左边加一列1，即模型包含常数项b,bint,r,rint,s=regress(y,xdata);%调用 regress函数作一元线性回归yhat=xdata*b;%计算 y 的估计值%定义元胞数组，以元胞数组形式显示系数的估计值和估计值得95%置信区间head1= ' 系数的估计值 '

25、 , ' 估计值的 95%置信下限 ' , ' 估计值的 95%置信上限 ' ; head1;num2cell(b,bint)%定义元胞数组，以元胞数组形式显示y 的真实值、 y 的估计值、残差和残差的95%置信区间head2= 'y 的真实值 ' , 'y 的估计值 ' , ' 残差 ' , ' 残差的 95%置信下限 ' , ' 残差的 95%置信上限 ' ; %同时显示 y 的真实值， y 的估计值、残差和残差的95%置信区间head2;num2cell(y,yhat,r,ri

26、nt)%定义元胞数组，以元胞数组形式显示判定系数、F统计量的观测值、检验的P值和误差方差的估计值head3='判定系数，F统计量的观测值,'检验的P值,'误差方差的估计值;head3;num2cell(s)(2)、程序运行结果ans =' 系数的估计值 ''估计值的95% 置信下限 ''估计值的95% 置信上限 ' ans ='y的真实值 ''y 的估计值 ''残差 ''残差的95%置信下限 ''残差的 95% 置信上限 0 ans =' 判定

27、系数 ''F 统计量的观测值''检验的 P 值 ''误差方差的估计值 (3)、结果分析从输出的结果看，常数项和回归系数的估计值分别为和，从而可以写出线性回归方程为 回归系数估计值的置信区间为 , 。对回归直线进行显着性检验，原假设和对立假设分别为检验P的值为X 10-4 <,可知显着性水品 a=下应拒绝原假设 H0,可认为y环流指数)与x (气温T的线 性关系是显着的。原始数据散点图与回归直线图>>plot(x,y, 'k.', 'Markersize' ,15)%原始数据散点图>>h

28、old on>>plot(x,yhat, 'linewidth' ,3)%回归直线图>>xlabel( '环流指标(x)') %标注x轴>>ylabel( '气温(y)') %标注 y 轴>>legend( ' 原始散点 ' , ' 回归直线 ' ) %加标注框实习五 对给定的海温数据进行 EOF 分析1、资料介绍给出海表温度距平数据资料，以及相应的数据描述文件，对其进行EOF 分析，资料的时空范围可以根据获知。数据在文件夹中单独给出，距平或者标准化距平处理后再进行

29、EOF给出了如何读取资料，为对距平或者标准化距平处理后的资料进行 EOF 分析2、要求实习报告中给出第一特征向量及其时间系数，并分析其时空特征3、实习结果：( 1)、 FORTRAN 源程序! prepare data for eof analysis! the program is to normalize sea surface temperature(SST)! mt: the length of time series;! mo: the month numbers;my: the year numbers;! sst: sea surface temperature data;!ss

30、t3: the work array;avf: the average of SST;! df: the variance of SST;program mainparameter(mo=12,my=43,nx=18,ny=12,mt=516)dimension avf(mo,nx,ny),df(mo,nx,ny)dimension sst(nx,ny,mt),sst3(nx,ny,mt) open(1,file='g:sstpx',form='unformatted',access='direct',recl=nx*ny) irec=1do i

31、t=1,mt read(1,rec=irec)(sst(i,j,it),i=1,nx),j=1,ny) irec=irec+1end do! averagedo j=1,nydo i=1,nxdo k=1,modo it=k,mt,12 avf(k,i,j)=avf(k,i,j)+sst(i,j,it) end doavf(k,i,j)=avf(k,i,j)/myend doend doend do! variancedo j=1,nydo i=1,nxdo k=1,modo it=k,mt,12df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j)*2 end

32、 dodf(k,i,j)=sqrt(df(k,i,j)/my)end doend doend do! standardizingdo j=1,nydo i=1,nxdo k=1,modo it=k,mt,12if(sst(i,j,it)=thensst3(i,j,it)=elsesst3(i,j,it)=(sst(i,j,it)-avf(k,i,j)/df(k,i,j)end ifend doend doend doend do! output file open(2,file='g:sstpx',form='unformatted',access='d

33、irect',recl=nx*ny) irec=1do it=1,mt write(2,rec=irec)(sst3(i,j,it),i=1,nx),j=1,ny) irec=irec+1end doclose(2)close(1)sstpx文件夹下的文件，特征end 分解后的时间系数写入文件中，空间场写入文件中，特征值和分析误差写入 向量写入文件。由 中的标准特征向量可得出一般特征值的前两个模态有效。用grads打开和，分别画岀海平面气温EOF分解后的空间场和时间序列。2)空间场和时间序列的 ctl 文件dset g:sstpxtitle Coads SSTA Eundefxdef

34、18 linear 120 10ydef 12 linear 5zdef 1 linear 1000 1tdef 2 linear 1jan1948 1month vars 1S 0 99 Coads SST anomaly interperated using endvars dset g:sstpx title Coads SSTA Tundefxdef 1linear 120 10ydef 1linear 5zdef 1linear 1000 1tdef 516 linear 1jan1948 1month vars 2a 0 99 time coefficient 1b 0 99 ti

35、me coefficient 2endvars3）运行结果第一模态空间场时间系数第二模态空间场时间系数第一特征向量文件给出了 EOF 分析的第一特征向量值的 216个值4）分析第一模态空间场时间系数此次试验 EOF 分析中的前两个特征向量最大限度地表征了海平面温度场的主要结构。第一特征向 量所描绘的第一经验正交函数的特征场(即第一模态)具有海表面气温 516 个样本的最相似的特征。若 其可以解释为 516 个月的标准化特征，它指示出海表温度变化的扰动。其对应的时间系数可以表示为第 一模态空间场的时间权重。 从第一模态的空间特征场可以看出，受到大尺度环流影响，整场的空间变化基本全为负值。而其值乘

36、以 时间权重后均变为负值。也就是大的时间系数乘以空间特征值对应海表温度的低值，而小的时间系数乘 以空间特征值则对应高值。海表温度的低值对应了气象上的拉尼娜年，而海表温度的高值对应了厄尔尼 诺年。厄尔尼诺现象泛指赤道附近的东部太平洋表层海水温度上升引起的气候异常现象，它是热带海洋洋 流与大气互作用的产物。其基本特征是太平洋沿岸的海面水温异常升高，海水水位上涨，并形成一股暖 流向南流动。它使原属冷水域的太平洋东部水域变成暖水域，结果引起海啸和暴风骤雨，造成一些地区 干旱，另一些地区又降雨过多的异常气候现象。所以，在空间特征场乘以时间系数后的高值表示厄尔尼 诺年。拉尼娜现象是指海洋中的赤道的中部和东

37、部太平洋，东西上万公里，南北跨度上千公里的范围 内，海洋温度比正常温度东部和中部海面温度偏低摄氏度，并持续半年(与现象正好相反)，东南信风 将表面被太阳晒热的海水吹向太平洋西部，致使西部比东部海平面增高将近 60 厘米，西部海水温度增 高，下降，潮湿空气积累形成和热带风暴，东部底层海水上翻，致使东太平洋海水变冷的现象。所以， 在空间特征场乘以时间系数后的低值表示拉尼娜年。实习三(附加) 计算给定数据的 11 年滑动平均和累积距 平1、资料介绍利用数据，编写11点滑动平均的程序，给出了阅读资料的fortran程序。 数据在文件夹中单独给出。2、要求实习报告中附出程序，并给出原数据和滑动后数据的图

38、形( 1张图) Matlab 程序 load 'g:' x=MA;year=1922:1:2006' year2=year(1+(ih-1)/2:length(x)-(ih-1)/2);ih=11;for i=1:length(x)-10 avr(i)=sum(x(i:i+10)/ih;endplot(year,x, 'b:' )hold onplot(year2,avr, 'r' )save ( '' , 'avr' , '-ascii' )3、实习结果 、 FORTRAN 程序滑动平均计

39、算值(已导入文件Exam_4 output ) 、原始数据及实验滑动平均拟合曲线图表（图表由matlab画岀） 、分析图中滑动平均值很好的描述了变量 x随时间的变化趋势。滑动平均值滤掉了较大的震荡，使趋 势更加明显。从图中看出， x的整体趋势被体现而震荡的极大值却在趋势线中没有显着表示。 所以，滑动平均在分析数据时可以更好的体现变化趋势但无法体现较大的异常值。在分析异常 时，不宜使用滑动平均。实习四（附加）求给定数据的多元线性回归方程1、说明X1-X4为四个预报因子，y为预报量；样本个数n=132、要求选取预报因子1、2、4,求预报量的标准化回归方程。i12345678910111213X17

40、111117113122111110X226295631525571315447406668X3615886917221842398X46052204733226442226341212y答案：标准化变量回归方程：? =0.5679人+0.4323x2 0.2613x4附加：利用上表资料，尝试编写逐步回归程序。3、实习结果：( 1) Matlab 源程序ClimateData=xlsread('F:-0 O -O伍气象统计方法实验数据气象统计实验四数据(附加题).xls' );%从Excel文件读取数据X=ClimateData(2:5,:);%提取 ClimateData 的

41、第二到 5行数据，即自变量观测矩阵 Xy=ClimateData(6,:);%提取 ClimateData的第六行，即预报量reglm(y,X)%调用 reglm 函数做 4重线性回归，显示回归分析的方差分析表和参数估计表function stats = reglm(y,X,model,varnames)% 多重线性回归分析或广义线性回归分析% reglm(y,X)，产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上 . 参%数X是自变量观测值矩阵，它是n行p列的矩阵.y是因变量观测值向量，它是n行1列的列向量.% stats = reglm(y,X) stats.%，还返回

42、一个包括了回归分析的所有诊断统计量的结构体变量%stats = reglm(y,X,model)，用可选的 model 参数来控制回归模型的类型 . model曰是一个字符串，%其可用的字符串如下%'linear'带有常数项的线性模型(默认情况)%'interaction'带有常数项、线性项和交叉项的模型%'quadratic'带有常数项、线性项、交叉项和平方项的模型%'purequadratic'带有常数项、线性项和平方项的模型%stats = reglm(y,X,model,varnames)，用可选的 varnames 参数

43、指定变量标签varnames%可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行%数或元胞数应与X的列数相同.默认情况下，用X1,X2,作为变量标签.%例 :%x = 215 250 180 250 180 215 180 215 250 215 215%'%y = '%reglm(y,x,'quadratic')%方差分析表% 方差来源 自由度 平方和 p值% 回归% 残差% 总计% 均方根误差 (Root MSE)% 因变量均值 (Dependent Mean)% %变量估计值p值%常数项%X1%X2% X1*X2% X

44、1*X1% X2*X2% Copyright 2009 - 2010 xiezhh.if nargin < 2error( ' 至少需要两个输入参数 ' ); end均方 F判定系数 (R-Square) 调整的判定系数 (Adj R-Sq)参数估计标准误t值p = size(X,2);% X 的列数，即变量个数if nargin < 3 | isempty(model)model = 'linear' ;% model 参数的默认值end% 生成变量标签 varnamesif nargin < 4 | isempty(varnames)var

45、name1 = strcat('X',num2str(1:p');varnames = makevarnames(varname1,model);else% 默认的变量标签if ischar(varnames) varname1 = cellstr(varnames);elseif iscell(varnames) varname1 = varnames(:);elseerror( 'varnames 必须是字符矩阵或字符串元胞数组 ' ); endif size(varname1,1) = perror('变量标签数与X的列数不一致);elsevarnames = makevarnames(varname1,model);% 指定的变量标签endendST = regstats(y,X,model);% 调用 regstats 函数进行线性回归分析，返回结构体变量 STf = ;% F 检验相关结果% 显示方差分析表fprintf(fprintf('n' );11fprintf('n' );fprintf('%s%7s%15s%15s%15s%12s'' , 'p 值 ' );fpr