当前位置:首页 >> 能源/化工 >>

大气统计作业


上机实践
一. 相关回归分析 (1)计算北京夏季降水与气温的相关系数及秩相关系数,并进行显著性检验。
定义计算秩相关系数函数:

{function coeff = mySpearman(X,Y); % 本函数用于实现斯皮尔曼等级相关系数的计算操作 % % 输入: % X:输入的数值序列 % Y:输入的数值序列 % % 输出: % coef

f:两个输入数值序列 X,Y 的相关系数

if length(X)~=length(Y) error('两个数值数列的维数不相等'); return; end N = length(X); %得到序列的长度 Xrank = zeros(1 , N); %存储 X 中各元素的排行 Yrank = zeros(1 , N); %存储 Y 中各元素的排行 %计算 Xrank 中的各个值 for i = 1 : N cont1 = 1; %记录大于特定元素的元素个数 cont2 = -1; %记录与特定元素相同的元素个数 for j = 1 : N if X(i) < X(j) cont1 = cont1 + 1; elseif X(i) == X(j) cont2 = cont2 + 1; end end Xrank(i) = cont1 + mean([0 : cont2]); end %计算 Yrank 中的各个值 for i = 1 : N cont1 = 1; %记录大于特定元素的元素个数

cont2 = -1; %记录与特定元素相同的元素个数 for j = 1 : N if Y(i) < Y(j) cont1 = cont1 + 1; elseif Y(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.4332 1.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 准则(绘出 Error Bar 图) ,说明需要选择前几个模态进行分析? 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; for n=1:336; c=a(j+1:j+616); b(n,:)=reshape(c,1,616);j=j+616;%b1=fread(fid,[616 336]);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; for i=1:336; count=count+lambda(i); G1(i)=count/sum_d;%累计方差贡献率 end for i=1:336; G2(i)=lambda(i)/sum_d;%方差贡献率 end

方差贡献率 第一模态 第二模态 第三模态 第四模态 第五模态 第六模态 第七模态 第八模态 第九模态 第十模态 0.716516015743107 0.147473083981738 0.050392162180530 0.024518677922468 0.014502739245485 0.012702830971315 0.007148842058117 0.006019078625006 0.004340892184153 0.003048921846347

累计方差贡献率 0.716516015743107 0.863989099724845 0.914381261905375 0.938899939827844 0.953402679073329 0.966105510044644 0.973254352102761 0.979273430727767 0.983614322911921 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('north errorbar');

north errorbar 140

120

100

eigenvalue

80

60

40

20

0

0

2

4

6 number

8

10

12

north errorbar 70

60

50

eigenvalue

40

30

20

10

0

3

4

5 number

6

7

8

只有对第一、二、三、四、五模态进行分析。因为第五模态后误差重合了,无法区分开来。 (2)给出需要分析的前几个模态的空间分布和时间序列。 x=1:336; plot(x,PC(1,:),'r'); xlabel('1979-2006 年 336 个月'); ylabel('INDEX'); title('赤道中东太平洋第一模态 1979-2006 年的 SST 时间序列'); grid on;

赤 道 中 东 太 平 洋 第 一 模 态 1979-2006年 的 SST时 间 序 列 50 40 30 20
INDEX

10 0 -10 -20 -30 -40

0

50

100

150 200 1979-2006年 336个 月

250

300

350

lon=[160:2:270]; lat=[-10:2:10]; m_proj('Equidistant Cylindrical','lat',[-10 10],'lon',[160 270]); 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','comic sans ms'); ylabel('latitude','fontsize',15,'fontname','comic sans ms'); title('北太平洋第 1 模态填色图','fontsize',15,'fontname','幼圆');%MATLAB 安装 m_map 绘图工 具才能调用以上程序,如果没有,直接调用 eof1=reshape(eof(:,n),56,11); eof1=eof1';contour(eof1);n 表示不同的模态

北太平洋第1模态填色图
latitude
8oN 4oN 0o 4oS o 0.06 0.04 0.02 0 180oW 160oW

8 S 160oE

longitude

140oW

120oW

100oW

赤 道 中 东 太 平 洋 第 二 模 态 1979-2006年 的 SST时 间 序 列 20

15

10

INDEX

5

0

-5

-10

-15

0

50

100

150 200 1979-2006年 336个 月

250

300

350

北太平洋第2模态填色图
latitude
8oN 4oN 0o 4oS o 0.05 0 -0.05 180oW 160oW

8 S 160oE

longitude

140oW

120oW

100oW

赤 道 中 东 太 平 洋 第 三 模 态 1979-2006年 的 SST时 间 序 列 10

5

0
INDEX

-5

-10

-15

0

50

100

150 200 1979-2006年 336个 月

250

300

350

北太平洋第3模态填色图
latitude
8oN 4oN 0o 4oS o 0.05 0 -0.05 -0.1 180oW 160oW

8 S 160oE

longitude

140oW

120oW

100oW

赤 道 中 东 太 平 洋 第 四 模 态 1979-2006年 的 SST时 间 序 列 8 6 4 2
INDEX

0 -2 -4 -6 -8 -10

0

50

100

150 200 1979-2006年 336个 月

250

300

350

北太平洋第4模态填色图
latitude
8oN 4oN 0o 4oS o 0.05 0 -0.05 180 W
o

8 S 160oE

160 W

o

longitude

140 W

o

120 W

o

100 W

o

赤 道 中 东 太 平 洋 第 五 模 态 1979-2006年 的 SST时 间 序 列 6

4

2
INDEX

0

-2

-4

-6

0

50

100

150 200 1979-2006年 336个 月

250

300

350

北太平洋第5模态填色图
latitude
8oN 4oN 0o 4oS o 0.1 0.05 0 -0.05 180oW 160oW

8 S 160oE

longitude

140oW

120oW

100oW

三、聚类分析
x= 1.0e+003 * 0.0162 1.4290 2.0000 -0.0082 0.0062 0.0157 0.9700 2.2090 -0.0206 0.0019 0.0163 1.2600 2.0850 -0.0173 0.0028 0.0172 1.4420 1.7260 -0.0095 0.0046 0.0188 1.8740 1.7090 -0.0049 0.0080 0.0179 1.6980 1.8480 -0.0045 0.0075 0.0163 0.9760 1.2390 -0.0046 0.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('Sunspot Data');title('Sunspot Data')
Sunspot Data 200 180 160 140
Sunspot Data

120 100 80 60 40 20 0 1700

1750

1800

1850 Years

1900

1950

2000

figure plot(year(1:50),wolfer(1:50),'b.-');

xlabel('Years');ylabel('Sunspot Data');title('At the first 50 years')
At the first 50 years 140

120

100
Sunspot Data

80

60

40

20

0 1700

1705

1710

1715

1720

1725 Years

1730

1735

1740

1745

1750

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')
2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 x 10
7

Periodogram

0

0.05

0.1

0.15

0.2 0.25 0.3 cycles/year

0.35

0.4

0.45

0.5

figure period=1./freq; stem(period,power); axis([0 50 0 2e+7]);

ylabel('Power'); xlabel('Period(Years/cycle)'); hold on; index=find(power==max(power)); mainPeriodStr=num2str(period(index)); plot(period(index),power(index),'r.','MarkerSize',25); text(period(index)+2,power(index),['Period=',mainPeriodStr]); hold off
2 1.8 1.6 1.4 1.2
Power

x 10

7

Period=11.0385

1 0.8 0.6 0.4 0.2 0

0

5

10

15

20 25 30 Period(Years/cycle)

35

40

45

50

主周期 11.0385 年 显著性检验:%不确定对否 y1=var(wolfer); y2=max(power); F=(0.5*y2)/2/((y1-0.5*y2)/(143-2-1)); abs(F)= ans = 70.0119 显著水平为 0.05 时,Fα =3.59,|F|>Fα 通过检验

(2)利用落后自相关方法计算太阳黑子数的连续功率谱密度,绘出谱图并给
出 95%信度检验曲线。 (参考程序 spectrum.m)
[num,Period,Frequency,Density,CL95]=spectrum007(x2,145); plot(Density,Period,'r'); plot(Period,Density,'b');

hold on index=find(Density==max(Density)); mainPeriodStr=num2str(Period(index)); plot(Period(index),Density(index),'r','MarkerSize',25); text(Period(index)+2,Density(index),['Period=',mainPeriodStr]); hold off xlabel('Period');ylabel('Density');

0.16 0.14 0.12 0.1
Density

Period=11.1538

0.08 0.06 0.04 0.02 0 -0.02

0

50

100

150 Period

200

250

300

hold on plot(Period,CL95,'r'); hold off

0.2 95% 置 信 度 曲 线 0.15 Period=11.1538

0.1
Density

0.05

0

-0.05

0

50

100

150 Period

200

250

300

五、小波分析 (1)请采用 Morlet 小波对太阳黑子活动进行小波分析,绘出小波能量谱图。
太阳黑子数年变化 200 180 160 140 120
sunspot

100 80 60 40 20 0 1700

1750

1800

1850 Year/a

1900

1950

2000

x=importdata('sunspot.dat'); stem(x(:,1),x(:,2),'.') xlabel('Year/a');ylabel('sunspot'); title('太阳黑子数年变化')

>> x1=x(:,1);

>> x2=x(:,2); >> coefs=cwt(x2,1:144,'morl','plot');%将尺度参数设为[1:144]
Absolute of Values of Ca,b for a = 1 2 3 4 5 ... 137 129 121 113 105 97 89 81 73 65 57 49 41 33 25 17 9 1 50 100 150 Time (or Space) b 200 250

小波系数等值线图:
figure;[c,h]=contour(x1,[1:144],coefs,'-'); colorbar xlabel('year/a');ylabel('period/a');
250 200 120 150 100 50
period/a

Scales a

140

100

80 0 60 -50 -100 -150 20 -200 1750 1800 1850 year/a 1900 1950 -250

40

1700

能量谱图
for i=1:144 coefs1(i,:)=coefs(i,:)./i^2;

end ;%小波能量 contour(coefs1(2:20,:)); colorbar; xlabel('year/a','fontsize',15,'fontname','Times New Roman'); ylabel('period/a','fontsize',15,'fontname','Times New Roman');title('小波能量图');
小波能量图 3 18 16 14 2

1

period/a

12 0 10 -1 8 6 4 2 50 100 150 200 250 -4 -2

-3

year/a

(2)根据小波方差判断多长周期的振动最强。
y=sum(abs(coefs).^2,2);%计算小波方差 x3=[1:144]; plot(x3,y); xlabel('period scales/a','fontsize',15); ylabel('Var','fontsize',15); title('小波方差图','fontsize',15);

6

x 10

6

小波方差图

5

4

Var

3

2

1

0

0

50

100

150

period scales/a
%取前 50 个数据 y=sum(abs(coefs).^2,2);%计算小波方差 x3=[1:144]; plot(x3(1:50),y(1:50)); xlabel('period scales/a','fontsize',15); ylabel('Var','fontsize',15); title('小波方差图','fontsize',15);

hold on; index=find(y==max(y(1:50))); mainx3Str=num2str(x3(index)); plot(x3(index),y(index),'r.','MarkerSize',25); text(x3(index)+2,y(index),['Period=',mainx3Str]); hold off

3.5

x 10

6

小波方差图
Period=9

3

2.5

2

Var
1.5 1 0.5 0 0

5

10

15

20

25

30

35

40

45

50

period scales/a

所用 morl 小波的周期为 2*pi/5, T ? 9 ?
图分析可知,11.30 年时,振动最强。

2? ? 11.30 5


相关文章:
大气统计作业
大气统计作业_能源/化工_工程科技_专业资料。中科大大气统计结课作业上机实践一. 相关回归分析 (1)计算北京夏季降水与气温的相关系数及秩相关系数,并进行显著性检验...
统计作业模板201301
1 《中医药统计学与软件应用[070103]》 论文评价作业汇编授课教师 魏高文 班级...测试期间大气平均温度 9"C。大气平均 压力 825 hPa,地表 5 cm 深度土层平均...
统计学作业
统计作业_数学_自然科学_专业资料。南医大 13 级护理本科班《医学统计学》书面...乙两种方法测定大气中 SO2 日平均浓度(μ g/m ) ,测定值如下表 所示,问...
统计法规概论》作业(一)
统计法规概论》平时作业(一) 一、单项选择题 1、我国历史上第一部现代意义上的统计法律是( A )。 A、 《中华民国统计法》 B、 《中华人民共和国统计法》...
5统计学作业——我国城市空气质量的状况分析
5统计作业——我国城市空气质量的状况分析_经济学_高等教育_教育专区。天津商业...2、工业污染源:随着经济的迅猛发展,各类工厂迅速增加,而大气污染的来源也迅速...
16秋北航《统计学》在线作业二
16秋北航《统计学》在线作业二_工学_高等教育_教育专区。北航《统计学》在线作业二 一、单选题(共 12 道试题,共 48 分。 ) 1. 编制数量指标综合指数的一般...
16秋北航《统计学》在线作业三
16秋北航《统计学》在线作业三_工学_高等教育_教育专区。北航《统计学》在线作业三 一、单选题(共 12 道试题,共 48 分。 ) 1. 权数对算术平均数的影响...
数理统计判别分析作业
数理统计判别分析作业_工程科技_专业资料。数理统计判别分析作业判别分析大气环境质量评价是环境质量评价的一项重要内容。 对空气环境质量的充分 认识对我国社会的可持续...
生物统计学作业
生物统计作业_工学_高等教育_教育专区。简单的生物统计作业,软件SPSS ...表 1 2004 至 2013 年孝感市大气主要污染物年均值统计数据 二氧化硫年均值 ...
大气降雨统计表
大气降水及影响 19页 1财富值如要投诉违规内容,请到百度文库投诉中心;如要提出功能问题或意见建议,请点击此处进行反馈。 大气降雨统计表 隐藏>> 付村煤业公司大气...
更多相关标签:
北航数理统计大作业 | 数理统计大作业 | 统计学作业答案 | 统计学原理作业1 | 教育统计学作业答案 | 北师大概率统计作业 | 华工统计学原理作业 | 概率论与数理统计作业 |