大气统计作业Word文件下载.docx
《大气统计作业Word文件下载.docx》由会员分享，可在线阅读，更多相关《大气统计作业Word文件下载.docx（20页珍藏版）》请在冰豆网上搜索。
fori=1:
N
cont1=1;
%记录大于特定元素的元素个数
cont2=-1;
%记录与特定元素相同的元素个数
forj=1:
ifX（i）<
X（j）
cont1=cont1+1;
elseifX（i）==X（j）
cont2=cont2+1;
Xrank（i）=cont1+mean（[0:
cont2]）;
%计算Yrank中的各个值
ifY（i）<
Y（j）
elseifY（i）==Y（j）
Yrank（i）=cont1+mean（[0:
%利用差分等级（或排行）序列计算斯皮尔曼等级相关系数
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（:
%相关系数
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）;
拟合曲线及散点图'
那么降水每增加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:
count=count+lambda（i）;
G1（i）=count/sum_d;
%累计方差贡献率
G2（i）=lambda（i）/sum_d;
%方差贡献率
方差贡献率
累计方差贡献率
第一模态
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）;
number'
eigenvalue'
northerrorbar'
只有对第一、二、三、四、五模态进行分析。
因为第五模态后误差重合了，无法区分开来。
（2）给出需要分析的前几个模态的空间分布和时间序列。
x=1:
plot（x,PC（1,:
）,'
1979-2006年336个月'
INDEX'
赤道中东太平洋第一模态1979-2006年的SST时间序列'
gridon;
lon=[160:
2:
270];
lat=[-10:
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'
longitude'
fontsize'
15,'
fontname'
comicsansms'
latitude'
北太平洋第1模态填色图'
幼圆'
%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
4
5
6
四、功率谱分析
（1）太阳黑子数的离散功率谱，绘出谱图。
sunspot=importdata（'
sunspot.dat'
year=sunspot（:
1）;
wolfer=sunspot（:
figure
plot（year,wolfer）
xlabel（'
Years'
SunspotData'
）
plot（year（1:
50）,wolfer（1:
50）,'
b.-'
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）;
cycles/year'
Periodogram'
period=1./freq;
stem（period,power）;
axis（[05002e+7]）;
Power'
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
升级会员 