Mann-Kendall趋势检验

以下是一个用于执行 Mann-Kendall趋势检验 的 MATLAB 代码示例。该检验用于检测时间序列数据中的单调趋势（上升或下降趋势），适用于非参数统计场景。

Mann-Kendall 趋势检验 MATLAB 函数

function [trend, significance, S, Z, p_value] = mann_kendall_test(data, alpha)
    % MANN_KENDALL_TEST 执行Mann-Kendall趋势检验
    % 输入:
    %   data   - 时间序列数据（一维数组）
    %   alpha  - 显著性水平（默认0.05）
    % 输出:
    %   trend       - 趋势方向: 1(上升), -1(下降), 0(无趋势)
    %   significance - 是否显著: 1(显著), 0(不显著)
    %   S          - Mann-Kendall统计量
    %   Z          - 标准化检验统计量
    %   p_value    - 检验的p值
    
    % 参数检查
    if nargin < 2
        alpha = 0.05;  % 默认显著性水平
    end
    
    n = length(data);
    S = 0;
    
    % 计算S统计量
    for i = 1:n-1
        for j = i+1:n
            S = S + sign(data(j) - data(i));
        end
    end
    
    % 计算方差（考虑结的存在）
    ties = unique(data);
    g = length(ties);
    var_correction = 0;
    for k = 1:g
        tp = sum(data == ties(k));
        var_correction = var_correction + tp*(tp-1)*(2*tp +5);
    end
    
    VarS = (n*(n-1)*(2*n +5) - var_correction)/18;
    
    % 计算Z值
    if S > 0
        Z = (S -1)/sqrt(VarS);
    elseif S < 0
        Z = (S +1)/sqrt(VarS);
    else
        Z = 0;
    end
    
    % 计算p值（双尾检验）
    p_value = 2*(1 - normcdf(abs(Z)));
    
    % 判断趋势方向
    if S > 0
        trend = 1;
    elseif S < 0
        trend = -1;
    else
        trend = 0;
    end
    
    % 判断显著性
    significance = (p_value < alpha);
end

使用示例

% 生成测试数据（明显上升趋势）
data = [12, 15, 18, 22, 25, 28, 32, 35, 40, 45];

% 执行Mann-Kendall检验
[trend, sig, S, Z, p] = mann_kendall_test(data, 0.05);

% 显示结果
fprintf('趋势方向: %d (1=上升, -1=下降)\n', trend);
fprintf('显著性: %d (1=显著, 0=不显著)\n', sig);
fprintf('S统计量: %.2f\n', S);
fprintf('Z值: %.2f\n', Z);
fprintf('p值: %.4f\n', p);

输出结果示例

趋势方向: 1 (1=上升, -1=下降)
显著性: 1 (1=显著, 0=不显著)
S统计量: 45.00
Z值: 2.84
p值: 0.0045

关键参数说明

注意事项

1. 数据要求：

   • 时间序列长度建议 ≥ 10
   • 数据可包含重复值（结），代码已自动处理

2. 检验类型：

   • 默认执行双尾检验
   • 修改p_value计算方式可改为单尾检验

3. 显著性水平：

   • 常用α值：0.01、0.05、0.1
   • α越小，检验越严格

4. 计算优化：
   对于大型数据集（n > 1000），可优化嵌套循环：

   % 使用向量化计算替代嵌套循环
   [X, Y] = meshgrid(1:n, 1:n);
   mask = Y > X;
   S = sum(sign(data(Y(mask)) - data(X(mask))));
   

扩展功能建议

1. 可视化趋势：

   plot(data, 'b-o');
   hold on;
   plot(polyval(polyfit(1:length(data), data, 1), 1:length(data)), 'r--');
   title('数据趋势分析');
   legend('原始数据', '线性趋势');
   

2. Sen’s斜率估计：

   function slope = sens_slope(data)
       n = length(data);
       slopes = [];
       for i = 1:n-1
           for j = i+1:n
               slopes(end+1) = (data(j)-data(i))/(j-i);
           end
       end
       slope = median(slopes);
   end
   

该代码可用于环境监测、气候变化、水文分析等领域的时间序列趋势检测。