用Matlab進(jìn)行MK檢驗(yàn)

1、精品文檔用 Matlab 進(jìn)行 Mann-Kendall 分析、M-K趨勢(shì)分析法與M-K突變檢驗(yàn)的原理介紹1、Mann-Kendall趨勢(shì)分析法非參數(shù)檢驗(yàn)，又稱(chēng)任意分布檢驗(yàn)， 它不對(duì)變量的分布做嚴(yán)格假定，檢驗(yàn)不針對(duì)特定的參數(shù)，而是模糊地對(duì)變量分布的中心位置或分布狀態(tài)做檢驗(yàn),由于其不對(duì)總體分布做嚴(yán)格假定,因而適用性強(qiáng)12。因此，本文采用非參數(shù)的 Mann-Kendall檢驗(yàn)法對(duì)昌馬河流域近 50年的氣候水文要素時(shí)間序列顯著性檢驗(yàn)，定量反映變化趨勢(shì)的顯著性。計(jì)算公式如下:S -1Jvar(S)Zc=0,S +1v'var( S)S. 0S =0S:二 0(1.1)在公式中,S =2 2 s

2、gn(xki4 k 4書(shū)-Xi)1,0 >0sgn(8)=«0,9=0-1,0 <0n J n(1.2)(1.3)varS = |n(n T)(2n+5) -Z t(t 1)(2t + 5) ,.18t(1.4)式中：xk, xi為連續(xù)的氣候、水文等數(shù)據(jù)序列，n為數(shù)據(jù)集合的總長(zhǎng)度，t為每個(gè)單位的寬度，2表示所有單位的總和。衡量趨勢(shì)大小的指標(biāo)為：f、RXi Xj-M = Median ,Vj <i' i - j J(1.5)式中：1<j<i<n , 3所代表的是斜率，正的則是上升，負(fù)的代表下降，值的大大小代表趨勢(shì)明顯與否。Mann-Kenda

3、ll檢驗(yàn)如下所示：. 。1歡迎下載精品文檔5歡迎下載零假設(shè)H: 3 =0當(dāng)Zc >Z12,拒絕H0假設(shè)。式中：= Z142為標(biāo)準(zhǔn)正態(tài)方差，”為顯著性檢驗(yàn)水平。2、Mann-kendall 突變檢驗(yàn)氣候系統(tǒng)變化是一個(gè)不穩(wěn)定且不連續(xù)的變化過(guò)程，而檢驗(yàn)其變化的常用方法之一就是 Mann-kendall突變檢驗(yàn)方法13,該方法對(duì)于變化要素從一個(gè)相對(duì)穩(wěn)定狀態(tài)變化到另一個(gè)狀 態(tài)的變化檢驗(yàn)非常有效。且廣泛應(yīng)用于水文，氣候，化學(xué)，礦物成分檢驗(yàn)等各個(gè)方面。 Mann-kendall突變檢驗(yàn)方法如下：對(duì)于具有n個(gè)樣本量的時(shí)間序列 x,構(gòu)造一個(gè)秩序列：kSk =£ ri (k =2,3，n)(1.6

4、)i 1其中1,Xi Xjni°,x。( j=1,2川,八(i.7)可見(jiàn)，秩序列Sk是第i時(shí)刻數(shù)值大于j時(shí)刻數(shù)值個(gè)數(shù)的累計(jì)數(shù)。在時(shí)間序列隨機(jī)獨(dú)立的假定下，定義統(tǒng)計(jì)量：sk-E(sk).UF k (k=1，2，n) (1.8)%Var )其中UF1=°, E(Sk) , Var(鼠)是累計(jì)數(shù)Sk的均值和方差，在X1，x2=，xn相互獨(dú)立,且有相同連續(xù)分布時(shí)，它們可由下式算出：n(n 一1)(1.9)(1.1°)E(s k) = 一一4n(n -1)(2n 5)Var(s k)=-72UFi為標(biāo)準(zhǔn)正態(tài)分布，它是按時(shí)間序列x順序x1,4,，Xn計(jì)算出的統(tǒng)計(jì)量序列，給定顯

