用Matlab进行MK检验

1、精品文档用 Matlab 进行 Mann-Kendall 分析、M-K趋势分析法与M-K突变检验的原理介绍1、Mann-Kendall趋势分析法非参数检验，又称任意分布检验， 它不对变量的分布做严格假定，检验不针对特定的参数，而是模糊地对变量分布的中心位置或分布状态做检验,由于其不对总体分布做严格假定,因而适用性强12。因此，本文采用非参数的 Mann-Kendall检验法对昌马河流域近 50年的气候水文要素时间序列显著性检验，定量反映变化趋势的显著性。计算公式如下:S -1Jvar(S)Zc=0,S +1v'var( S)S. 0S =0S:二 0(1.1)在公式中,S =2 2 s

2、gn(xki4 k 4书-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为连续的气候、水文等数据序列，n为数据集合的总长度，t为每个单位的宽度，2表示所有单位的总和。衡量趋势大小的指标为：f、RXi Xj-M = Median ,Vj <i' i - j J(1.5)式中：1<j<i<n , 3所代表的是斜率，正的则是上升，负的代表下降，值的大大小代表趋势明显与否。Mann-Kenda

3、ll检验如下所示：. 。1欢迎下载精品文档5欢迎下载零假设H: 3 =0当Zc >Z12,拒绝H0假设。式中：= Z142为标准正态方差，”为显著性检验水平。2、Mann-kendall 突变检验气候系统变化是一个不稳定且不连续的变化过程，而检验其变化的常用方法之一就是 Mann-kendall突变检验方法13,该方法对于变化要素从一个相对稳定状态变化到另一个状 态的变化检验非常有效。且广泛应用于水文，气候，化学，矿物成分检验等各个方面。 Mann-kendall突变检验方法如下：对于具有n个样本量的时间序列 x,构造一个秩序列：kSk =£ ri (k =2,3，n)(1.6

4、)i 1其中1,Xi Xjni°,x。( j=1,2川,八(i.7)可见，秩序列Sk是第i时刻数值大于j时刻数值个数的累计数。在时间序列随机独立的假定下，定义统计量：sk-E(sk).UF k (k=1，2，n) (1.8)%Var )其中UF1=°, E(Sk) , Var(鼠)是累计数Sk的均值和方差，在X1，x2=，xn相互独立,且有相同连续分布时，它们可由下式算出：n(n 一1)(1.9)(1.1°)E(s k) = 一一4n(n -1)(2n 5)Var(s k)=-72UFi为标准正态分布，它是按时间序列x顺序x1,4,，Xn计算出的统计量序列，给定显

5、著性水平a,查正态分布表，若 UFi >Ua,则表明序列存在明显的趋势变化。把此方法 引用到时间序列的逆序序列中，按 xn, xn-1 ,x1,再重复上述过程，同时使UFk=-UBk,k=n, n-1 ,1, UB=Q给定显著性水平a ,将 UFk和UBk两个统计量曲线和显著性水平线 绘在同一个图上，若UFk和UBk的值大于°,则表明序列呈上升趋势，小于°则呈下降趋势。当超过临界直线时，表明上升或下降趋势显著，超过临界线的范围确定为突变的时间区域。如果UFk和UBk两条曲线出现交点，且交点在临界线之间，那么交点对应的时刻便是突变开始的时间。二、M-K程序介绍1、M-K

6、趋势检验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突变检验clc;y=data(1:50,6);%输入数据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,