Mann-Kendall（MK）检验

最新推荐文章于 2024-04-24 11:12:43 发布

水文康

最新推荐文章于 2024-04-24 11:12:43 发布

阅读量5.4k

点赞数 25

文章标签： python 机器学习开发语言数据分析概率论

本文链接：https://blog.csdn.net/weixin_51956867/article/details/136263436

版权

Mann-Kendall方法是一种用于检验时间序列趋势变化的有效方法。它的基本思想是，通过检验时间序列中的每一个数据点，来检验它们之间的相关性，从而判断时间序列是否存在趋势变化。

MK检验的原理

Mann-Kendall秩次检验方法是一种对时间序列进行趋势性检验的非参数统计检验方法。它通常用于检验趋势的显著性，因而被广泛适用于水文趋势检验研究。假定一组时间序列 $X=(x_1,x_2,x_3,...,x_n)$ ( $n$ 为变量的个数且 $n\geq{10}$ )，建立标准正态分布统计量Z:
$Z=\begin{cases} \frac{S-1}{\sqrt{Var(S)}},S>0 \\ 0,S=0 \\ \frac{S+1}{\sqrt{Var(S)}},S<0 \end{cases}$
其中：
$Var(S)=[n(n-1)(2n+5)-\sum_{p=1}^{q}{t_p(t_{p}-1)(2t_{p}+5)}]$

$S=\sum_{i=1}^{n-1}{\sum_{j=i+1}^{n}{sgn(x_{j}-x_{i})}}$

$sgn(x_{j}-x_{i})=\begin{cases} {-1},(x_j-x_i)>0 \\{0},(x_j-x_i)=0 \\-1,(x_j-x_i)<0 \end{cases}$

式中： $Va r (S)$ 为方差， $S$ 为近似服从正态分布， $q$ 为数据相同的组数， $t_p$ 为与第 $p$ 组的数据相同的个数。在趋势检验中，原假设 $H_0$ 表示数据集样本独立同分布，没有趋势存在。给定显著性水平，如果 $|Z|\geq{Z_{1-\frac{\alpha}{2}}}$ ,那么原假设 $H_0$ 将不成立，说明序列具有显著的变化趋势；如果 $Z > 0$ ,说明序列具有显著的上升趋势，反之，如果 $Z < 0$ ,说明序列具有显著的下降趋势。反之，原假设 $H_0$ 被接受。

MK检验的Python实现

MK突变检验的python检验

根据原理实现的过程如下：

导入相关的库：

from scipy import stats
import numpy as np
from matplotlib import pyplot as plt
import matplotlib

写出具体实现函数（为方便绘图，将绘图代码一并写入）：

def sk(data):
    n=len(data)
    Sk     = [0]
    UFk    = [0]
    s      =  0
    E      = [0]
    Var    = [0]
    for i in range(1,n):
        for j in range(i):
            if data[i] > data[j]:
                s = s+1
            else:
                s = s+0
        Sk.append(s)
        E.append((i+1)*(i+2)/4 )                     # Sk[i]的均值
        Var.append((i+1)*i*(2*(i+1)+5)/72 )            # Sk[i]的方差
        UFk.append((Sk[i]-E[i])/np.sqrt(Var[i]))
    UFk=np.array(UFk)
    return UFk
#a为置信度
def MK(data,a):
    ufk=sk(data)          #顺序列
    ubk1=sk(data[::-1])   #逆序列
    ubk=-ubk1[::-1]        #逆转逆序列
    #输出突变点的位置
    p=[]
    u=ufk-ubk
    for i in range(1,len(ufk)):
        if u[i-1]*u[i]<0:
            p.append(i)
    if p:
        print("突变点位置：",p)
    else:
        print("未检测到突变点")
    #画图
    conf_intveral = stats.norm.interval(a, loc=0, scale=1)   #获取置信区间
    plt.rcParams['font.family'] = ['sans-serif']
    plt.rcParams['font.sans-serif'] = ['SimHei']#解决中文不显示问题
    plt.rcParams['axes.unicode_minus']=False # 解决负号不显示问题
    plt.figure(figsize=(10,5))
    plt.plot(range(len(data)),ufk,label = 'UFk',color = 'r')
    plt.plot(range(len(data)),ubk,label = 'UBk',color = 'b')
    plt.ylabel('UFk-UBk',fontsize=25)
    x_lim = plt.xlim()
    plt.ylim([-6,7])
    plt.plot(x_lim,[conf_intveral[0],conf_intveral[0]],'m--',color='r',label='95%显著区间')
    plt.plot(x_lim,[conf_intveral[1],conf_intveral[1]],'m--',color='r')
    plt.axhline(0,ls="--",c="k")
    plt.legend(loc='upper center',frameon=False,ncol=3,fontsize=20) # 图例
    plt.xticks(fontsize=25)
    plt.yticks(fontsize=25)
    plt.show()

实例3:生成三组随机数据进行MK突变检验

#输入数据和置信度即可
# 生成3组数据
data1= np.random.random(120)
data2 = np.arange(120)
data3 = np.arange(100,0,-1)

for data in [data1,data2,data3]:
    MK(data,0.95)

结果如下：

对于data1，可得出突变点位置为 [116, 118, 119]。

在这里插入图片描述

对于data2，可得出突变点位置为[64]。

在这里插入图片描述

对于data3，可得出突变点位置为[46]。

在这里插入图片描述

实例4：基于已经写好的函数，用郑州市1956年到2016年的日降雨量数据在年尺度上进行MK突变检验，并画出日降雨量与年降雨量的时间序列图。（数据文件为4zhengzhou57083.xlsx，文件中包含温度、湿度和降雨量等多个sheet文件）

导入相关的库：

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

先读取数据

