零门槛学习——Mannkendall趋势检测算法

背景介绍

厉害的算法有很多,但是找到最合适数据情况的算法,才是前期工作中最重要的一步。对于单条时序数据的趋势判断,本文将介绍Mannkendall算法,作为一种非参数统计检验方法,

通过对数据之间的符号差值进行统计分析,可以识别其显著趋势。

使用MK 算法检验时序数据大致趋势时,趋势分为无明显趋势(稳定)、趋势上升、趋势下降。该检验比较样本数据的相对量级,而不是数据值本身,因此基础数据不需要遵循特定分布,也不要求变化趋势是线性的。使用比较每一对数据点来计算总体的统计量S,并将其标准化为Z 值,从而判断时间序列数据是否存在显著的趋势。

算法原理

1.假设检验

前提: H 0 H_0 H0为真,或者在拒绝 H 0 H_0 H0,接受 H 1 H_1 H1时,数据的统计值需要达到一定的置信度

零假设 H 0 H_0 H0:没有单调趋势

H 1 H_1 H1:有单调趋势

2.根据给定的时间序列,得到对应数据{ x 1 , x 2 , . . . x n x_1,x_2,...x_n x1,x2,...xn}

对应这个算法的计算数据总量为 n ( n − 1 ) 2 \frac{n(n-1)}{2} 2n(n1)

3.定义差值函数

