MK趋势检验和MK突变检验

原理

这两个B站教程讲解的很清楚，而且还有实例演示
B站教程：MK趋势检验
B站教程：MK突变检验
MK检验反应的应该是数据的整体趋势，而不是某一处时间节点的趋势

代码

代码中的测试数据就是B站教程中的测试数据

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


class MK(object):
    def __init__(self, s: list, p=0.05) -> None:
        self.s = s
        self.p = p

    def mkQs(self):
        """mk趋势检验
        Return: return_description
        1 -- 有递增趋势
        -1 -- 有递减趋势
        0 -- 没有明显趋势
        """
        s = self.s
        p = self.p
        S = 0
        n = len(s)
        dsL = []
        for i in range(0, n-1):
            si = s[i]
            for j in range(i+1, n):
                sj = s[j]
                ds = sj-si
                if ds < 0:
                    S -= 1
                elif ds > 0:
                    S += 1
                dsL.append((sj-si)/(j-i))
        sen = np.median(dsL)
        sDict = dict()
        for v in s:
            if v in sDict:
                sDict[v] += 1
            else:
                sDict[v] = 1
        vS1 = n*(n-1)*(2*n+5)
        vS2 = 0
        for v, vn in sDict.items():
            if vn == 1:
                continue
            vS2 += vn*(vn-1)*(2*vn+5)
        vS = (vS1-vS2)/18
        if S > 0:
            Z = (S-1)/np.sqrt(vS)
        elif S == 0:
            Z = 0
        else:
            Z = (S+1)/np.sqrt(vS)
        print(S, vS, Z)
        pV = norm.ppf([1-p/2])
        if np.abs(Z) > pV:
            if sen < 0:
                return -1
            else:
                return 1
        else:
            return 0

    def mkTb(self):
        """Mk突变检验
        Keyword arguments:
        argument -- description
        s -- 时序列表
        p -- 显著水平
        Return: return_description
        UFkList 正序
        UBkList 逆序
        """
        UFkList = [0]
        s = self.s
        n = len(s)
        # 顺序计算
        for k in range(2, n+1):
            sk = 0
            for i in range(0, k):
                si = s[i]
                for j in range(0, i+1):
                    sj = s[j]
                    if si > sj:
                        sk += 1
            usk = k*(k-1)/4
            vsk = k*(k-1)*(2*k+5)/72
            UFk = (sk-usk)/np.sqrt(vsk)
            UFkList.append(UFk)
        # 逆序计算
        s.reverse()
        UBkList = [0]
        for k in range(2, n+1):
            sk = 0
            for i in range(0, k):
                si = s[i]
                for j in range(0, i+1):
                    sj = s[j]
                    if si > sj:
                        sk += 1
            usk = k*(k-1)/4
            vsk = k*(k-1)*(2*k+5)/72
            UBk = (sk-usk)/np.sqrt(vsk)
            UBkList.append(-UBk)
        UBkList.reverse()
        return UFkList, UBkList


if __name__ == '__main__':
    data = [6.8, 5.9, 5.7, 5.5, 5.5, 4.5, 4, 5.1, 4.5, 4.5, 3.3, 4.8]
    mk = MK(data)
    # 趋势性检验
    flag = mk.mkQs()
    if flag == -1:
        print("下降趋势！")
    elif flag == 0:
        print("无明显趋势！")
    else:
        print("上升趋势！")
    # 突变型检验
    UFks, UBks = mk.mkTb()
    plt.figure(figsize=(12, 10))
    plt.plot(UFks, label='UFK')
    plt.plot(UBks, label='UBK')
    # 显著水平p=0.05
    pV = norm.ppf([1-0.05/2])
    plt.hlines([pV, -pV], xmin=0, xmax=12)
    plt.legend()
    plt.show()
 

检验结果

和教程中的完全一致