环境：Windows 10 专业版
归类：algorithm/mk
简介：整理记录Mann-Kendall趋势检验算法，主要是趋势分析和突变点分析。
重点：魏凤英老师的《现代气候统计诊断与预测技术》中关于Mann-Kendall的突变点计算E[Sk]公式与许多论文中的公式不一致，且按照书中数据和公式绘制的图与书中的图不一致。

Mann-Kendall检验是一种非参数检验（无分布检验），其优点是不要求样本遵从一定的分布，也不受少数异常值的干扰。常用于对降水、径流、气温和水质等要素时间序列变化趋势和突变点分析。

1. 趋势分析

MK检验是检验是否拒绝零假设（null hypothesis：H₀），并接受替代假设（alternative hypothesis：H₁）：

H₀：没有单调趋势

H₁：存在单调趋势

最初的假设是：H₀为真，在拒绝H₀并接受H₁之前，数据必须要超出合理怀疑——要到达一定的置信度。

• 在MK检验中，原假设H₀为时间序列数据（X₁,X₂,...X_n），是n个独立的、随机变量同分布的样本；备选假设H₁是双边检验，对于所有的i,j≤n，且i≠j，X_i和X_j的分布是不相同的，检验的统计量S计算如下：
  

  其中，X_i、X_j分别为第i、j时间序列对应的观测值，且i<j，sgn()是符号函数：

• 当n≥8时，统计量S大致的服从正态分布，在不考虑序列存在等值数据点的情况下，其均值E(S)=0，方差V_ar(S)=n(n-1)(2n+5)/18。标准化后的检验统计量Z计算如下：
  

在双边趋势检验中，对于给定的置信水平（显著性水平）α，若|Z|≥Z_1-α/2，则原假设H₀是不可接受的，即在置信水平α（显著性检验水平）上，时间序列数据存在明显的上升或下降趋势。Z为正值表示上升趋势，负值表示减少趋势，Z的绝对值在大于等于1.645，1.96，2.576时表示分别通过了置信度90%，95%，99%的显著性检验。计算过程：以α=0.1为例，Z_1-α/2=Z_0.95，查询标准正态分布表Z_0.95=1.645，故Z≥1.645时通过90%的显著性检验，H₀假设不成立，Z>0，序列存在上升趋势。

• 方差Var简化前计算公式：
  

  减法后面方程式含义：将观测到的数据按照相同元素进行分组，g为分组数，tp为每一组元素个数，之后分别计算求和。
  举例：如一组数据（1，1，2，2，2，3，4，4，4，4），其中可以分为4组，以（元素，元素个数）格式显示，分别是（1，2个）、（2，3个）、（3，1个）、（4，4个），则g=4，依次带入求和：2(2-1)(2×2+5)+3(3-1)(2×3+5)+1(1-1)(2×1+5)+4(4-1)(2×4+5)，由公式可知，若序列中每个元素只出现一次，则求和部分结果为0，方差方程式就简化为V_ar(S)=n(n-1)(2n+5)/18。

• 衡量趋势大小的指标，用倾斜度β表示为：
  

  median表示中位值，β为正值表示“上升趋势”，β为负值表示“下降趋势”。

2.突变检测

• 设要素时间序列为X₁,X₂,...X_n，S_k表示第j个样本X_j>X_i（1≤i≤j）的累计数，定义统计量Sk：
  

• 在时间序列随机独立的假定下，Sk的均值和方差分别为：
  

  魏凤英老师的《现代气候统计诊断与预测技术》中关于Mann-Kendall的突变点计算E[Sk]公式是：k(k+1)/4。与很多发表的论文中公式都不一致，采用论文中公式。

• 将S_k标准化：
  

  其中UF₁=0，给定显著性水平α，若|UF_k|>U_α，则表明序列存在明显的趋势变化。所有UF_k可组成一条曲线。将此方法引用到反序列，将反序列X_n,X_n-1,....,X₁表示为X'₁,X'₂,....,X'_n。_j表示第j个样本X_j大于X_i（j≤i≤n）的累计数。当j'=n+1-j时，_j=r'_j，则反序列的UB_k为：

其中UB₁=0。UB_k不是简单的等于UF_k负值，而是进行了倒置再取负，此处UF_k是根据反序列算出来的。

给定显著性水平，若α=0.05，那么临界值为±1.96，绘制UFk和UBk曲线图和±1.96俩条直线再一张图上，若UFk得值大于0，则表明序列呈现上升趋势，小于0则表明呈现下降趋势，当它们超过临界直线时，表明上升或下降趋势显著。超过临界线的范围确定为出现突变的时间区域。如果UFk和UBk两条曲线出现交点，且交点在临界线内，那么交点对应的时刻便是突变开始的时间。

3.举例说明

利用经典数据：用Mann-Kendall法检测1900-1990年上海年平均气温序列，给出趋势及突变点分析，给定的显著性水平α=0.05，即U_0.05=±1.96。

20210820&003

3.1 使用EXCEL分析

3.1.1 趋势分析

步骤：（具体请参考附录3中Excel表格）

1.初始化数据

A列A2“年份”，A3-A93放<年份值>，B2“平均气温”，B3-B93放<气温值>

B列B1“年份”，C1-CO1放<年份值>，B2“平均气温”，C2-CO2放<气温值>

简单讲，就是AB列空一行放年份和平均气温，然后通过EXCEL粘贴转置成空一列在12行放置年份和平均气温，形成网格状。

2.D3插入公式：=D$2-$B3然后在网状格沿着D3对角线，全部粘贴复制，公式会自动变化。即在E4,F5,...,CO92上直接将公式粘贴复制。

