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

1、上機實踐一相關(guān)回歸分析(1)計算北京夏季降水與氣溫的相關(guān)系數(shù)及秩相關(guān)系數(shù)，并進行顯著性檢驗定義計算秩相關(guān)系數(shù)函數(shù)：functioncoeff=mySpearman(X,Y);%本函數(shù)用于實現(xiàn)斯皮爾曼等級相關(guān)系數(shù)的計算操作%輸入：%X：輸入的數(shù)值序列%Y：輸入的數(shù)值序列%輸出：%coeff:兩個輸入數(shù)值序列X,Y的相關(guān)系數(shù)iflength(X)=length(Y)error('兩個數(shù)值數(shù)列的維數(shù)不相等')；return;endN=length(X)；%得到序列的長度Xrank=zeros(1,N)；%存儲X中各元素的排行Yrank=zeros(1,N)；%存儲Y中各元素的排行%計

2、算Xrank中的各個值fori=1:Ncont1=1；%記錄大于特定元素的元素個數(shù)cont2=-1；%記錄與特定元素相同的元素個數(shù)forj=1:NifX(i)<X(j)cont1=cont1+1；elseifX(i)=X(j)cont2=cont2+1；endendXrank(i)=cont1+mean(0:cont2)；end%計算Yrank中的各個值fori=1:Ncont1=1；%記錄大于特定元素的元素個數(shù)cont2=-1;%記錄與特定元素相同的元素個數(shù)forj=1:NifY(i)<Y(j)cont1=cont1+1;elseifY(i)=Y(j)cont2=cont2+1;

