用Matlab進(jìn)行MK趨勢(shì)分析與突變檢驗(yàn)

《用Matlab進(jìn)行MK趨勢(shì)分析與突變檢驗(yàn)》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在應(yīng)用文檔-天天文庫(kù)。

1、%M-K趨勢(shì)檢定clear,closeall,clc%clear:清變數(shù)closeall:清圖面clc:清畫面%defineandassignthefullfilepathusing"fileopen"dialog[filenamefilepath]=uigetfile('data1.xls');full_filepath=[filepathfilename];[X,TXTX,RAWX]=xlsread(full_filepath,1);%數(shù)據(jù)全部讀入,數(shù)據(jù)缺失不影響結(jié)果x=X(:,1);%x時(shí)間y=X(:,2);%y數(shù)據(jù)%計(jì)算

2、Sn=size(y,1);%數(shù)據(jù)個(gè)數(shù)S=0;fori=1:n-1S=S+sum(sign(y(i+1:n)-y(i)));%S計(jì)算式end%計(jì)算VarSVarS=n*(n-1)*(2*n+5)/18;%計(jì)算ZifS>0Z=(S-1)/sqrt(VarS);elseZ=(S+1)/sqrt(VarS);end%計(jì)算Zabsalpha1=0.05;%信度95%的顯著水平alpha2=0.01;%信度99%的顯著水平PZ1=norminv(1-alpha1/2,0,1);PZ2=norminv(1-alpha2/2,0,1);H=0;%虛無假設(shè)Zabs=

3、abs(Z);ifZabs>=PZ1H=1;elseH=0;endP_value=2*(1-normcdf(abs(Z),0,1));%若P_value比alpha1小,則否定虛無假設(shè)%計(jì)算傾斜度ndash=n*(n-1)/2;%對(duì)稱矩陣上半部slope1=zeros(ndash,1);%起始?xì)w零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);%中位數(shù)%歷線繪圖yd=max(y)-min(y

4、);figureplot(x,y,'b-o','linewidth',1.5);axis([min(x),max(x),min(y)-0.2*yd,max(y)+0.2*yd]);%全距外擴(kuò)20%xlabel('時(shí)間','FontName','TimesNewRoman','FontSize',12);ylabel('數(shù)據(jù)','FontName','TimesNewRoman','Fontsize',12);title('數(shù)據(jù)歷線圖')%添加標(biāo)題gridonoutput='數(shù)據(jù)歷線圖';saveas(gcf,output,'jpg')%M-K突變檢

5、定Sk=zeros(size(y));%起始?xì)w零UFk=zeros(size(y));%起始?xì)w零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;%起始?xì)w零y2=zeros(size(y));Sk2=zeros(size(y));UBk=zeros(size(y));s2=0;fori=1:ny2(i)=y(n-i+1)

6、;%逆序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可為非連續(xù)時(shí)間序列,有缺失數(shù)據(jù)無所謂x1,非x1=[1:n]'r=co

7、rrcoef(x1,y)%相關(guān)系數(shù)R2=r(1,2)^2C=polyfit(x1,y,1)%C(1):一次項(xiàng)系數(shù)C(2):常數(shù)項(xiàng)系數(shù)%畫UFk,UBkM-K統(tǒng)計(jì)量曲線圖dFB=max(max(UFk)-min(UFk),max(UBk2)-min(UBk2));dFB1=min(min(UFk),min(UBk2))-0.2*dFB;%全距外擴(kuò)20%dFB2=max(max(UFk),max(UBk2))+0.2*dFB;ifdFB1>-PZ2dFB1=-5;endifdFB2

8、'linewidth',1.5);holdonplot(x,UBk2,'b-','linewidth',1.5);plot