大氣統(tǒng)計(jì)作業(yè)

上機(jī)實(shí)踐相關(guān)回歸分析（1）計(jì)算北京夏季降水與氣溫的相關(guān)系數(shù)及秩相關(guān)系數(shù)，并進(jìn)行顯著性檢驗(yàn)。定義計(jì)算秩相關(guān)系數(shù)函數(shù)：{functioncoeff=mySpearman(X,Y);%本函數(shù)用于實(shí)現(xiàn)斯皮爾曼等級相關(guān)系數(shù)的計(jì)算操作% %輸入： %X：輸入的數(shù)值序列 %Y：輸入的數(shù)值序列 % %輸出： %coeff：兩個(gè)輸入數(shù)值序列X，Y的相關(guān)系數(shù) iflength(X)~=length(Y) error('兩個(gè)數(shù)值數(shù)列的維數(shù)不相等'); return; end N=length(X);%得到序列的長度 Xrank=zeros(1,N);%存儲(chǔ)X中各元素的排行 Yrank=zeros(1,N);%存儲(chǔ)Y中各元素的排行 %計(jì)算Xrank中的各個(gè)值 fori=1:N cont1=1;%記錄大于特定元素的元素個(gè)數(shù) cont2=-1;%記錄與特定元素相同的元素個(gè)數(shù) forj=1:N ifX(i)<X(j)cont1=cont1+1; elseifX(i)==X(j) cont2=cont2+1; end end Xrank(i)=cont1+mean([0:cont2]); end %計(jì)算Yrank中的各個(gè)值 fori=1:N cont1=1;%記錄大于特定元素的元素個(gè)數(shù) cont2=-1;%記錄與特定元素相同的元素個(gè)數(shù)forj=1:N ifY(i)<Y(j) cont1=cont1+1; elseifY(i)==Y(j) cont2=cont2+1; endendYrank(i)=cont1+mean([0:cont2]); end %利用差分等級(或排行)序列計(jì)算斯皮爾曼等級相關(guān)系數(shù) fenzi=6*sum((Xrank-Yrank).^2); fenmu=N*(N^2-1); coeff=1-fenzi/fenmu; end%函數(shù)mySpearman結(jié)束}計(jì)算相關(guān)系數(shù)：A=load('RBJsum.txt');B=load('TBJsum.txt');A1=A(:,2);B1=B(:,2);%相關(guān)系數(shù)R=corrcoef(A1,B1);R=1.0000-0.4332-0.43321.0000%秩相關(guān)系數(shù)X=A1;Y=B1;coeff=mySpearman(X,Y);coeffcoeff=-0.4693(2)顯著性檢驗(yàn)：相關(guān)系數(shù)自由度n=52-2查表可知，rc=0.273<0.4332計(jì)算得到的相關(guān)系數(shù)大于0.273，則在顯著水平0.05是顯著的。秩相關(guān)系數(shù)：同理，對秩相關(guān)系數(shù)檢驗(yàn)可得：rc=0.273<0.4693，因此秩相關(guān)系數(shù)在顯著水平0.05上是顯著的，通過檢驗(yàn)。（2）將降水量作為預(yù)報(bào)因子，氣溫作為預(yù)報(bào)量，試給出回歸方程，并說明降水每增加100mm，氣溫大致下降多少oC？legend('降水量-氣溫散點(diǎn)圖');scatter(x,y)xlabel('降水量');ylabel('氣溫');title('降水量-氣溫散點(diǎn)圖');%建立回歸方程：stepwise(x,y)；b=25.995-0.0019776x1=[ones(52,1),x];y1=x1*b;plot(x,y,'*',x1,y1);xlabel('降水量');ylabel('氣溫');title('擬合曲線及散點(diǎn)圖');那么降水每增加100mm，氣溫大致下降0.19776oC二、EOF分析（1）進(jìn)行EOF分析，給出前十個(gè)模態(tài)的方差解釋率，并根據(jù)North準(zhǔn)則（繪出ErrorBar圖），說明需要選擇前幾個(gè)模態(tài)進(jìn)行分析？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’;endb=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;%累計(jì)方差貢獻(xiàn)率endfori=1:336;G2(i)=lambda(i)/sum_d;%方差貢獻(xiàn)率end方差貢獻(xiàn)率累計(jì)方差貢獻(xiàn)率第一模態(tài)0.7165160157431070.716516015743107第二模態(tài)0.1474730839817380.863989099724845第三模態(tài)0.0503921621805300.914381261905375第四模態(tài)0.0245186779224680.938899939827844第五模態(tài)0.0145027392454850.953402679073329第六模態(tài)0.0127028309713150.966105510044644第七模態(tài)0.0071488420581170.973254352102761第八模態(tài)0.0060190786250060.979273430727767第九模態(tài)0.0043408921841530.983614322911921第十模態(tài)0.0030489218463470.986663244758267%North準(zhǔn)則f=lambda(1:10);errorf=f*sqrt(2/(336-2));number=1:1:10;errorbar(number,f,errorf);xlabel('number');ylabel('eigenvalue');title('northerrorbar');只有對第一、二、三、四、五模態(tài)進(jìn)行分析。因?yàn)榈谖迥B(tài)后誤差重合了，無法區(qū)分開來。（2）給出需要分析的前幾個(gè)模態(tài)的空間分布和時(shí)間序列。x=1:336;plot(x,PC(1,:),'r');xlabel('1979-2006年336個(gè)月');ylabel('INDEX');title('赤道中東太平洋第一模態(tài)1979-2006年的SST時(shí)間序列');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模態(tài)填色圖','fontsize',15,'fontname','幼圓');%MATLAB安裝m_map繪圖工具才能調(diào)用以上程序，如果沒有，直接調(diào)用eof1=reshape(eof(:,n),56,11);eof1=eof1';contour(eof1);n表示不同的模態(tài)三、聚類分析x=1.0e+003*0.01621.42902.0000-0.00820.00620.01570.97002.2090-0.02060.00190.01631.26002.0850-0.01730.00280.01721.44201.7260-0.00950.00460.01881.87401.7090-0.00490.00800.01791.69801.8480-0.00450.00750.01630.97601.2390-0.00460.0056x1=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);c1c1=0.5540TT=3122456四、功率譜分析（1）太陽黑子數(shù)的離散功率譜，繪出譜圖。sunspot=importdata('sunspot.dat');year=sunspot(:,1);wolfer=sunspot(:,2);figureplot(year,wolfer)xlabel('Years');ylabel('SunspotData');title('SunspotData')figureplot(year(1:50),wolfer(1:50),'b.-');xlabel('Years');ylabel('SunspotData');title('Atthefirst50years')figuren=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')figureperiod=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=',mainPeriodStr]);holdoff主周期11.0385年顯著性檢驗(yàn)：%不確定對否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時(shí)，F(xiàn)α＝3.59，|F|＞Fα通過檢驗(yàn)（2）利用落后自相關(guān)方法計(jì)算太陽黑子數(shù)的連續(xù)功率譜密度，繪出譜圖并給出95%信度檢驗(yàn)曲線。（參考程序spectrum.m）[num,Period,Frequency,Density,CL95]=spectrum007(x2,145);plot(Density,Period,'r');plot(Period,Density,'b');holdonindex=find(Density==max(Density));mainPeriodStr=num2str(Period(index));plot(Period(index),Density(index),'r','MarkerSize',25);text(Period(index)+2,Density(index),['Period=',mainPeriodStr]);holdoffxlabel('Period');ylabel('Density');holdonplot(Period,CL95,'r');holdoff五、小波分析（1）請采用Morlet小波對太陽黑子活動(dòng)進(jìn)行小波分析，繪出小波能量譜圖。x=importdata('sunspot.dat');stem(x(:,1),x(:,2),'.')xlabel('Year/a');ylabel('sunspot');title('太陽黑子數(shù)年變化')>>x1=x(:,1);>>x2=x(:,2);>>coefs=cwt(x2,1:144,'morl','plot');%將尺度參數(shù)設(shè)為[1:144]小波系數(shù)等值線圖:figure;[c,h]=contour(x1,[1:144],coefs,'-');colorbarxlabel('year/a');ylabel('period/a');能量譜圖fori=1:144coefs1(i,:)=coefs(i,:)./i^2;end;%小波能量contour(coefs1(2:20,:));colorbar;xlabel('year/a','fontsize',15,'f