3、endendYrank(i)=cont1+mean(0:cont2);end%利用差分等級(或排行)序列計算斯皮爾曼等級相關(guān)系數(shù)fenzi=6*sum(Xrank-Yrank).人2);fenmu=N*(NA2-1);coeff=1-fenzi/fenmu;end%函數(shù)mySpearman結(jié)束計算相關(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;c

4、oeff=mySpearman(X,Y);coeffcoeff=-0.4693(2)顯著性檢驗：相關(guān)系數(shù)自由度n=52-2查表可知，rc=0.273v0.4332計算得到的相關(guān)系數(shù)大于0.273,則在顯著水平0.05是顯著的。相關(guān)系集總譽性槍驗表自由度(n-2)顯著性朮口自由度顯著性水平自由度(n-2)顯善性水平0.0c0.010.030.010.Cc0.0110.9971,160.46S0.59350.32o0.41S20.950.991'0.4560.；炕400.3040.39330.9780.959L80.44561閔0.2S80.372=0.S1：0.917190.4330.“

5、49d00.2730.3545C.7340.8742C0.4230.537600.2S0.32580.7070.834210.526700.2320.30270.B660.798220.4040.515SOD.2L70.283SC.6320.765230.3%0.503900.2050.2679C.G020.735240.33S6=961.0061950.254.00.3760.70S250.3810.1.250.228110.0.684260.370.781300.1590.208二0.3C20.661270.3670.472000.380.SI13C.5140.641280.3610,46

6、33000,1130.143140.4970.623290.麺0,4564000,0980.12815C.4820,60630C.S490.44910000,0620.081秩相關(guān)系數(shù)：同理，對秩相關(guān)系數(shù)檢驗可得：rc=0.273v0.4693，因此秩相關(guān)系數(shù)在顯著水平0.05上是顯著的，通過檢驗。(2)將降水量作為預(yù)報因子，氣溫作為預(yù)報量，試給出回歸方程，并說明降水每增加100mm，氣溫大致下降多少°C?legend('降水量-氣溫散點圖');scatter(x,y)xlabel('降水量');ylabel(氣溫');title('降

7、水量-氣溫散點圖');%建立回歸方程：stepwise(x,y)；b=25.995-0.0019776xl=ones(52,l),x;y1=x1*b;plot(x,y,'*',x1,y1);xlabel('降水量');ylabel(氣溫');title('擬合曲線及散點圖')；那么降水每增加100mm,氣溫大致下降0.19776°C二、EOF分析(1)進行EOF分析，給出前十個模態(tài)的方差解釋率，并根據(jù)North準(zhǔn)則(繪出ErrorBar圖)，說明需要選擇前幾個模態(tài)進行分析？fid=fopen('ssta.dat&

8、#39;,'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

9、(EOF);PC=flipud(PC);sum_d=sum(lambda);count=0;fori=1:336;count=count+lambda(i);Gl(i)=count/sum_d;%累計方差貢獻率endfori=l:336;G2(i)=lambda(i)/sum_d;%方差貢獻率end方差貝獻率累計方差貢獻率第一模態(tài)0.7165160157431070.716516015743107第一模態(tài)0.1474730839817380.863989099724845第三模態(tài)0.0503921621805300.914381261905375第四模態(tài)0.0245186779224680.9

10、38899939827844第五模態(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(l:l0);errorf=f*sqrt(2/(336-2);number=l:l:

11、l0;errorbar(number,f,errorf);xlabel('number');ylabel('eigenvalue');title('northerrorbar');eulavnegie°。northerrorbar1401201008060402024810126numbernumbereulavnegi只有對第一、二、三、四、五模態(tài)進行分析。因為第五模態(tài)后誤差重合了，無法區(qū)分開來(2)給出需要分析的前幾個模態(tài)的空間分布和時間序列。x=1:336;plot(x,PC(1,:),'r');xlabel(&#

12、39;1979-2006年336個月')；ylabel('INDEX');title('赤道中東太平洋第一模態(tài)1979-2006年的SST時間序列')；gridon；赤道中東太平洋第一模態(tài)1979-2006年的SST時間序列l(wèi)on=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),&#

13、39;linestyle','none');colorbar;m_grid('box','fancy','tickdir','in');xlabel('longitude','fontsize',15,'fontname','comicsansms');ylabel('latitude','fontsize',15,'fontname','comicsansms');title(&

14、#39;北太平洋第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)ed北太平洋第1模態(tài)填色圖180oW120oW100oW160oW140oWlongitude0.06|0.04-0-028oN4oN0o4oS8oS160oE赤道中東太平洋第二模態(tài)1979-2006年的SST時間序列20151050-5-105030035

15、01001502002501979-2006年336個月XEDNI-150北太平洋第2模態(tài)填色圖0.05-0.05oNoNo0oSoS84048tudEo0180oW160oW140oW120oW100oWlongitude赤道中東太平洋第三模態(tài)1979-2006年的SST時間序列1050-5-10503003501001502002501979-2006年336個月XEDNI-150北太平洋第3模態(tài)填色圖tudoNoNo0oSoS840480.05-0.05-0.1o60180oW160oW140oW120oW100oWlongitudetud8oN4oN0o4oS8oStud8赤道中東太平

16、洋第四模態(tài)1979-2006年的SST時間序列-6-8-10501001502002503003501979-2006年336個月64202-XEDNI0北太平洋第4模態(tài)填色圖160oE180oW160oW140oWlongitude120oW100oW0XEDNI8oN4oN0o4oS8oS0.050-0.05642赤道中東太平洋第五模態(tài)1979-2006年的SST時間序列-4501001502002503003501979-2006年336個月-2-60北太平洋第5模態(tài)填色圖0.10.05-0.05160oE180oW160oW140oWlongitude120oW100oW三、聚類分析1

17、.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');z

18、1=linkage(y1);c1=cophenet(z1,y1);T=cluster(z1,6);H=dendrogram(z1);c1c1=0.5540TT=3122456四、功率譜分析(D太陽黑子數(shù)的離散功率譜，繪出譜圖。sunspot=importdata('sunspot.dat');year=sunspot(:,l);wolfer=sunspot(:,2);figureplot(year,wolfer)xlabel('Years');ylabel('SunspotData');title('SunspotData')ak

19、tutMDSnus1750180018501900Years19502000200SunspotDataooooooooo00642086421700figureplot(year(1:50),wolfer(1:50),'b.-');xlabel('Years');ylabel('SunspotData');title('Atthefirst50years')ataDtopsnuS20017001705171017151720172517301735174017451750Years66figuren=length(Y);powe

20、r=abs(Y(l:n/2).人2;nyquist=l/2;%取采樣頻率為0.5；freq=(l:n/2)/(n/2)*nyquist;stem(freq,power);xlabel('cycles/year');title('Periodogram')cycles/yearfigureperiod=l./freq;stem(period,power);axis(05002e+7);ylabel('Power');xlabel('Period(Years/cycle)');holdon;index=find(power=max(p

21、ower);mainPeriodStr=num2str(period(index);plot(period(index),power(index),'r.','MarkerSize',25);text(period(index)+2,power(index),'Period=',mainPeriodStr);holdoff主周期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

22、時，F(xiàn)a=3.59,IFI>Fa通過檢驗(2)利用落后自相關(guān)方法計算太陽黑子數(shù)的連續(xù)功率譜密度，繪出譜圖并給出95%信度檢驗曲線。(參考程序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),Den

23、sity(index),'r','MarkerSize',25);text(Period(index)+2,Density(index),'Period=',mainPeriodStr);holdoffxlabel('Period');ylabel('Density');holdonplot(Period,CL95,'r');holdoff五、小波分析(D請采用Morlet小波對太陽黑子活動進行小波分析，繪出小波能量譜圖。topsnus太陽黑子數(shù)年變化002601800402017501900195

24、0f°70O1850Year/a2000x=importdata('sunspot.dat');stem(x(:,1),x(:,2),'.')xlabel('Year/a');ylabel('sunspot');title('太陽黑子數(shù)年變化')>>x1=x(:,1);>>x2=x(：,2);>>coefs=cwt(x2,l:144,'morl','plot');%將尺度參數(shù)設(shè)為1:14450250AbsoluteofValuesofCa

25、,bfora=12345.9135791357913579221098876544321aselacS100150200Time(orSpace)b小波系數(shù)等值線圖：figure;c,h=contour(x1,1:144,coefs,'-');colorbarxlabel('year/a');ylabel('period/a');150100250200500-50-100-150-2006017501900195018001850year/aa/doirep1700能量譜圖fori=1:144coefsl(i,:)=coefs(i,:)./iA2;end;%小波能量contour(coefs1(2:20,:);colorbar;xlabel('year/a','fontsize',15,'fontname'