M-K(Mann-Kendall)法是一种气候诊断与预测技术,可以判断气候序列中是否存在气候突变,如果存在,可确定出突变发生的时间。Mann-Kendall检验法也经常用于气候变化影响下的降水、干旱频次趋势检测。由于最初由曼(H.B.Mann)和肯德尔(M.G.Kendall)提出了原理并发展了这一方法,故称其为曼—肯德尔(Mann-Kendall)法。

1 原理

  1. 对于一个含有 n 个样本的时间序列 x,构造秩序列 sks_ksk​:
    sk=∑i=1kri,ri={1,xi>xj0,xi≤xj,j=1,2,...,i(1){s_k} = \displaystyle\sum_{i=1}^{k}r_i, \quad r_i = \begin{cases} 1,\quad {x_i} > {x_j}\\ 0, \quad {x_i} ≤ {x_j} \end{cases}, \quad j = 1,2,...,i \tag{1} sk​=i=1∑k​ri​,ri​={1,xi​>xj​0,xi​≤xj​​,j=1,2,...,i(1)
  2. 定义统计量UFkUF_kUFk​:
    UFk=sk−E(sk)Var(sk),k=1,2,...n(2)UF_k = \frac{s_k - E(s_k)}{\sqrt{Var(s_k)}},\quad k = 1,2,...n \tag{2} UFk​=Var(sk​)​sk​−E(sk​)​,k=1,2,...n(2)
    其中:
    E(sk)=n×(n−1)4(2-1)E(s_k) = \frac{n \times (n - 1)}{4} \tag{2-1} E(sk​)=4n×(n−1)​(2-1)
    Var(sk)=n×(n−1)×(2n+5)72(2-2)Var(s_k) = \frac{n \times (n - 1) \times (2n + 5)}{72} \tag{2-2} Var(sk​)=72n×(n−1)×(2n+5)​(2-2)
    规定 UF1=0UF_1 = 0UF1​=0。
  3. 将序列 x 倒序,重复上述步骤,获得倒序下的 UFkrUF_krUFk​r,计算 UBkUB_kUBk​:
    UBk=−UFkrUB_k = -UF_kr UBk​=−UFk​r
  4. 分析 UFkUF_kUFk​ 和 UBkUB_kUBk​值
  • UFkUF_kUFk​<0,说明持续减少趋势,值在 0.05 显著性水平线内,说明通过0.05显著性检验

  • UFkUF_kUFk​和 UBkUB_kUBk​曲线的交点在置信水平区间[-1.96 1.96]内,并且确定交点具体年份,说明该年份参数呈现突变性增长状态;

  • 如果交点不位于检验范围内,说明交点没有通过0.05 的检验,所以该年份参数突变性不具有突变性

2 代码思路

这里建立一个随机序列来模拟代码数据处理过程。

import numpy as np
np.random.seed(0)
Data = np.random.uniform(size = 72)

公共参数准备:

根据公式 2-1 和 2-2,E 和 Var 仅与数据序号 n 有关,因此,这里先准备 E 和 Var。

参考网上的解释:i 从2开始,根据统计量UFk公式,i=1时,Sk(1)、E(1)、Var(1)均为0,此时UFk无意义,因此公式中,令UFk(1) = 0

# 对时序数据 X ,生成一个序号序列,次数列范围为 2 ~ n。
# i 从 2 开始(即第二个数,其编号为 1,此时 1 处没有必要进行计算,因为其之前没有数据,所以这里从 2 开始生成)
# 规定第一个结果为0,因此我们不考虑第一个位置的结果
NUMI = np.arange(1, len(Data))
# 计算 E
# E = NUMI * (NUMI - 1) / 4
E = (NUMI + 1) * NUMI / 4
# 计算 Var
# VAR = NUMI * (NUMI - 1) * (2 * NUMI + 5) / 72
VAR = (NUMI + 1)  * NUMI * (2 * (NUMI + 1) + 5) / 72

第一步:正序计算 UFkUF_kUFk​

# 1.计算 Ri。即:序列中的某一个值与此值之前的所有值以此相比,结果为大出现的次数。
Ri = [(Data[i] > Data[:i]).sum() for i in NUMI]# 2.计算 Sk。使用 numpy 累计求和函数 cumsum。
Sk = np.cumsum(Ri)# 3.计算 UFk。考虑到 i 从 1 开始,因此把未计算的两个位置填充 0 。
UFk = np.pad((Sk - E) / np.sqrt(VAR), (1,0))

第二步:倒序计算 UBkUB_kUBk​

# 思路参考第一步,这里进行简写。
## 对于倒序,由于 Python 支持传入负数表示倒序取值,这里利用此特性直接生成倒序(反向) Bk,不包含最后一个数(编号 -1)。
Bk = np.cumsum([(Data[i] > Data[i:]).sum() for i in -(NUMI+1)])
## 按照 UFk 的计算方法后取负数即为 UBk。由于本身未对 Data 进行倒序,这里计算完成后对数据进行倒序。
UBk = np.pad((-(Bk - E) / np.sqrt(VAR)), (1,0))[::-1]

3 绘图

import matplotlib.pyplot as plt
# 配置参数
PAR = {'font.sans-serif': 'Times New Roman','axes.unicode_minus': False}
plt.rcParams.update(PAR)plt.figure(figsize = (10, 5.5), dpi = 300)
plt.plot(range(1 ,len(Data)+1),UFk,label = 'UFk',color = 'orange')
plt.plot(range(1 ,len(Data)+1),UBk,label = 'UBk',color = 'cornflowerblue')
plt.grid(True, linestyle = (0,(6,6)), linewidth = 0.4)## 画出 0.05 置信区间边界
x_lim = plt.xlim()
plt.plot(x_lim,[-1.96,-1.96],linestyle = (0,(6,6)),color = 'r')
plt.plot(x_lim, [0,0],linestyle = (0,(6,6)))
plt.plot(x_lim,[1.96,1.96],linestyle = (0,(6,6)),color = 'r')plt.legend(frameon = False)
plt.show()

4 后续

  • gma 会在 1.0.12 版本合入此算法,并对其进行扩展。欢迎各位使用、反馈和讨论。
  • 反馈请联系:Luo_Suppe(微信号)。

基于 Python 的 M-K(Mann-Kendall)突变检验 的简单实现相关推荐

  1. 精简 opencv python_基于Python的OpenCV人脸检测!简直不要太简单!

    一.文章概述 注意:本文只是人脸检测,人脸识别的实现请参见本人另一篇博客:基于OpenCV+TensorFlow+Keras实现人脸识别 本文将要讲述的是Python环境下如何用OpenCV检测人脸, ...

  2. 基于Python的ADF单位根检验方法——时间序列平稳检验

    目录 [隐藏] Abstract 平稳随机过程 肉眼检验 单位根检验 ADF Test in Python Reference Abstract 在ARMA/ARIMA这样的自回归模型中,模型对时间序 ...

  3. 基于python的证件照_python证件照换底色原来这么简单,20行代码解决!

    本文的文字及图片来源于网络,仅供学习.交流使用,不具有任何商业用途,版权归原作者所有,如有问题请及时联系我们以作处理 以下文章来源于腾讯云 作者:Python进阶者 ( 想要学习Python?Pyth ...

  4. 一个基于 Python 的简单服务|Tips

    正文共: 1479字 3图 预计阅读时间: 3分钟 一.前言 在工作中,可能遇上各种情况需要使用一个临时的 Http Server.而如果你本机上,已经存在 Python 的环境了,这样的需求基本上一 ...

  5. 基于python滑动T检验的实现--结合MK突变确定突变点

    文章目录 滑动T检验 原理 python代码 结果 结果图: step==3 结果图: step==2 相关数据可查看 滑动T检验 前面做了基于python的mk突变检验发现UF,UB曲线有较多相交点 ...

  6. 基于Python的k均值聚类不同规格的商品名

    基于Python的k均值聚类不同规格的商品名 前言 聚类的目标是使得同一簇内的点之间的距离较短,而不同簇中点之间的距离较大.以此来区分不同的群体. 本篇讲述使用k均值算法对超市购物记录集中的商品名称进 ...

  7. 基于python的股票数据的读取及可视化(K线图)

    文章目录 1.读取数据 2.绘制股票走势图 3.绘制K线图 1.读取数据 TuShare是一个免费.开源的python财经数据接口包.主要实现对股票等金融数据从数据采集.清洗加工 到 数据存储的过程, ...

  8. 基于Python的多时相数据合成

    文中的示例代码及数据可关注公众号回复20230105下载,公众号二维码见文末. 1 多时相数据合成 由于云覆盖.季节积雪.传感器故障等多种因素的影响,导致从遥感数据中提取的地表参数存在空间分布上的数据 ...

  9. Python OpenCV应用K均值聚类进行颜色量化

    Python OpenCV应用K均值聚类进行颜色量化 1. 效果图 2. 颜色量化是什么? 3. MiniBatchKMeans & KMeans 4. 源码 参考 在这篇博客文章中,我将向您 ...

最新文章

  1. 热词统计发现算法3则
  2. Java虚拟机JVM常用的几种回收算法和垃圾回收器
  3. 文献记录(part94)--Clustering and outlier detection using isoperimetric number of trees
  4. 爱护身体之简易程序员健身操
  5. hdu1520 (树形dp)
  6. c+调用java编写mq_C语言实现mq收发数据的函数
  7. UnityParticle2:5x基础模块
  8. python输入时间_一文搞懂python日期时间处理
  9. 系统字体大小导致rem布局变大
  10. How to convert any valid date string to a DateTime.
  11. linux boa post方式失败,移植boa出现的错误及解决方法
  12. 《Implicit Class-Conditioned Domain Alignment for Unsupervised Domain Adaptation》
  13. NepCTF2022 Writeup
  14. linux skb机制,skb 的分配细节
  15. 再次登顶GitHub,阿里内网首次自曝炫彩版微服务响应式与K8S手册
  16. (最简单)Java 格式化数字每3位加逗号分隔(自己封装好的工具类,直接可用)
  17. 河南省 第十一届 ACM 省赛 试题
  18. 用户6.5亿 墨迹天气难舍现金贷广告:合作方仅小米贷款
  19. Linux内核驱动开发-USB热插拔信息调取
  20. 命题逻辑在计算机中的作用,离散数学的命题逻辑.ppt

热门文章

  1. 收藏:WBS任务分解法
  2. win10下程序无法录音或使用麦克风
  3. CC(Catmull–Clark) 细分
  4. 黄金原油持续弱势震荡,本周多空思路
  5. 数字旅游——智慧景区三维可视化综合运营平台
  6. 部署filebeat收集nginx日志
  7. 【报告分享】2021年抖音双十一数据报告-蝉妈妈(附下载)
  8. 【无标题】rstudio绘制的图形导出pdf不显示文字内容,导出png正常显示
  9. 数据分析 --- 数据分析的流程
  10. 关于刷新网页F5,Ctrl+F5amp;amp;Shift+F5