% Mann-Kendall突變檢測

% 數(shù)據(jù)序列y

% 結(jié)果序列UFk，UBk2

%讀取excel中的數(shù)據(jù)，賦給矩陣y

%獲取y的樣本數(shù)

%A為時(shí)間和降水?dāng)?shù)據(jù)列

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

y=降水(:,2);%降水?dāng)?shù)據(jù)列

N=length(x);

n=length(y);

% 正序列計(jì)算---------------------------------

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

Sk=zeros(size(y));

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

UFk=zeros(size(y));

% 定義Sk序列元素s

s = 0;

% i從2開始，因?yàn)楦鶕?jù)統(tǒng)計(jì)量UFk公式，i=1時(shí)，Sk(1)、E(1)、Var(1)均為0

% 此時(shí)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;

% ------------------------------正序列計(jì)算end

% 逆序列計(jì)算---------------------------------

% 構(gòu)造逆序列y2，長度=y，初始值=0

y2=zeros(size(y));

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

Sk2=zeros(size(y));

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

UBk=zeros(size(y));

% s歸0

s=0;

% 按時(shí)間序列逆轉(zhuǎn)樣本y

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

for i=1:n

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

end;

% i從2開始，因?yàn)楦鶕?jù)統(tǒng)計(jì)量UBk公式，i=1時(shí)，Sk2(1)、E(1)、Var(1)均為0

% 此時(shí)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)的方差

% 由于對逆序序列的累計(jì)量Sk2的構(gòu)建中，依然用的是累加法，即后者大于前者時(shí)s加1，

% 則s的大小表征了一種上升的趨勢的大小，而序列逆序以后，應(yīng)當(dāng)表現(xiàn)出與原序列相反

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

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

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

end;

% ------------------------------逆序列計(jì)算end

% 此時(shí)上一步的到UBk表現(xiàn)的是逆序列在逆序時(shí)間上的趨勢統(tǒng)計(jì)量

% 與UFk做圖尋找突變點(diǎn)時(shí)，2條曲線應(yīng)具有同樣的時(shí)間軸，因此

% 再按時(shí)間序列逆轉(zhuǎn)結(jié)果統(tǒng)計(jì)量UBk，得到時(shí)間正序的UBk2，做圖用

UBk2=zeros(size(y));

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

for i=1:n

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

end;

% 做突變檢測圖時(shí)，使用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統(tǒng)計(jì)量','UB統(tǒng)計(jì)量');

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

ylabel('統(tǒng)計(jì)量','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('統(tǒng)計(jì)量');

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。