sgn ⁡ ( x j − x i ) = { + 1 , x j − x i > 0 0 , x j − x i = 0 − 1 , x j − x i < 0 \operatorname{sgn}\left(x_{j}-x_{i}\right)=\left\{\begin{aligned}+1, & x_{j}-x_{i}>0 \\0, & x_{j}-x_{i}=0 \\-1, & x_{j}-x_{i}<0\end{aligned}\right. sgn(xjxi)= +1,0,1,xjxi>0xjxi=0xjxi<0

4.计算统计量

统计量 S S S 是所有符号差值的和:
S = ∑ i = 1 n − 1 ∑ j = i + 1 n sgn ⁡ ( x j − x i ) S=\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \operatorname{sgn}\left(x_{j}-x_{i}\right) S=i=1n1j=i+1nsgn(xjxi)

S S S的值很大或很小时,表明存在显著的趋势。

趋势的强度与 S S S 的大小成正比, S S S的绝对值越大,真实增加或减少趋势的证据越强

如果 S S S是一个正数,那么后一部分的观测值相比之前的观测值会趋向于变大(上涨),

如果 S S S是一个负数,那么后一部分的观测值相比之前的观测值会趋向于变小(下降)。

计算完 S S S其实已经可以得到一个趋势的判断了,为什么还要继续下面的计算?

  • 首先算法执行下面步骤的前提是样本量n>10,满足。(如果样本量<10,应在概率表中查找S对应的p值,然后比较α来确实是否有趋势)
  • 计算期望方差是为了得到z值,然后进一步计算P值。

P值的作用是?

MK 检验的计算概率(p 值)表示任何观察到的趋势纯属偶然发生的概率(考虑到数据集的可变性和样本大小)。

所以后面的所有计算量其实都是为了最终得到P值,确定算出来的结果不是偶然和意外。

  • 对此,科学的解释是,显著性水平 0.05(即 95% 置信度)用于检验数据中没有趋势的零假设。显著性水平是测试错误地检测到没有趋势的概率。MK 检验的结果可能是显著增加或减少的趋势,也可能是不显著的结果(不能拒绝无趋势的零假设)。
5.期望值和方差计算

在没有趋势的情况下,期望值 E ( S ) E(S) E(S)和方差 V a r ( S ) Var(S) Var(S)的计算公式如下:

E ( S ) = ∑ i = 1 n − 1 ∑ j = i + 1 n sgn ⁡ ( x j − x i ) E(S)=\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \operatorname{sgn}\left(x_{j}-x_{i}\right) E(S)=i=1n1j=i+1nsgn(xjxi)

对于一个没有趋势的时间序列,S 的期望值。这是因为在随机数据中,增加和减少的对数应该大致相等。

对于方差,计算的原始公式如下:
$Var(S)=\frac{n(n-1)(2n+5)}{18} $

6.标准化检验统计量Z 的计算

为了确定 Z Z Z 的显著性,我们计算标准化检验统计量(这里是统计学定义公式):

Z = { S − 1 Var ⁡ ( S ) , S > 0 0 S = 0 S + 1 Var ⁡ ( S ) , S < 0 Z=\left\{\begin{array}{cc}\frac{S-1}{\sqrt{\operatorname{Var}(S)}}, & S>0 \\0 & S=0 \\\frac{S+1}{\sqrt{\operatorname{Var}(S)}}, & S<0\end{array}\right. Z= Var(S) S1,0Var(S) S+1,S>0S=0S<0

根据 Z Z Z 的值,我们可以使用正态分布来评估趋势的显著性。

7.假设检验判断

计算置信度 α \alpha α p p p α \alpha α可以用来决定是否拒绝原假设。当时,即表示有5%的概率在原假设为真的情况下错误地拒绝原假设。

是基于样本数据计算出来的概率,用于评估观察到的结果在原假设为真时出现的概率。

P −  value  = 2 ( 1 − cd ⁡ f ( ∣ Z M K ∣ ) ) P-\text { value }=2\left(1-\operatorname{cd} f\left(\left|Z_{M K}\right|\right)\right) P value =2(1cdf(ZMK))

(由于是双边检测,就是它不管是上升还是下降都要看,所以要*2来跟, α \alpha α比较)

如果 p 值 < α p值<\alpha p<α,则拒绝原假设,即时间序列中存在显著趋势。

如果 p 值 > α p值>\alpha p>α,则接受原假设,即时间序列中不存在显著趋势。

该检验的零假设是没有趋势,备择假设是双侧检验中有趋势,或单侧检验中有上升趋势(或下降趋势)。

举个例子

比如说有10个数,用第二个减去第一个,第三个减去第一个,第四个减去第一个,这样分别减,都会得到他们对应的正负不同的结果(对应有增有降),只保留这些结果为 ± 1 \pm1 ±1或为0,一直进行这个过程到第10个数,最后就会得到9个数,他们分别代表着各自对应着第一天的趋势。然后把这几个正负求和的结果就是后面的数字对于第一天的判断结果,然后第二天,第三天也一样。得到他们每天的数字判断结果之后,进行s、z的统计量计算

假设我们有以下一组时间序列数据,表示某地区5年的年平均气温:

模拟数据=[12,14,13,16,15]

1.计算S统计量

我们需要比较每一对数据点。具体操作如下:

  • 比较第1年和后面的年份:

14 > 12,计数 +1

13 > 12,计数 +1

16 > 12,计数 +1

15 > 12,计数 +1

结果:4

  • 比较第2年和后面的年份:

13 < 14,计数 -1

16 > 14,计数 +1

15 > 14,计数 +1

结果:1

  • 比较第3年和后面的年份:

16 > 13,计数 +1

15 > 13,计数 +1

结果:2

  • 比较第4年和第5年:

15 < 16,计数 -1

结果:-1

把所有结果相加,得到 S: S=4+1+2−1=6

2.计算期望值和方差

E ( S ) = 0 E(S)=0 E(S)=0

$Var(S)=\frac{n(n-1)(2n+5)}{18} $

对于n=5,方差的计算如下:
$Var(S)=\frac{5 \times 4 \times 15}{18}=\frac{300}{18}=16.67 $

因此,对于特定值 n = 5 n=5 n=5,方差$Var(S)=16.67 $

3.计算Z值

给定随机变量 S S S的正负性,计算标准化值 Z Z Z
Z = S − 1 V a r ( S ) Z=\frac{S-1}{\sqrt{Var(S)}} Z=Var(S) S1

S = 6 S=6 S=6且$Var(S)=16.67 $时,计算步骤如下:
Z = 6 − 1 16.67 = 5 16.67 ≈ 5 4.08 ≈ 1.22 Z=\frac{6-1}{\sqrt{16.67}}=\frac{5}{\sqrt{16.67}}\approx\frac{5}{4.08}\approx1.22 Z=16.67 61=16.67 54.0851.22

因此标准化值 Z Z Z约为1.22

4.确定显著性

通过查找标准正态分布表,Z值为1.22对应的p值约为0.22。一般我们选择的显著性水平是0.05或者0.01,对应的临界值分别是1.96和2.58。

因为1.22 < 1.96,所以在0.05显著性水平下,我们不能拒绝无趋势假设。也就是说,这组数据没有显著的上升或下降趋势。

总结

通过这个简单的例子,可以看到MK-test如何使用比较每一对数据点来计算总体的统计量S,并将其标准化为Z值,从而判断时间序列数据是否存在显著的趋势。在这个例子中,尽管S值是正的,但Z值不够大,因此我们不能认为数据存在显著的上升趋势。

代码实现

使用python进行mk-test非常容易,因为已经有开源的package可以使用了,链接如下

https://github.com/mmhs013/pyMannKendall/blob/master/Examples/Example_pyMannKendall.ipynb

pip install pymannkendall

然后就可以调用了

import pymannkendall as mk

# 给定数据
data = [12, 14, 13, 16, 15]
# 进行 Mann-Kendall 趋势检验

result = mk.original_test(data)

# 输出结果

print(f"S: {result.S}")

print(f"Var(S): {result.VarS}")

print(f"Z: {result.Z}")

print(f"P-value: {result.p}")

print(f"Trend: {result.trend}")

print(f"H: {result.h}")

· S: 统计量S。

· Var(S): S的方差。

· Z: Z值。

· P-value: P值。

· Trend: 趋势方向(‘increasing’ 或 ‘decreasing’ 或 ‘no trend’)。

· H: 趋势的显著性(True表示趋势显著,False表示趋势不显著)。

改良算法

由于在此类任务中,对于一些特殊时间序列: [3,3,5,6,6,6,8],其中有以下平局组:

  • 平局组1:值为3的两个观测值 t 1 = 2 t_1=2 t1=2
  • 平局组2:值为6的三个观测值 t 2 = 3 t_2=3 t2=3

使用原始的 Mann-Kendall (MK) 算法在处理平局组 (ties) 时会遇到以下几个问题:

1. 偏差问题

平局组在序列中会引入偏差,影响统计结果的准确性。原始 MK 算法假设所有观测值都是独立的和连续的,当存在平局组时,这一假设被打破,导致统计量 SSS 的期望值和方差的计算出现偏差。

2. 统计量 S 的计算

在原始 MK 算法中,统计量 S 通过比较所有可能的数字对来计算:
S = ∑ i = 1 n − 1 ∑ j = i + 1 n sgn ⁡ ( x j − x i ) S=\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \operatorname{sgn}\left(x_{j}-x_{i}\right) S=i=1n1j=i+1nsgn(xjxi)

其中 s g n ( x j − x i ) sgn(x_j-x_i) sgn(xjxi)表示符号函数,当 x j > x i x_j>x_i xj>xi时,值为1;当 x j < x i x_j < x_i xj<xi时,值为-1;

x j = x i x_j=x_i xj=xi时,值为0。当存在平局组时,符号函数会产生大量的0,这会对 $S $的统计值产生显著影响。

3. 方差的计算

原始 MK 算法中方差的计算公式是基于没有平局组的假设:
$Var(S)=\frac{n(n-1)(2n+5)}{18} $

当存在平局组时,上述公式不再适用。为了修正这一问题,需要引入平局组的调整项:
V a r ( S ) = 1 18 [ n ( n − 1 ) ( 2 n + 5 ) − ∑ p − 1 g t p ( t p − 1 ) ( 2 t p + 5 ) ] Var(S)=\frac{1}{18}\left[n(n-1)(2n+5)-\sum_{p-1}^{g} t_p\left(t_{p}-1\right)\left(2 t_p+5\right)\right] Var(S)=181[n(n1)(2n+5)p1gtp(tp1)(2tp+5)]

其中 g g g是平局组的数量, t p t_p tp是第 p p p个平局组的大小。这一调整项可以有效地修正平局组对方差的影响。

4. p 值计算

在存在平局组的情况下,原始 MK 算法计算的统计量 Z 的正态近似可能不准确,导致 p 值计算出现偏差。调整后的算法通过修正方差来改善这一问题,但仍需要注意在实际应用中可能存在的其他影响因素。

除此之外,还有其他针对原始mk检验的改良算法,可以稍作了解:

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

Yue and Wang 改进的 MK 检验 ( yue_wang_modification_test ):这也是Yue, S., & Wang, CY (2004)提出的一种考虑序列自相关的方差校正方法。用户还可以设置所需的计算显著 n 滞后。

优势

个人使用感受:计算复杂度是 O ( n 2 ) O(n_2) O(n2),对于单条数据,总体上来说mk的判断结果是比较让人满意的,尤其是对于长期数据,整体的判断会更加直观。计算结果只需要通过 S , Z S,Z SZ就可以快速知道,而对于一些特殊时序数据类型,如本文所提到的平局组,只需要针对性进行调整公式,就可以更加贴近实际情况。

最后总结一下算法优势如下

  1. 非参数方法:无需假设数据的特定分布形式,适用于各种类型的数据。对于店铺商品日销量数据,这种灵活性特别重要,因为销量数据往往不符合正态分布等常见分布假设。

  2. 稳健性:因为实际销售数据可能会受促销活动、季节变化等多种因素影响,存在波动和异常值。该算法对异常值和缺失数据不敏感,能够提供相对可靠的检测结果。

  3. 运算开销低:计算速度快,且计算效率高。

4.有现成的开源包,方便调用

参考链接:

https://www.stats.uwo.ca/faculty/mcleod/2003/DBeirness/MannKendall.pdf

https://real-statistics.com/time-series-analysis/time-series-miscellaneous/mann-kendall-test/

https://github.com/mmhs013/pyMannKendall

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值