大气统计作业.docx
- 文档编号:1159760
- 上传时间:2022-10-18
- 格式:DOCX
- 页数:20
- 大小:617.42KB
大气统计作业.docx
《大气统计作业.docx》由会员分享,可在线阅读,更多相关《大气统计作业.docx(20页珍藏版)》请在冰豆网上搜索。
大气统计作业
上机实践
一.相关回归分析
(1)计算北京夏季降水与气温的相关系数及秩相关系数,并进行显著性检验。
定义计算秩相关系数函数:
{functioncoeff=mySpearman(X,Y);
%本函数用于实现斯皮尔曼等级相关系数的计算操作
%
%输入:
%X:
输入的数值序列
%Y:
输入的数值序列
%
%输出:
%coeff:
两个输入数值序列X,Y的相关系数
iflength(X)~=length(Y)
error('两个数值数列的维数不相等');
return;
end
N=length(X);%得到序列的长度
Xrank=zeros(1,N);%存储X中各元素的排行
Yrank=zeros(1,N);%存储Y中各元素的排行
%计算Xrank中的各个值
fori=1:
N
cont1=1;%记录大于特定元素的元素个数
cont2=-1;%记录与特定元素相同的元素个数
forj=1:
N
ifX(i) cont1=cont1+1; elseifX(i)==X(j) cont2=cont2+1; end end Xrank(i)=cont1+mean([0: cont2]); end %计算Yrank中的各个值 fori=1: N cont1=1;%记录大于特定元素的元素个数 cont2=-1;%记录与特定元素相同的元素个数 forj=1: N ifY(i) cont1=cont1+1; elseifY(i)==Y(j) cont2=cont2+1; end end Yrank(i)=cont1+mean([0: cont2]); end %利用差分等级(或排行)序列计算斯皮尔曼等级相关系数 fenzi=6*sum((Xrank-Yrank).^2); fenmu=N*(N^2-1); coeff=1-fenzi/fenmu; end%函数mySpearman结束} 计算相关系数: A=load('RBJsum.txt'); B=load('TBJsum.txt'); A1=A(: 2); B1=B(: 2); %相关系数 R=corrcoef(A1,B1); R= 1.0000-0.4332 -0.43321.0000 %秩相关系数 X=A1; Y=B1; coeff=mySpearman(X,Y); coeff coeff= -0.4693 (2)显著性检验: 相关系数自由度n=52-2 查表可知,rc=0.273<0.4332 计算得到的相关系数大于0.273,则在显著水平0.05是显著的。 秩相关系数: 同理,对秩相关系数检验可得: rc=0.273<0.4693,因此秩相关系数在显著水平0.05上是显著的,通过检验。 (2)将降水量作为预报因子,气温作为预报量,试给出回归方程,并说明降水每增加100mm,气温大致下降多少ºC? legend('降水量-气温散点图'); scatter(x,y) xlabel('降水量');ylabel('气温'); title('降水量-气温散点图'); %建立回归方程: stepwise(x,y); b=25.995 -0.0019776 x1=[ones(52,1),x]; y1=x1*b; plot(x,y,'*',x1,y1); xlabel('降水量');ylabel('气温'); title('拟合曲线及散点图'); 那么降水每增加100mm,气温大致下降0.19776ºC 二、EOF分析 (1)进行EOF分析,给出前十个模态的方差解释率,并根据North准则(绘出ErrorBar图),说明需要选择前几个模态进行分析? fid=fopen('ssta.dat','r'); a=fread(fid,’float’);%a=reshape(a,56*11,336);b=a’; b=zeros(336,616); c=zeros(1,616); j=0; forn=1: 336; c=a(j+1: j+616); b(n,: )=reshape(c,1,616);j=j+616;%b1=fread(fid,[616336]);b=b’; end b=b'; C=b*b'/616; [EOF,E]=eig(C); PC=EOF'*b; E=fliplr(flipud(E)); lambda=diag(E); EOF=fliplr(EOF); PC=flipud(PC); sum_d=sum(lambda); count=0; fori=1: 336; count=count+lambda(i); G1(i)=count/sum_d;%累计方差贡献率 end fori=1: 336; G2(i)=lambda(i)/sum_d;%方差贡献率 end 方差贡献率 累计方差贡献率 第一模态 0.716516015743107 0.716516015743107 第二模态 0.147473083981738 0.863989099724845 第三模态 0.050392162180530 0.914381261905375 第四模态 .024********* 0.938899939827844 第五模态 0.014502739245485 0.953402679073329 第六模态 0.012702830971315 0.966105510044644 第七模态 0.007148842058117 0.973254352102761 第八模态 0.006019078625006 0.979273430727767 第九模态 0.004340892184153 0.983614322911921 第十模态 0.003048921846347 0.986663244758267 %North准则 f=lambda(1: 10); errorf=f*sqrt(2/(336-2)); number=1: 1: 10; errorbar(number,f,errorf); xlabel('number');ylabel('eigenvalue'); title('northerrorbar'); 只有对第一、二、三、四、五模态进行分析。 因为第五模态后误差重合了,无法区分开来。 (2)给出需要分析的前几个模态的空间分布和时间序列。 x=1: 336; plot(x,PC(1,: ),'r'); xlabel('1979-2006年336个月'); ylabel('INDEX'); title('赤道中东太平洋第一模态1979-2006年的SST时间序列'); gridon; lon=[160: 2: 270]; lat=[-10: 2: 10]; m_proj('EquidistantCylindrical','lat',[-1010],'lon',[160270]); m_contourf(lon,lat,rot90(reshape(EOF(: 1),56,11)',2),'linestyle','none'); colorbar; m_grid('box','fancy','tickdir','in'); xlabel('longitude','fontsize',15,'fontname','comicsansms'); ylabel('latitude','fontsize',15,'fontname','comicsansms'); title('北太平洋第1模态填色图','fontsize',15,'fontname','幼圆');%MATLAB安装m_map绘图工具才能调用以上程序,如果没有,直接调用eof1=reshape(eof(: n),56,11); eof1=eof1';contour(eof1);n表示不同的模态 三、聚类分析 x= 1.0e+003* 0.01621.42902.0000-0.00820.0062 0.01570.97002.2090-0.02060.0019 0.01631.26002.0850-0.01730.0028 0.01721.44201.7260-0.00950.0046 0.01881.87401.7090-0.00490.0080 0.01791.69801.8480-0.00450.0075 0.01630.97601.2390-0.00460.0056 x1=zscore(x); y1=pdist(x2); y1=pdist(x1); y1=pdist(x1,'mahalanobis'); z1=linkage(y1); c1=cophenet(z1,y1); T=cluster(z1,6); H=dendrogram(z1); c1 c1= 0.5540 T T= 3 1 2 2 4 5 6 四、功率谱分析 (1)太阳黑子数的离散功率谱,绘出谱图。 sunspot=importdata('sunspot.dat'); year=sunspot(: 1); wolfer=sunspot(: 2); figure plot(year,wolfer) xlabel('Years');ylabel('SunspotData');title('SunspotData') figure plot(year(1: 50),wolfer(1: 50),'b.-'); xlabel('Years');ylabel('SunspotData');title('Atthefirst50years') figure n=length(Y); power=abs(Y(1: n/2)).^2; nyquist=1/2;%取采样频率为0.5; freq=(1: n/2)/(n/2)*nyquist; stem(freq,power); xlabel('cycles/year');title('Periodogram') figure period=1./freq; stem(period,power); axis([05002e+7]); ylabel('Power'); xlabel('Period(Years/cycle)'); holdon; index=find(power==max(power)); mainPeriodStr=num2str(period(index)); plot(period(index),power(index),'r.','MarkerSize',25); text(period(index)+2,power(index),['Period=',mainPeri
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大气 统计 作业