本帖最后由 vb1987 于 2013-6-12 23:27 編輯

%最近寫論文需要用到MK檢驗法，網上收集到大量的matlab代碼，但是沒有一個代碼能夠

%完全正確運行或者分析信息不全，結合多位網友編寫的MK檢驗法，經過我的改編，順利得到

%正確的運行結果，謝謝各位網友，希望對有需要的盆友有幫助

% Mann-Kendall突變檢測

% 數據序列y

% 結果序列UFk，UBk2

%--------------------------------------------

%讀取excel中的數據，賦給矩陣y

%獲取y的樣本數

%A為時間和徑流數據列

A=xlswrite('數據.xls')

x=A(:,1);%時間序列

y=A(:,2);%徑流數據列

N=length(y);

n=length(y);

% 正序列計算---------------------------------

% 定義累計量序列Sk，長度=y，初始值=0

Sk=zeros(size(y));

% 定義統(tǒng)計量UFk，長度=y，初始值=0

UFk=zeros(size(y));

% 定義Sk序列元素s

s = 0;

% i從2開始，因為根據統(tǒng)計量UFk公式，i=1時，Sk(1)、E(1)、Var(1)均為0

% 此時UFk無意義，因此公式中，令UFk(1)=0

for i=2:n

for j=1:i

if y(i)>y(j)

s=s+1;

else

s=s+0;

end;

end;

Sk(i)=s;

E=i*(i-1)/4; % Sk(i)的均值

Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差

UFk(i)=(Sk(i)-E)/sqrt(Var);

end;

% ------------------------------正序列計算end

% 逆序列計算---------------------------------

% 構造逆序列y2，長度=y，初始值=0

y2=zeros(size(y));

% 定義逆序累計量序列Sk2，長度=y，初始值=0

Sk2=zeros(size(y));

% 定義逆序統(tǒng)計量UBk，長度=y，初始值=0

UBk=zeros(size(y));

% s歸0

s=0;

% 按時間序列逆轉樣本y

% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);

for i=1:n

y2(i)=y(n-i+1);

end;

% i從2開始，因為根據統(tǒng)計量UBk公式，i=1時，Sk2(1)、E(1)、Var(1)均為0

% 此時UBk無意義，因此公式中，令UBk(1)=0

for i=2:n

for j=1:i

if y2(i)>y2(j)

s=s+1;

else

s=s+0;

end;

end;

Sk2(i)=s;

E=i*(i-1)/4; % Sk2(i)的均值

Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差

% 由于對逆序序列的累計量Sk2的構建中，依然用的是累加法，即后者大于前者時s加1，

% 則s的大小表征了一種上升的趨勢的大小，而序列逆序以后，應當表現出與原序列相反

% 的趨勢表現，因此，用累加法統(tǒng)計Sk2序列，統(tǒng)計量公式(S(i)-E(i))/sqrt(Var(i))

% 也不應改變，但統(tǒng)計量UBk應取相反數以表征正確的逆序序列的趨勢

UBk(i)=0-(Sk2(i)-E)/sqrt(Var);

end;

% ------------------------------逆序列計算end

% 此時上一步的到UBk表現的是逆序列在逆序時間上的趨勢統(tǒng)計量

% 與UFk做圖尋找突變點時，2條曲線應具有同樣的時間軸，因此

% 再按時間序列逆轉結果統(tǒng)計量UBk，得到時間正序的UBk2，做圖用

UBk2=zeros(size(y));

% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);

for i=1:n

UBk2(i)=UBk(n-i+1);

end;

% 做突變檢測圖時，使用UFk和UBk2

% 寫入目標xls文件：f:\test2.xls

% 目標表單：Sheet1

% 目標區(qū)域：UFk從A1開始，UBk2從B1開始

xlswrite('f:\test2.xls',UFk,'Sheet1','A1');

xlswrite('f:\test2.xls',UBk2,'Sheet1','B1');

figure(3)%畫圖

plot(x,UFk,'r-','linewidth',1.5);

hold on

plot(x,UBk2,'b-.','linewidth',1.5);

plot(x,1.96*ones(N,1),':','linewidth',1);

axis([min(x),max(x),-5,5]);

legend('UF統(tǒng)計量','UB統(tǒng)計量','0.05顯著水平');

xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);

ylabel('統(tǒng)計量','FontName','TimesNewRoman','Fontsize',12);

%grid on

hold on

plot(x,0*ones(N,1),'-.','linewidth',1);

plot(x,1.96*ones(N,1),':','linewidth',1);

plot(x,-1.96*ones(N,1),':','linewidth',1);

總結

以上是生活随笔為你收集整理的mk突变点检测_Mann-Kendall突变检测(mk突变检测)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。