任務(wù)描述：對時間序列進行MK突變檢驗：

將MK突變檢驗的代碼封裝為函數(shù)，直接調(diào)用即可，代碼如下：

%% MK突變檢驗 %% 修改日期 2022/7/29function [UF,UB] = MKbreak(time_series)n = length(time_series);%% ---------------------------------正序列計算-------------------------------- % 定義統(tǒng)計量UF，長度=n，初始值=0； UF=zeros(size(time_series));E = n*(n-1)/4; % E是累計數(shù)s(k)的均值 Var = n*(n-1)*(2*n+5)/72; % Var是累計數(shù)s(k)的方差% 定義秩序列，r(i)記錄的是第i個時刻，其數(shù)值大于j時刻(其中j<=i)數(shù)值的個數(shù)的 for i= 1:nr1(i) = sum(time_series(i)>time_series(1:i)); end% 定義累計量序列s，s(k)記錄的是第i個時刻(其中i<=k)，其數(shù)值大于j時刻(其中j<=i)數(shù)值個數(shù)的累計數(shù) s = zeros(size(time_series)); % k從2開始，因為根據(jù)統(tǒng)計量UF(k)公式，k=1時，s(1)、E(1)、Var(1)均為0，此時UF無意義，因此公式中，令UFk(1)=0 for k = 2:ns(k) = sum(r1(1:k));E = k*(k-1)/4; % s(k)的均值Var = k*(k-1)*(2*k+5)/72; % s(k)的方差UF(k) = (s(k)-E)/sqrt(Var); end%% ---------------------------------逆序列計算-------------------------------- % 按時間序列逆轉(zhuǎn)樣本,構(gòu)造逆序列time_series2和 for i=1:ntime_series2(i)=time_series(n-i+1); end % 定義逆序統(tǒng)計量UB，長度=n，初始值=0 UB = zeros(size(time_series2)); for i= 1:nr2(i) = sum(time_series2(i)>time_series2(1:i)); end% 定義逆序累計量序列s2，長度=n，初始值=0 s2 = zeros(size(time_series2)); % i從2開始，因為根據(jù)統(tǒng)計量UB(k)公式，i=1時，s2(1)、E(1)、Var(1)均為0，此時UB無意義，因此公式中，令UB(1)=0 for k = 2:ns2(k) = sum(r2(1:k));E = k*(k-1)/4; % s2(k)的均值Var = k*(k-1)*(2*k+5)/72; % s2(k)的方差 % 由于對逆序序列的累計量S2的構(gòu)建中，依然用的是累加法，即后者大于前者時r加1， % 則r的大小表征了一種上升的趨勢的大小，而序列逆序以后，應(yīng)當表現(xiàn)出與原序列相反 % 的趨勢表現(xiàn)，因此，用累加法統(tǒng)計S2序列，統(tǒng)計量公式(S(i)-E(i))/sqrt(Var(i)) % 也不應(yīng)改變，但統(tǒng)計量UBk應(yīng)取相反數(shù)以表征正確的逆序序列的趨勢UB(k) = -(s2(k)-E)/sqrt(Var); end%% ---------------------------------繪圖 % 此時上一步的到UB表現(xiàn)的是逆序列在逆序時間上的趨勢統(tǒng)計量 % 與UF做圖尋找突變點時，2條曲線應(yīng)具有同樣的時間軸，因此 % 再按時間序列逆轉(zhuǎn)結(jié)果統(tǒng)計量UB，得到時間正序的UB，做圖用 for i=1:nUB2(i)=UB(n-i+1); end x = 1:n;plot(x,UF,'g-','linewidth',1.5); hold on plot(x,UB2,'m-','linewidth',1.5); % plot(x,UB,'m-','linewidth',1.5); plot(x,1.96*ones(n,1),'-.r','linewidth',1); plot(x,0*ones(n,1),'-','color',[0.2,0.2,0.2],'linewidth',1); plot(x,-1.96*ones(n,1),'-.r','linewidth',1); grid(gca,'minor')% 畫出次格網(wǎng)axis([min(x),max(x),-5,5]); legend('UF統(tǒng)計量','UB統(tǒng)計量','0.05顯著水平'); set(gca,'Fontsize',12) set(gca,'ytick',-5:1:5) xlabel('{\itt} (year)','FontSize',15); ylabel('統(tǒng)計量','Fontsize',15);end

注：代碼并非全部原創(chuàng)，是修改自Mann-Kendall突變檢測(mk突變檢測)，優(yōu)化了部分語法，修正了部分錯誤，更正了部分描述。

突變檢驗結(jié)果如下：

圖像解讀方法如下，來自mk突變點檢測_（案例）時間序列分析之MK突變檢測

Hope this may bring you some inspirations!

后記：

寫博客的初衷是分享經(jīng)驗，同時是算是自己對思路和代碼的整理，方便日后處理數(shù)據(jù)，應(yīng)該可以幫到很多人。
我已免費分享我的心得，如果看官還有其他問題的，那么：知識付費，我的時間和經(jīng)驗正好可以解決你的問題。
咨詢問題請?zhí)砑観Q：819369354
2022年4月20日

總結(jié)

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

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