基于python滑动T检验的实现--结合MK突变确定突变点

滑动T检验

前面做了基于python的mk突变检验发现UF,UB曲线有较多相交点，于是通过查询相关代码，最终基于python实现了滑动T检验方法，用来检验MK趋势检验中突变点的显著性，确定突变点。

由于该结果用于论文，故最终结果图整体配色为黑色。

如果代码有什么不足之处，请大家指正。

原理

python代码

# coding:utf-8
from matplotlib import pyplot as plt
from tqdm import tqdm
import pandas as pd
import numpy as np
import os
import xlwt


# v:自由度
def get_tvalue(v, sig_level):
    t_values = pd.read_excel(r'../t_values.xlsx')
    return t_values[t_values['n'] == v][sig_level].values


# 传入时间序列以及待检验数据
# data:数据， step：步长， sig_level：显著水平
def huaTTest(data, step, sig_level):
    datacount = len(data)
    v = 2*step-2  # 自由度
    t_value = get_tvalue(v, sig_level)  # t检验值
    n1 = step
    n2 = step
    t = np.zeros(len(data))
    c = 1.0 / n1 + 1.0 / n2

    for i in range(step-1, datacount-step):
        data1 = data[i-step+1:i+1]
        data2 = data[i+1:i+step+1]
        # 计算均值
        x1_mean = data1.mean()
        x2_mean = data2.mean()
        # 计算方差
        s1 = data1.var()
        s2 = data2.var()
        sw2 = (n1*s1 + n2*s2)/(n1+n2-2.0)
        t[i - step + 1] = (x1_mean - x2_mean) / np.sqrt(sw2 * c)
    return t, t_value

plt.rcParams['font.family'] = ['SimHei']  # 设置字体
plt.rcParams['axes.unicode_minus'] = False

# 获取数据
path = ur'../单年均值.xlsx'
data_xls = pd.read_excel(path)

# 获取时间列以及casa模拟值
time = np.array(data_xls['year'])
casa = np.array(data_xls['casa'])
datacount = len(casa)

# 计算步长范围
if datacount % 2 == 0:
    steps = range(2, datacount/2)
else:
    steps = range(2, datacount/2+1)

dict = {}  # 存储显著变化的点信息
# 按照可取步长计算
for step in tqdm(steps, ncols=80, desc=u'滑动T检验进度'):
    t, t_value = huaTTest(casa, step, 0.05)

    t = t[t != 0]

    sig_values = []  # 存储显著变化值
    for i in range(len(t)):
        if np.abs(t[i]) > t_value:
            sig_values.append(time[i + step - 1])
    if sig_values:
        # 将结果存入词典
        dict_key = 'step'+str(step)
        dict[dict_key] = sig_values

    fig, ax = plt.subplots(figsize=(9, 5), sharey=True)

    fontdict = {'size': 12, 'color': 'k'}

    x = time[step-1:datacount-step]
    ax.set_xticks(x[::2])
    # ax.set_xlim((time[0], time[datacount-step]))
    ax.plot(x, t, c='k', linewidth=1.5, label=u'滑动T值', zorder=1)
    # 添加显著水平线
    ax.hlines(y=t_value, xmin=np.min(x), xmax=np.max(x), linewidth=1.5, linestyles='--', label=u'0.05显著水平')
    ax.hlines(y=-t_value, xmin=np.min(x), xmax=np.max(x), linewidth=1.5, linestyles='--', label=None)
    # 修改刻度线指向
    ax.tick_params(left=True, bottom=True, direction='in', labelsize=12)
    # 添加坐标轴标签
    ax.set_xlabel(u'年份 Year', fontdict=fontdict)
    ax.set_ylabel('T', fontdict=fontdict)
    # 设置图例
    ax.legend(bbox_to_anchor=(0.95, 0.95), facecolor='w', frameon=False)

    # 保存文件
    basepath = r'../IMG/huaT'
    filename = r'/huaTstep%d.png'%step
    savepath = basepath + filename

    # 检查保存目录是否存在
    if not os.path.exists(basepath):
        print('创建保存目录')
        os.mkdir(basepath)
        print('创建完成！')

    plt.savefig(savepath)

# 写入数据
# 创建workbool和sheet对象
workbook = xlwt.Workbook()
sheet1 = workbook.add_sheet('sheet1', cell_overwrite_ok=True)

# 写入头部信息
for head in range(len(dict.keys())):
    sheet1.write(0, head, label=dict.keys()[head])

# 写入数据
# 获取最大行数
max_row = max([len(dict.values()[i]) for i in range(len(dict.values()))])
for row in range(max_row):
    for col in range(len(dict.keys())):
        if row < len(dict.values()[col]):
            sheet1.write(row+1, col, label=int(dict.values()[col][row]))  # 将np.int32转换为int

# 保存结果
workbook.save(r'..\IMG\huaT\huaTSig.xls')

结果

结果图： step==3

结果图： step==2

相关数据可查看

基于python的滑动T检验及相应数据