时间序列数据的Mann-Kendall突变分析

stwen031

function [UF,UB2,U_Alpha] = MannKendallMutation1(X,Alpha,OutIndex)
% 时间序列数据的Mann-Kendall突变分析
%
% Input
% X                             - 时间序列数据，向量,Time Series of Data
% Alpha                         - 显著性水平
% OutIndex                      - 画图指示，默认OutIndex=1
%
% Output
% UF                            - UF统计量序列,1*N
% UB2                            - UB统计量序列,1*N
% U_Alpha                       - 显著性水平Alpha下的正态分布分位数
%
%
% 作者：***
% E_mail:******
% 时间：2017年10月18日
% 修改：2018年05月02日
%
%  -----------------------------------------------------

if (nargin < 3) || (isempty(OutIndex)), OutIndex = 1; end
if (nargin < 2) || (isempty(Alpha)), Alpha = 0.05; end

load('Y.mat')
N = length(X);

%%
% 正序列计算
% 在时间序列X为独立随机的假设下，UF服从标准正态分布
% 计算UF统计量
UF = zeros(1,N);                       % 时间序列X的UF统计量
Sk = zeros(1,N);                       % 时间序列X的秩序列(累计量序列)
s = 0;

for ii = 2:N
    for jj = 1:(ii-1)
        if X(ii) > X(jj)
            s = s + 1;
        end
    end
    Sk(ii) = s;

    % ESk = ii * (ii + 1) / 4;           % 文献中公式为加号：徐建华
    ESk = ii * (ii - 1) / 4;
    VarSk = ii * (ii - 1) * (2 * ii + 5) / 72;

    UF(ii) = (Sk(ii) - ESk) / sqrt(VarSk);
end

%%
% 由于对逆序序列的累计量Sk2的构建中依然用的是累加法
% 即后者大于前者时s加1，则s的大小表征了一种上升的趋势的大小
% 而序列逆序以后，应当表现出与原序列相反的趋势表现
% 因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))也不应改变
% 但统计量UB应取相反数以表征正确的逆序序列的趋势
%
% 在时间序列X为独立随机的假设下，UB服从标准正态分布
% 计算UB统计量
Y = X(N:-1:1);                         % 取时间序列X的逆序
% Y = zeros(1,N);                        % 取时间序列X的逆序
% for ii = 1:N
%     Y(ii) = X(N - ii + 1);
% end

UB = zeros(1,N);                       % 时间序列X的UB统计量
Sk2 = zeros(1,N);                      % 时间序列Y的秩序列
s = 0;

for ii = 2:N
    for jj = 1:(ii-1)
        if Y(ii) > Y(jj)
            s = s + 1;
        end
    end
    Sk2(ii) = s;

    % ESk = ii * (ii + 1) / 4;           % 文献中公式为加号：徐建华
    ESk = ii * (ii - 1) / 4;
    VarSk = ii * (ii - 1) * (2 * ii + 5) / 72;

    UB(ii) = 0 - (Sk2(ii) - ESk) / sqrt(VarSk);
end

%%
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
% 与UF作图寻找突变点时，2条曲线应具有同样的时间轴
% 因此再按时间序列逆转结果统计量UB，得到时间正序的UB2，作图用
%
UB2 = UB(N:-1:1);

% 也可以使用UB2=flipud(UB);或者UB2=flipdim(UB,2);
% UB2 = zeros(1,N);
% for ii = 1:N
%     UB2(ii) = UB(N - ii + 1);
% end

%%
% 若UF > U_Alpha, 表明序列存在明显的趋势变化
U_Alpha = norminv(1 - Alpha / 2);

xlswrite('MK_UFUB.xlsx',UF','Sheet1','A1');
xlswrite('MK_UFUB.xlsx',UB2','Sheet1','B1');

if OutIndex == 1
    YearID = 1:N;
    figure
    plot(YearID,UF,'r-','linewidth',1.5);
    % plot(YearID,UF,'r-');
    hold on
    plot(YearID,UB2,'b-.','linewidth',1.5);
    % plot(YearID,UB2,'b-.');
    legend('UF统计量','UB统计量');
    % xlabel('t(year)','FontName','TimesNewRoman','FontSize',12);
    xlabel('t(year)');
    ylabel('统计量');
    % ylabel('统计量','FontName','TimesNewRoman','FontSize',12);
    hold on
    plot(YearID,U_Alpha*ones(1,N),':','linewidth',1);
    plot(YearID,-U_Alpha*ones(1,N),':','linewidth',1');
    plot(YearID,zeros(1,N),'-.','linewidth',1);
    axis([1,N,-6,6]);
    % grid on
end

% end of function MannKendallMutation
%%