python中的Mann-Kendall单调趋势检验--及原理说明

https://blog.csdn.net/liuchengzimozigreat/article/details/87931248
在python中使用Mann-kendall：
Q: Using Mann Kendall in python with a lot of data
在python中使用mann-Kendall，可以用scipy.stats.kendalltau，该函数返回两个值：tau-反映两个序列的相关性，接近1的值表示强烈的正相关，接近-1的值表示强烈的负相关；p_value：p值反映的是假设检验的双边p值，其零假设为无关联——即通常所谓的显著性水平，一般取p<0.05为显著。
其实上面重要的是显著性水平，至于正相关还是负相关我们可以用斜率来表示。斜率可用线性回归的斜率，也可以是Sen’s slope(python中可以用stat包提供的函数来计算)。

---

??想了解关于风速历年变化趋势，被推荐用Mann-Kendall趋势检验，查到一篇相对详细的介绍，转译成中文以飨诸君。原文：Mann-Kendall Test For Monotonic Trend
??水平有限，如有纰漏，还望斧正。

背景知识

??Mann-Kendall(MK)检验(test)(Mann 1945, Kendall 1975, Gilbert 1987) 的目的是统计评估我们所感兴趣的变量，随着时间变化，是否有单调上升或下降的趋势。单调上升（下降）的趋势意味着该变量随时间增加（减少），但此趋势可能是、也可能不是线性的。MK test可替代参数线性回归分析——线性回归可检验线性拟合直线的斜率是否不为零。回归分析要求拟合回归线的残差是正态分布的，MK检验不需要这种假设，MK检验是非参数检验（不要求服从任何分布-distribution free）
??Hirsch, Slack and Smith (1982, page 107)表明MK检验最好被视作探索性分析，最适合用于识别变化显着或幅度较大的站点，并量化这些结果。

相关前提（assumption）

以下这些假设是MK检验的基础：

• 当没有趋势时，随时间获得的数据是独立同分布的。独立的假设是说数据随着时间不是连续相关的。
• 所获得的时间序列上的数据代表了采样时的真是条件。（样本具有代表性）
• 样本的采集、处理和测量方法提供了总体样本中的无偏且具有代表性的观测值。

MK检验不要求数据是正态分布，也不要求变化趋势——如果存在的话——是线性的。如果有缺失值或者值低于一个或多个检测限制，是可以计算MK检测的，但检测性能会受到不利影响。独立性假设要求样本之间的时间足够大，这样在不同时间收集的测量值之间不存在相关性。

计算

MK检验是检验是否拒绝零假设(null hypothesis: H0)，并接受替代假设（alternative hypothesis: Ha）：
??H0：没有单调趋势
??Ha：存在单调趋势
最初的假设是：H0为真，在拒绝H0并接受Ha之前，数据必须要超出合理怀疑——要到达一定的置信度。
MK检验的流程 (from Gilbert 1987, pp. 209-213)：
??1.将数据按采集时间列出：x1,x2,…,xn,，即分别在时间1,2,…,n得到的数据。

??2.确定所有n(n-1)/2个xj - xk差值的符号，其中j > k，这些差值是：x2 - x1，x3 - x1, … , xn - x1，x3 - x2，x4 - x2，…，xn - xn-2，xn - xn-1

