【森气杂谈】python利用pymannkendall包进行MK（Mann-Kendall）趋势检验

【森气杂谈】python利用pymannkendall包进行MK（Mann-Kendall）趋势检验

气象学中常用的Mann-Kendall趋势检验，是一种非参数统计检验方法。该方法可用于分析中心趋势不稳定的时间序列，基于数据的秩，而不是数据本身。Mann-Kendall趋势检验适用于分析持续增长或下降趋势(单调趋势)的时间序列数据。适用于所有的分布(即数据不需要满足正态分布的假设)，但数据应该没有序列相关性。如果数据具有序列相关性，则会在显著性水平(p值)上产生影响。
更进一步的原理和公式可以参考这个博客。
今天要介绍的是一个python第三方包----pymannkendall，该包集成了MK（Mann-Kendall）趋势检验，而不需要自己写公式定义函数。

安装pymannkendall

首先需要安装pymannkendall

平台	方法
Windows	pip install pymannkendall
conda	conda install -c conda-forge pymannkendall
Linux	sudo pip install pymannkendall

实例----普通MK趋势检验

调用pymannkendall

# -*- encoding: utf-8 -*-
'''
@File    :   mk_demo.py
@Time    :   2021/11/27 12:01:23
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''
# here put the import lib
import pymannkendall as mk
import numpy as np

# 生成3组数据
data = np.random.random(100)
data2 = np.arange(100)
data3 = np.arange(100,0,-1)

for x in [data,data2,data3]:
    result = mk.original_test(x,alpha=0.05)#alpha默认为0.05
    print(result)
    

结果

Mann_Kendall_Test(trend='no trend', h=False, p=0.823256299016955, z=-0.22335875907380243, Tau=-0.015353535353535354, s=-76.0, var_s=112750.0, slope=-0.0002038833935962383, intercept=0.4745895001589582)
Mann_Kendall_Test(trend='increasing', h=True, p=0.0, z=14.73869998208331, Tau=1.0, s=4950.0, var_s=112750.0, slope=1.0, intercept=0.0)
Mann_Kendall_Test(trend='decreasing', h=True, p=0.0, z=-14.73869998208331, Tau=-1.0, s=-4950.0, var_s=112750.0, slope=-1.0, intercept=100.0)

当然我们也可以直接print我们需要的参数即可。

print(f'trend:{res.trend}','p_value:{:.2f}'.format(res.p),'slope:{:.2f}'.format(res.slope),sep = ',')

trend:no trend,p_value:0.74,slope:0.00
trend:increasing,p_value:0.00,slope:1.00
trend:decreasing,p_value:0.00,slope:-1.00

除了普通的mk趋势检验外，针对空间条件开发了其他几种修正的Mann-Kendall检验，如多元MK检验、区域MK检验、相关MK检验、部分MK检验等。具体请参考pymannkendall官方文档。

自定义的mk趋势检验函数

这边是我参考这篇博客修改后添加了趋势拟合后的代码

# -*- encoding: utf-8 -*-
'''
@File    :   mk_trend.py
@Time    :   2021/08/15 20:31:49
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''

# here put the import lib
from scipy.stats import norm
import numpy as np
from sklearn.linear_model import LinearRegression
 
 
def mk(x, alpha=0.1):
    n = len(x)

    # 计算趋势slope
    model = LinearRegression()
    model.fit(np.arange(1,n+1).reshape(-1,1),x)
    slope = model.coef_[0]
 
    # 计算S的值
    s = 0
    for j in range(n - 1):
        for i in range(j + 1, n):
            s += np.sign(x[i] - x[j])
 
    # 判断x里面是否存在重复的数，输出唯一数队列unique_x,重复数数量队列tp
    unique_x, tp = np.unique(x, return_counts=True)
    g = len(unique_x)
 
    # 计算方差VAR(S)
    if n == g:  # 如果不存在重复点
        var_s = (n * (n - 1) * (2 * n + 5)) / 18
    else:
        var_s = (n * (n - 1) * (2 * n + 5) - np.sum(tp * (tp - 1) * (2 * tp + 5))) / 18
 
    # 计算z_value
    if n <= 10:  # n<=10属于特例
        z = s / (n * (n - 1) / 2)
    else:
        if s > 0:
            z = (s - 1) / np.sqrt(var_s)
        elif s < 0:
            z = (s + 1) / np.sqrt(var_s)
        else:
            z = 0
 
    # 计算p_value，可以选择性先对p_value进行验证
    p = 2 * (1 - norm.cdf(abs(z)))
 
    # 计算Z(1-alpha/2)
    h = abs(z) > norm.ppf(1 - alpha / 2)
 
    # 趋势判断
    if (z < 0) and h:
        trend = -1#'decreasing'
    elif (z > 0) and h:
        trend = 1#'increasing'
    else:
        trend = 0#'no trend'
 
    return trend,p,slope


if __name__ == '__main__':
    data = np.random.random(100)
    data2 = np.arange(100)
    data3 = np.arange(100,0,-1)
    for x in [data,data2,data3]:
        trend,p,slope = mk(x,alpha = 0.05)
        print(f'trend:{trend}','p_value:{:.2f}'.format(p),'slope:{:.2f}'.format(slope),sep = ',')

结果

trend:0,p_value:0.23,slope:0.00
trend:1,p_value:0.00,slope:1.00
trend:-1,p_value:0.00,slope:-1.00

对比之下可以发现利用pymannkendall，MK趋势检验更简单方便。

参考文献

趋势检验方法（二）MK趋势检验
Mann-Kendall’ test 曼－肯德尔趋势检验法
《现在气候统计诊断与预测技术（第2版）魏凤英编》