用Matlab进行MK趋势分析与突变检验

%M-K趋势检定clear,closeall,clc%clear:清变数closeall:清图面clc:清画面%defineandassignthefullfilepathusingfileopendialog[filenamefilepath]=uigetfile('data1.xls');full_filepath=[filepathfilename];[X,TXTX,RAWX]=xlsread(full_filepath,1);%数据全部读入,数据缺失不影响结果x=X(:,1);%x时间nx1y=X(:,2);%y数据nx1%计算Sn=size(y,1);%数据个数S=0;fori=1:n-1S=S+sum(sign(y(i+1:n)-y(i)));%S计算式end%计算VarSVarS=n*(n-1)*(2*n+5)/18;%计算ZifS0Z=(S-1)/sqrt(VarS);elseZ=(S+1)/sqrt(VarS);end%计算Zabsalpha1=0.05;%信度95%的显著水平alpha2=0.01;%信度99%的显著水平PZ1=norminv(1-alpha1/2,0,1);PZ2=norminv(1-alpha2/2,0,1);H=0;%虚无假设Zabs=abs(Z);ifZabs=PZ1H=1;elseH=0;endP_value=2*(1-normcdf(abs(Z),0,1));%若P_value比alpha1小,则否定虚无假设%计算倾斜度ndash=n*(n-1)/2;%对称矩阵上半部slope1=zeros(ndash,1);%起始归零m=0;fork=1:n-1,forj=k+1:n,m=m+1;slope1(m)=(y(j)-y(k))/(x(j)-x(k));%分母非(j-k)end;end;slope=median(slope1);%中位数%历线绘图yd=max(y)-min(y);figureplot(x,y,'b-o','linewidth',1.5);axis([min(x),max(x),min(y)-0.2*yd,max(y)+0.2*yd]);%全距外扩20%xlabel('时间','FontName','TimesNewRoman','FontSize',12);ylabel('数据','FontName','TimesNewRoman','Fontsize',12);title('数据历线图')%添加标题gridonoutput='数据历线图';saveas(gcf,output,'jpg')%M-K突变检定Sk=zeros(size(y));%起始归零UFk=zeros(size(y));%起始归零s1=0;fori=2:nforj=1:iify(i)y(j)s1=s1+1;elses1=s1+0;end;end;Sk(i)=s1;E=i*(i-1)/4;%均值Var=i*(i-1)*(2*i+5)/72;%方差UFk(i)=(Sk(i)-E)/sqrt(Var);end;%起始归零y2=zeros(size(y));Sk2=zeros(size(y));UBk=zeros(size(y));s2=0;fori=1:ny2(i)=y(n-i+1);%逆序end;fori=2:nforj=1:iify2(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)的方差UBk(i)=-(Sk2(i)-E)/sqrt(Var);end;UBk2=zeros(size(y));fori=1:nUBk2(i)=UBk(n-i+1);%逆序end;%线性回归x1=x-x(1)+1;%x1可为非连续时间序列,有缺失数据无所谓x1nx1,非x1=[1:n]'r=corrcoef(x1,y)%相关系数R2=r(1,2)^2C=polyfit(x1,y,1)%C(1):一次项系数C(2):常数项系数%画UFk,UBkM-K统计量曲线图dFB=max(max(UFk)-min(UFk),max(UBk2)-min(UBk2));dFB1=min(min(UFk),min(UBk2))-0.2*dFB;%全距外扩20%dFB2=max(max(UFk),max(UBk2))+0.2*dFB;ifdFB1-PZ2dFB1=-5;endifdFB2PZ2dFB2=5;endfigureplot(x,UFk,'r-','linewidth',1.5);holdonplot(x,UBk2,'b-','linewidth',1.5);plot(x,PZ1*ones(n,1),'k-','linewidth',2);plot(x,PZ2*ones(n,1),'g-','linewidth',2);axis([min(x),max(x),dFB1,dFB2]);legend('UF统计量','UB统计量','0.05显著水平','0.01显著水平');xlabel('时间','FontName','TimesNewRoman','FontSize',12);ylabel('UFkUBk统计量','FontName','TimesNewRoman','Fontsize',12);title('M-K统计量曲线图')%添加标题plot(x,0*ones(n,1),'k-.','linewidth',2);plot(x,-PZ1*ones(n,1),'k-','linewidth',2);plot(x,-PZ2*ones(n,1),'g-','linewidth',2);gridon%画出主格网holdoffoutput='M-K统计量曲线图';saveas(gcf,output,'jpg')[filenamefilepath]=uigetfile('test.xls');full_filepath=[filepathfilename];xlswrite(full_filepath,x,'Sheet1','A1');xlswrite(full_filepath,UFk,'Sheet1','B1');xlswrite(full_filepath,UBk2,'Sheet1','C1');