??3.令sgn(xj - xk)作为指示函数，依据xj - xk的正负号取值为1,0或-1，即：
s g n ( x j ? x k ) = { 1 , x j ? x k > 0 0 , x j ? x k = 0 或 者 x j ? x k 的 值 因 没 检 测 （ 数 据 缺 失 ） 而 不 能 确 定 ； ? 1 , x j ? x k < 0 ( 3 ) \boldsymbol{sgn(x_j-x_k)} =

?????1,xj?xk>00,xj?xk=0或者xj?xk的值因没检测（数据缺失）而不能确定；?1,xj?xk<0{1,xj?xk>00,xj?xk=0或者xj?xk的值因没检测（数据缺失）而不能确定；?1,xj?xk<0

\quad\quad\quad\boldsymbol{(3)}sgn(xj?xk)=?????1,xj?xk>00,xj?xk=0或者xj?xk的值因没检测（数据缺失）而不能确定；?1,xj?xk<0(3)
??例如：如果xj - xk > 0，就意味着 j 时刻的观测值——用xj表示，大于 k 时刻的观测值——用xk表示.

??4.计算

????????S=${\sum }_{k=1}^{\mathbf{n}?1}{\sum }_{\mathbf{j}=\mathbf{k}+1}^{\mathbf{n}}$sgn(xj - xk) (1)
??即差值为正的数量减去差值为负的数量。如果S是一个正数，那么后一部分的观测值相比之前的观测值会趋向于变大；如果S是一个负数，那么后一部分的观测值相比之前的观测值会趋向于变小。

??5.如果n $\le$ 10，依据Gilbert (1987, page 209, Section 16.4.1)中所描述的程序，接下来要在概率表 (Gilbert 1987, Table A18, page 272) 中查找S。如果此概率小于$\alpha$(认为没有趋势时的截止概率)，那就拒绝零假设，认为趋势存在。如果在概率表中找不到n(存在结数据——tied data values——会发生此情况)，就用表中远离0的下一个值。比如S=12，如果概率表中没有S=12，那么就用S=13来处理也是一样的。
??如果n > 10，则依以下步骤6-10来判断有无趋势。这里遵循的是Gilbert (1987, page 211, Section 16.4.2)中的程序。

??6.计算S的方差如下：
??$\mathbf{V}\mathbf{A}\mathbf{R}(\mathbf{S})=\frac{1}{18}[\mathbf{n}(\mathbf{n}?1)(2\mathbf{n}+5)?{\sum }_{\mathbf{p}?1}^{\mathbf{g}}{\mathbf{t}}_{\mathbf{p}}({\mathbf{t}}_{\mathbf{p}}?1)(2{\mathbf{t}}_{\mathbf{p}}+5)](2)$
??其中g是结组（tied groups）的数量，$t_p$是第p组的观测值的数量。例如：在观测值的时间序列{23, 24, 29, 6, 29, 24, 24, 29, 23}中有g = 3个结组，相应地，对于结值(tiied value)23有$t_1$= 2、结值24有$t_2$= 3、结值29有$t_3$= 3。当因为有相等值或未检测到而出现结时，$\mathbf{V}\mathbf{A}\mathbf{R}(\mathbf{S})$可以通过Helsel (2005, p. 191)中的结修正方法来调整。

??7.计算MK检验统计量$Z_{MK}$，如下：

Z M K = { S ? 1 V A R ( S )   , S > 0   0    , S = 0   S + 1 V A R ( S ) , S < 0 ( 3 ) \quad \quad \quad \quad \boldsymbol{Z_{MK}} =

?????????S?1VAR(S)√ ,S>0 0  ,S=0 S+1VAR(S)√,S<0{S?1VAR(S) ,S>0 0  ,S=0 S+1VAR(S),S<0

\quad\quad\quad\boldsymbol{(3)}ZMK=???????VAR(S)S?1 ,S>0 0  ,S=0 VAR(S)S+1,S<0(3)
??正（负）的Z M K \boldsymbol{Z_{MK}}ZMK表明数据随着时间有增大（减小）的趋势。

??8.设想我们要测试零假设
??${\mathbf{H}}_0$：没有单调趋势
????对比替代假设
??${\mathbf{H}}_{\mathbf{a}}$：有单调增趋势
??其1型错误率为$\alpha$，$0<\alpha <0.5$（注意$\alpha$是MK检验错误地拒绝了零假设时可容忍的概率——即MK检验拒绝了零假设是错误地，但这个事情发生概率是$\alpha$，我们可以容忍这个错误）。如果$Z_{MK}\ge Z_{1?\alpha }$，就拒绝零假设$H_0$，接受替代假设$H_a$，其中$Z_{1?\alpha }$是标准正态分布的$100(1?\alpha {)}^{th}$百分位——percentile(不是很懂他说的这些是什么，需要看一下参考文献)。这些百分位在许多统计书(比如 Gilbert 1987, Table A1, page 254)和统计软件包中都有提供。标准正态分布表——百度文库

??9.测试上面的$H_0$与
??${\mathbf{H}}_{\mathbf{a}}$：无单调递减趋势
??其1型错误率为$\alpha$，$0<\alpha <0.5$，如果$Z_{MK}\le ?Z_{1?\alpha }$，就拒绝零假设$H_0$，接受替代假设$H_a$

??10.测试上面的$H_0$与
??${\mathbf{H}}_{\mathbf{a}}$：有单调递增或递减趋势
??其1型错误率为$\alpha$，$0<\alpha <0.5$，如果$|Z_{MK}|\ge Z_{1?\frac{\alpha }{2}}$，就拒绝零假设$H_0$，接受替代假设$H_a$，其中竖线代表绝对值。

缺失数据

??假设时间序列中会有一些缺失数据。比如，数据在每月的第一天被搜集，但是三月一日和七月一日的数据丢失了。在这种情况下$VSP$会以更小的数据集，以通常所用的方法来计算MK检验，适当地减小n的值。
??
未完待续。。。