利用Matlab實現Mann-Kendall（MK）突變檢驗函數

一、MK突變檢驗

????????1、一般取顯著性水平α=0.05，那么臨界值U0.05= ±1.96?。將UFk和UBk兩個統計量序列曲線和±1.96 兩條直線均繪在一張圖上。

?????????2、若UFk和UBk的值大于0，則表明序列呈上升趨勢，小于0?則表明呈下降趨勢。?當它們超過臨界直線時，表明上升或下降趨勢顯著，超過臨界線的范圍確定為出現突變的時間區域。

????????3、如果UFk和UBk兩條曲線出現交點，且交點在臨界線之間，那么交點對應的時刻便是突變開始的時間。

?????????4、Mann-Kendall突變檢測方法的簡要計算步驟:

（1）計算順序時間序列的秩序列，按照上述公式計算UFx；

（2）計算逆序時間序列的秩序列，按照上述公式計算UBx;

（3）給定顯著性水平，如a=0.05，對于臨界值為U0.05 = ±1.96，將UFx與UB，兩個二個統計量序列曲線與U0.05 =±1.96兩條直線繪制在一個平面直角坐標α=0.10對應U0.10 = ±1.28，1=0.01對應U0.01= ±2.32。

（4）分析繪制出的UFR與UB,曲線圖，若UF,或UBx的值大于0，則表明序列呈上升趨勢，小于0則呈下降趨勢。當它們超出臨界直線時，表明上升或下降趨勢顯著。超過臨界線的范圍確定為出現突變的時間區域。若UFx與UB,兩條曲線出現交叉點，且交叉點在臨界線之間，它們交叉點對應的時刻便是突變開始的數據。

二、在Matlab上的兩種實現方法

代碼一

data1=load('D:\1\2.txt');x=data1;year=1987:2021; %% 突變檢驗for i=2:length(x)r(i)=0;for j=1:iif x(i)>x(j)r(i)=r(i)+1;endendendfor k=2:length(x)S(k)=sum(r(1:k));E(k)=k*(k-1)/4;Var(k)=k*(k-1)*(2*k+5)/72;UF(k)=(S(k)-E(k))./sqrt(Var(k));endx1=x(end:-1:1);for i=2:length(x)r1(i)=0;for j=1:iif x1(i)>x1(j)r1(i)=r1(i)+1;endendendfor k=2:length(x)S1(k)=sum(r1(1:k));E1(k)=k*(k-1)/4;Var1(k)=k*(k-1)*(2*k+5)/72;UB(k)=-(S1(k)-E1(k))./sqrt(Var1(k));end%% 繪圖figure(1)plot(year,data1)xlabel('Year','FontSize',12);ylabel('Sunspot','FontSize',12);set(gca,'FontSize',12);figure(2)plot(year,UF,'r-','MarkerSize',2,'linewidth',1.5);hold onplot(year,UB(end:-1:1),'b-','MarkerSize',2,'linewidth',1.5);plot(year,1.96*ones(length(x),1),'k--','linewidth',1);plot(year,-1.96*ones(length(x),1),'k--','linewidth',1);xlabel('Year','FontSize',12);ylabel('UF&UB','FontSize',12);set(gca,'FontSize',12);legend('UF','UB');

代碼二：

[filename,pathname] = uigetfile('*.txt','請選擇打開的數據文件'); file = [pathname, filename]; data = importdata(file); x=data(:,1);%時間序列 y=data(:,2);%數據列N=length(y);n=length(y);% 正序列計算---------------------------------% 定義累計量序列Sk，長度=y，初始值=0Sk=zeros(size(y));% 定義統計量UFk，長度=y，初始值=0UFk=zeros(size(y));% 定義Sk序列元素ss = 0;% i從2開始，因為根據統計量UFk公式，i=1時，Sk(1)、E(1)、Var(1)均為0% 此時UFk無意義，因此公式中，令UFk(1)=0for i=2:nfor j=1:iif y(i)>y(j)s=s+1;elses=s+0;endendSk(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，初始值=0y2=zeros(size(y));% 定義逆序累計量序列Sk2，長度=y，初始值=0Sk2=zeros(size(y));% 定義逆序統計量UBk，長度=y，初始值=0UBk=zeros(size(y));% s歸0s=0;% 按時間序列逆轉樣本y% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);for i=1:ny2(i)=y(n-i+1);end% i從2開始，因為根據統計量UBk公式，i=1時，Sk2(1)、E(1)、Var(1)均為0% 此時UBk無意義，因此公式中，令UBk(1)=0for i=2:nfor j=1:iif y2(i)>y2(j)s=s+1;elses=s+0;endendSk2(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:nUBk2(i)=UBk(n-i+1);end% 做突變檢測圖時，使用UFk和UBk2% 寫入目標xls文件：f:\test2.xls% 目標表單：Sheet1% 目標區域：UFk從A1開始，UBk2從B1開始xlswrite('D:\12.xls',UFk,'Sheet1','A1');xlswrite('D:\12.xls',UBk2,'Sheet1','B1');figure(3)%畫圖plot(x,UFk,'r-','linewidth',1.5);hold onplot(x,UBk2,'b-.','linewidth',1.5);plot(x,1.96*ones(N,1),':','linewidth',1);axis([min(x),max(x),-5,5]);legend('UF統計量','UB統計量','0.05顯著水平');xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);ylabel('統計量','FontName','TimesNewRoman','Fontsize',12);%grid onhold onplot(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);

三、參考文獻：

Mann-Kendall突變檢驗原理及實現

總結

以上是生活随笔為你收集整理的利用Matlab实现Mann-Kendall（MK）突变检验函数的全部內容，希望文章能夠幫你解決所遇到的問題。

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