5、著性水平a,查正態(tài)分布表，若 UFi >Ua,則表明序列存在明顯的趨勢(shì)變化。把此方法 引用到時(shí)間序列的逆序序列中，按 xn, xn-1 ,x1,再重復(fù)上述過(guò)程，同時(shí)使UFk=-UBk,k=n, n-1 ,1, UB=Q給定顯著性水平a ,將 UFk和UBk兩個(gè)統(tǒng)計(jì)量曲線(xiàn)和顯著性水平線(xiàn) 繪在同一個(gè)圖上，若UFk和UBk的值大于°,則表明序列呈上升趨勢(shì)，小于°則呈下降趨勢(shì)。當(dāng)超過(guò)臨界直線(xiàn)時(shí)，表明上升或下降趨勢(shì)顯著，超過(guò)臨界線(xiàn)的范圍確定為突變的時(shí)間區(qū)域。如果UFk和UBk兩條曲線(xiàn)出現(xiàn)交點(diǎn)，且交點(diǎn)在臨界線(xiàn)之間，那么交點(diǎn)對(duì)應(yīng)的時(shí)刻便是突變開(kāi)始的時(shí)間。二、M-K程序介紹1、M-K

6、趨勢(shì)檢驗(yàn)function slope,zc,za,sign=MannKendall(x)n n%十算 s % S = "'、sgn(Xk 一 Xi )i"1s=0;len=size(x,2);for m=1:len-1for n=m+1:lenif x(n)>x(m)s=s+1;elseif x(n)=x(m)s=s+0;elses=s-1;endendendvarS = n(n -1)(2n 5)，t(t -1)(2t 5) 18刈算vars和e% -t二vars=len*(len-1)*(2*len+5)/18;S 0S =0S :二 0S -1Jvar(

7、 S)Zc=<0,S-1%十算zc %了兩if s>0zc=(s-1)/sqrt(vars);elsezc=(s+1)/sqrt(vars);end%十算zaza=var(x);sign=0;zc1=abs(zc);if zc1 >= zasign=1;elsesign=0;end% 十算傾斜度ndash = len * ( len - 1 ) / 2; slope1= zeros( ndash, 1 ); m=1;for k = 1:len-1,for j = k+1:len,slope1(m) = ( x(j) - x(k) ) / ( j - k ) ; % m = m

8、 + 1;end;end;Xi -XjM = Median , Vj < iI i-j Jslope= median( slopel );2、M-K突變檢驗(yàn)clc;y=data(1:50,6);%輸入數(shù)據(jù)n=length(y);Sk=zeros(size(y); UFk=zeros(size(y);kSk =、 ri (k =2,3,n)i =1s1=0;for i=2:n%for j=1:iif y(i)>y(j) s1=s1+1;elses1=s1+0;%end;亞；：1,Xi Xj ri10, x - xj(j=1,2,lll,i)end;Sk(i尸s1;%E=i*(i-1)

9、/4;Var=i*(i-1)*(2*i+5)/72;UFk(i)=(Sk(i)-E)/sqrt(Var);E(Sk)=，4UFk = %;(:; (k=1,2,n) end;y2=zeros(size(y);Sk2=zeros(size(y); UBk=zeros(size(y); s2=0;for i=1:ny2(i)=y(n-i+1);end;for i=2:nfor j=1:iif y2(i)>y2(j)s2=s2+1;elses2=s2+0;end;end;Sk2(i)=s2;E=i*(i-1)/4; % Sk2(i) 的均值Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差 % Var(s k)=n (n -1)(2 n 5)72UBk(i)=0-(Sk2(i)-E)/sqrt(Var);end;UBk2=zeros(size(y);for i=1:nUBk2(i)=UBk(n-i+1);end;year=1961:2010;year=year'xlswrite('D:test.xls',year,'Sheet1','A1');xlswrite('D:test.xls',UFk,