#读取数据
df = pd.read_excel('4zhengzhou57083.xlsx', sheet_name='降雨量', header=0)
dayrain=df["value"]
year=df["year"]

结果：

在这里插入图片描述

绘图：

#绘图
plt.figure()
plt.xlabel('Date')
plt.ylabel('the rainfall of the day')
plt.title('Day Rainfall Time Plot')
plt.grid(True)
plt.plot(year,dayrain)
plt.show()

在这里插入图片描述

将日数据合成为年数据

yearrain= df.groupby(df['year'])['value'].sum().reset_index(drop=True)
year1 = np.arange(1956,2017)
data = pd.concat([yearrain],axis=1)
data['year'] = year 
data

结果：

在这里插入图片描述

绘图：

plt.figure()
plt.xlabel('Date')
plt.ylabel('the rainfall of the day')
plt.title('Year Rainfall Time Plot')
plt.grid(True)
plt.plot(year1,yearrain)
plt.show()

在这里插入图片描述

进行检验：

#进行检验
MK(yearrain,0.95)

在这里插入图片描述

结果显示未找到突变点。

MK趋势检验的Python实现

Mann-Kendall趋势检验（又被称为MK检验），用于分析时间序列数据的增加或减少的趋势。它是一个非参数方法，可以适用于所有的分布（即数据不必满足正态性假设），但数据不能存在自相关性。如果数据有自相关，它可能会影响到显著水平（P值）。为了克服这个问题，研究人员提出了多个改进的Mann-Kendall检验（Hamed和Rao的修正MK检验，Yue和Wang的修正MK检验，使用预白化方法的修正MK检验）以及季节性Mann-Kendall检验。

pyMannkendal 是基于Python实现的非参数 Mann-Kendall 趋势模块，这个包有 11 个 Mann-Kendall 检验和 2 个 sen 的斜率估计函数。

1.原始 Mann-Kendall 检验 (original_test)：原始 Mann-Kendall 检验是一种非参数检验，不考虑序列相关或季节性影响。

2.Hamed 和 Rao 修正 MK 检验 (hamed_rao_modification_test)：Hamed 和 Rao (1998) 提出的这种修正 MK 检验用于解决序列自相关问题。他们建议采用方差校正方法来改进趋势分析。用户可以通过在此函数中插入滞后数来考虑前 n 个显着滞后。默认情况下，它考虑了所有重要的滞后。

3.岳和王修改的MK检验 (yue_wang_modification_test) ：这也是由 Yu, S., & Wang, C. Y. (2004) 提出的考虑序列自相关的方差校正方法。用户还可以为计算设置所需的显着 n 滞后。

4.Modified MK test using Pre-Whitening method (pre_whitening_modification_test)：Yue 和 Wang (2002) 建议在应用趋势测试之前对时间序列使用 Pre-Whitening。

5.Modified MK test using Trend free Pre-Whitening method (trend_free_pre_whitening_modification_test)：这个测试也是Yue和Wang (2002)提出的，去除趋势成分，然后在应用趋势测试之前对时间序列进行预白化。

6.Multivariate MK Test (multivariate_test)：这是Hirsch (1982)提出的多参数的MK测试。他将这种方法用于季节性 mk 测试，他将每个月视为一个参数。

7.季节性 MK 测试 (seasonal_test)：对于季节性时间序列数据，Hirsch, R.M.、Slack, J.R. 和 Smith, R.A. (1982) 提出这个检验来计算季节性趋势。

8.区域 MK 测试 (regional_test)：基于 Hirsch (1982) 提出的季节性 MK 测试，Helsel, D.R.和 Frans, L.M. (2006) 建议使用区域 mk 检验来计算区域范围内的总体趋势。

9.相关多元 MK 检验 (correlated_multivariate_test)：Hipel (1994) 提出的这种多元 MK 检验，其中参数是相关的。

10.Correlated Seasonal MK Test (correlated_seasonal_test)：Hipel (1994) 提出的这种方法，当时间序列与前一个或多个月/季显着相关时使用。

11.部分 MK 检验 (partial_test)：在真实事件中，许多因素都在影响主要研究的响应参数，这可能会使趋势结果产生偏差。为了克服这个问题，Libiseller (2002) 提出了这个部分 mk 检验。它需要两个参数作为输入，其中一个是响应参数，另一个是独立参数。

12.Theil-sen’s Slope Estimator (sens_slope)：该方法由 Theil (1950) 和 Sen (1968) 提出，用于估计单调趋势的幅度。

13.Seasonal sen’s Slope Estimator (seasonal_sens_slope)：Hipel (1994) 提出的这种方法，用于在数据具有季节性影响时估计单调趋势的幅度。

实例:生成三组随机数据基于pyMannkendal包进行MK趋势检验

import pymannkendall as mk
import numpy as np

#生成3组数据
data = np.random.random(120)
data2 = np.arange(120)
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.5417580078707189, z=0.6101567029830647, Tau=0.037815126050420166, s=270.0, var_s=194366.66666666666, slope=0.00046242202528962487, intercept=0.5047975643680629)
#Mann_Kendall_Test(trend=‘increasing’, h=True, p=0.0, z=16.192969154632337, Tau=1.0, s=7140.0, var_s=194366.66666666666, 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)

实例6：用实例2的年数据基于pyMannkendal包进行MK趋势检验

print(mk.original_test(yearrain,alpha=0.05))
print(mk.sens_slope(yearrain))

结果：

#Mann_Kendall_Test(trend=‘no trend’, h=False, p=0.7042439956814543, z=-0.37959779792429826, Tau=-0.033879781420765025, s=-62.0, var_s=25823.333333333332, slope=-0.5062500000000025, intercept=643.8875000000002)

可以发现，基于1956-2016年的数据表明，逐年的平均降雨量无显著的趋势。