3.将第3行，D3后面的行直接应用D3公式，通过EXCEL往后填充。同样，E4按行往后鼠标往后拖，自动填充。F5...，一直到CO92

4.最终填充完，会填充n(n-1)/2=91*90/2=4095个格子。C3和CO93是空的。

5.在CO列后面增加两列，CP,CQ，分别填入公式：=countif(D3:CO3,">0")，=countif(D3:CO3,"<0")，应用整列。

6.计算S=CP列和-CQ列和=2561-1306=1255

7.令n=91，（从1900-1990共91年的数据），计算V_ar(S)=n(n-1)(2n+5)/18=85085

8.由于S>0，且n≥8，计算Z=(S-1)/sqrt(V_ar(s))=4.299

9.由于Z＞2.576>0，所以趋势向上，且通过了99%的显著性检验。

倾斜度β分析

1.由上面计算出来D3-CO92区域的（X_j-X_i）值，将其分别除以（j-i），然后计算中位数即可

2.新增CR列，CR2填入“β值”,CS2-GD2填入1,2,3,....90的序列

3.在CS3填入公式：=D3/CS2，向后应用

4.在CS3对角线分别填入=E4/CS2 =F5/CS2 =G6/CS2 ...暂时不能粘贴复制，然后直接按行填充

5.在CR3填入公式：=MEDIAN(CS3:GD3) 向下填充

6.根据AB列年份、平均气温和CR列β值绘制双坐标曲线。

7.β最低点和最高点分别表示气温变化剧烈的点，并非突变点。

3.1.2 突变分析

1.在A-M列分别存储年份，平均气温，k，r，S_k，E(S_k)，V_ar(S_k)，UF(k)，逆序列，r'，S'_k，UF'(k)，UB(k)数据，以及辅助画线列N、O列“α=0.05上限（1.96）”、“α=0.05下限（-1.96）”。

2.在k列一次输入1,2,3,...,91

3.在r列输入公式=COUNTIF($B$2:B2,"<"&B3)，r₁=0，向下填充

4.在S_k列输入公式=D3+E2，S₁=0，向下填充。

5.在E(S_k)列输入公式=C3*(C3-1)/4 ，E(S₁)=0，向下填充。

6.在V_ar(S_k)列输入公式=C3*(C3-1)*(2*C3+5)/72 ，V_ar(S₁)=0，向下填充。

7.在UF(k)列输入公式=(E3-F3)/SQRT(G3) ，UF(1)=0，向下填充。

8.在逆序列中输入与平均气温顺序相反的数据

9.在r'列输入公式=COUNTIF($I$2:I2,"<"&I3) ，r'₁=0，向下填充。

10.在S'_k列输入公式=J3+K2 ，S'₁=0，向下填充。

11.在UF'(k)列输入公式=(K3-F3)/SQRT(G3) ，UF'(1)=0，向下填充。

12.将UF'(k)列的数据顺序颠倒然后加负号填入UB(k)列，UB(1)=-UF'(92)

13.选择UF(k)、UB(k)、α四列数据，绘制曲线。

3.2 使用PYTHON分析

3.2.1 趋势分析

1.需要python3.0以上

2.pip install pymannkendall

3.编写脚本

import pymannkendall as mk
data=[] #数据
mk.original_test(data)

4.结果解读

Mann_Kendall_Test(trend='increasing', h=True, p=0.020044668622627437, z=2.3255106965997814, 
      Tau=0.6, s=27.0, var_s=125.0, slope=25.75, intercept=3034.125)

trend：是否存在趋势，no trend 无趋势，increasing 上升，decreasing 下降

h：是否拒绝原假设，True 拒绝原假设H₀，False 接受原假设

p：显著性检验水平，p<0.05，越小表示接受H₁假设概率越高，即存在趋势越显著。

z：标准化后的统计量Z

tau：Kendall's tau，反映两个序列的相关性，接近1的值表示强烈的正相关，接近-1的值表示强烈的负相关

【引用网上对tau的解释，供大家参考】 Kendall's tau是数学统计中一个常用的系数，用来描述两个序列的相关系数。如果两个序列完全一致，则Kendall's tau值为1，两个毫不相关的序列的Kendall's tau值为0，而两个互逆的序列的Kendall's tau系数为-1. <br /> 具 体的计算方式为: 1 - 2 * symDif / (n * (n -1)),其中n为排列的长度(两个序列的长度相同)，symDif为对称距离。对称距离的计算方式如下: 对于两个给定的序列S1 = {a,b,c,d}; S2 = {a, c, b, d}. 分别找出两个序列的二元约束集。在这个例子中S1的所有二元约束集为{(a,b), (a,c), (a,d), (b,c), (b,d), (c,d)}, S2的所有二元约束集为{(a,c), (a,b), (a,d), (c,b), (c,d), (b,d)},比较两个二元约束集，其中不同的二元约束有两个(b,c)和(c,b),所以对称距离为2。代入上面的计算公式可以得到这两个序列的相关系 数为: 1 - 2 * 2 / (4 * 3) = 2 / 3 = 0.667

这是一个很有用的参数，可以用来比较两个序列的相似性，例如可以用于搜索引擎的排序结果的好坏。比较一个序列与一个类似标准答案的排序序列的相似性（人工评价），得出排序序列的有效性。

s：统计量S

var_s:方差V_ar(S)

slope：sen's slope，就是上文的倾斜率β，正值表示上升，负值表示下降

3.2.2 突变分析

暂无。留待有时间再研究

附录：

1. 标准正态分布表

   
   
   20210820&001
2. 置信水平
   
   20210820&002
   
   注意附录中的置信水平是1-α，等于文档中的置信水平α
3. 举例EXCEL文件存储
   202108020Mann-Kendall举例