气象数据分析之突变检验及python的实现：MK突变、Pettitt方法、滑动T检验

---

前言：什么是突变？

常见的气候突变是把它定义为气候从一个平均值到另 一个平均值的急剧变化, 它表现为气候变化的不连续性（符淙斌，1992）。

下图总结了四种常见的突变：
(a)均值突变：从一个均值到另一个均值的变化，表现气候变化的不连续性
(b)变率突变：平均值没有变但是方差变了
©跷跷板突变
(d)转折突变：某一 时段持续减少 ( 增加 ) , 然后突然在某点开始持续增加 (减少 )

检验突变的方法有很多，但是啊每种方法都有优缺，有可能不同方法的检验结果不同，所以建议使用多种方法进行比较。
另外，要指定严格的显著性水平进行检验。

本文介绍几种常用的方法，内容均来自《现代气候统计诊断与预测技术(第二版）》（魏凤英 著）。

1. MK突变分析

Mann-Kendall法是一种非参数统计检验方法，该类型方法亦称为五分部检验，其优点是不需要样本遵从一定的分布，也不受到少数异常值的干扰，更实用于类型变量和顺序变量，计算也比较简便。但是不适用于检测有多个突变点的序列。


例图如下，置信区间内（红色虚线内）的交点就是突变点，正值代表增长，负值反之

from scipy import stats
import numpy as np
from matplotlib import pyplot as plt 
 
def sk(data):
    n=len(data)
    Sk     = [0]
    UFk    = [0]
    s      =  0
    E      = [0]
    Var    = [0]
    for i in range(1,n):
        for j in range(i):
            if data[i] > data[j]:
                s = s+1
            else:
                s = s+0
        Sk.append(s)
        E.append((i+1)*(i+2)/4 )                     # Sk[i]的均值
        Var.append((i+1)*i*(2*(i+1)+5)/72 )            # Sk[i]的方差
        UFk.append((Sk[i]-E[i])/np.sqrt(Var[i]))
    UFk=np.array(UFk)
    return UFk

#a为置信度
def MK(data,a):
    ufk=sk(data)          #顺序列
    ubk1=sk(data[::-1])   #逆序列
    ubk=-ubk1[::-1]        #逆转逆序列
    
    #输出突变点的位置
    p=[]
    u=ufk-ubk
    for i in range(1,len(ufk)):
        if u[i-1]*u[i]<0:
            p.append(i)            
    if p:
        print("突变点位置：",p)
    else:
        print("未检测到突变点")
    
    #画图
    conf_intveral = stats.norm.interval(a, loc=0, scale=1)   #获取置信区间

    plt.figure(figsize=(10,5))
    plt.plot(range(len(data)),ufk,label = 'UFk',color = 'r')
    plt.plot(range(len(data)),ubk,label = 'UBk',color = 'b')
    plt.ylabel('UFk-UBk',fontsize=25)
    x_lim = plt.xlim()
    plt.ylim([-6,7])
    plt.plot(x_lim,[conf_intveral[0],conf_intveral[0]],'m--',color='r',label='95%显著区间')
    plt.plot(x_lim,[conf_intveral[1],conf_intveral[1]],'m--',color='r')
    plt.axhline(0,ls="--",c="k")
    plt.legend(loc='upper center',frameon=False,ncol=3,fontsize=20) # 图例
    plt.xticks(fontsize=25)
    plt.yticks(fontsize=25)
    plt.show()

#输入数据和置信度即可
MK(data,0.95)

2. Pettitt方法

是一种与MK方法相似的非参数检验方法。


例图如下，“突变点位置:32, 显著”

def Pettitt(data):
    data = np.array(data)
    n = len(data)
    sk     = [0]
    
    for i in range(1,n):
        s = 0
        for j in range(i):
            if data[i] > data[j]:
                s = s+1
            if data[i] < data[j]:
                s = s-1
            else:
                s = s+0
        sk.append(s)
        
    k = np.max(sk)
    kt = sk.index(max(sk))    
    print(k,kt)
    p = 2 * np.exp((-6 * (k**2))/(n**3 + n**2))
    
    if p <= 0.05:
        a = '显著'
    else:
        a = '不显著'
    print('突变点位置:%s, %s'%(kt,a))
    
    #画图
    plt.plot(range(len(data)),data)
    plt.plot([0,kt],[np.mean(data[0:kt]),np.mean(data[0:kt])],'m--',color='r')
    plt.plot([kt,len(data)],[np.mean(data[kt::]),np.mean(data[kt::])],'m--',color='r')
    plt.axvline(x=kt,ls="--",c="g")
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    
    return k,kt #,Pettitt_result

Pettitt(data)

3. 滑动T检验(Moving T test , MTT)

滑动t检验是通过考察两组样本平均值的差异是否显著来检验突变。
基本思想是把一气候序列中两段子序列均值有无显著差异看作来自两个总体均值有无显著差异的问题来检验。如果两段子序列的均值差异超过了一定的显著性水平，则可以认为有突变发生。

顺便叨叨一句怎么查t分布表
从左一列找到自由度，第一行表头表示的就是显著性，比如我选的步长为5、显著性95%，则自由度v=5+5-2=8，8和0.05对应的就是2.31。

例图如下，超过显著区间的点就是突变点，正（负）值代表增长（下降）

import numpy as np
from matplotlib import pyplot as plt 

def MTT(data, step):
    n = len(data)
    v=step+step-2 #自由度
    t=np.zeros((n-step-step+1))
    ss=np.sqrt(1/step+1/step)  
    
    ttest = 2.31  #step=5,alpha=0.05，这个需要根据需要查表做改动
    
    for i in range(len(t)):
        n1=data[i:i+step]
        n2=data[i+step:i+step+step]
        x1=np.mean(n1)#平均值
        x2=np.mean(n2)
        s1=np.std(n1)#方差
        s2=np.std(n2)
        s=np.sqrt((step*s1*s1+step*s2*s2)/v)
        t[i]=(x1-x2)/(s*ss) 
    
    plt.plot(t,label="滑动T值(step=%s)"%step)
    plt.axhline(0,ls="--",c="k")
    plt.axhline(ttest,ls="--",c="r",label='95%显著区间')
    plt.axhline(-ttest,ls="--",c="r")
    plt.legend(loc='upper center',frameon=False,ncol=2,fontsize=14) # 图例
    
    return t

写得乱七八糟的先更新到这里啦还有一些方法我日后随缘加上~
祝大家磕盐顺利、身心健康！
有什么错误的地方欢迎在评论区批评指正~