% Mann-Kendall突變檢測

% 數據序列y

% 結果序列UFk，UBk2

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

%獲取y的樣本數

%A為時間和降水數據列

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

y=降水(:,2);%降水數據列

N=length(x);

n=length(y);

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

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

Sk=zeros(size(y));

% 定義統計量UFk，長度=y，初始值=0

UFk=zeros(size(y));

% 定義Sk序列元素s

s = 0;

% i從2開始，因為根據統計量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));

% 定義逆序統計量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開始，因為根據統計量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的大小表征了一種上升的趨勢的大小，而序列逆序以后，應當表現出與原序列相反

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

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

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

end;

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

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

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

% 再按時間序列逆轉結果統計量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

figure(3)%畫圖

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

hold on

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

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

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

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

legend('UF統計量','UB統計量');

xlabel('年份','FontName','TimesNewRoman','FontSize',9);

ylabel('統計量','FontName','TimesNewRoman','Fontsize',9);

%grid on

hold on

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

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

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

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

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

time=1951:6:2014

Xlabel('年份');

Ylabel('統計量');

總結

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

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