R语言Mk检验+python滑动T检验

R语言Mk检验+python滑动T检验

MK检验

#####设置文件路径######
setwd("D:\\station\\站点SPEI计算\\spei3\\")
#####原始数据#####
library(openxlsx)
data<-read.xlsx("MKwinter.xlsx",colNames = FALSE)
data<-as.matrix(data)
#####对数据求倒序#####
data_dx<-rev(data)
data_dx<-matrix(data_dx)
n<-length(data)
#######MK检验#######
Q<-matrix(0,1,n)
UF<-matrix(0,1,n)
UB<-matrix(0,1,n)
h<-matrix(0,n,n)
for (i in 1:n) {
  for (j in 1:n) {
    if(data[i,1]>data[j,1]){
      h[i,j]<-1
    }
    else
      h[i,j]<-0
  }
  Q[1,i]<-sum(h[lower.tri(h)])
  UF[1,i]<-(Q[1,i]-(i*(i-1)/4))/sqrt((i*(i-1)*(2*i+5))/72)
}
#####计算UB#####
h<-matrix(0,n,n)
for (i in 1:n) {
  for (j in 1:n) {
    if(data_dx[i,1]>data_dx[j,1]){
      h[i,j]<-1
    }
    else
      h[i,j]<-0
  }
  Q[1,i]<-sum(h[lower.tri(h)])
  UB[1,i]<-(Q[1,i]-(i*(i-1)/4))/sqrt((i*(i-1)*(2*i+5))/72)
}
#####绘图#####
UF[1,1]<-0
UB[1,1]<-0
plot(x=1990:2020,y=UF,ylim = c(-4,4),type = "l",ylab = "",xaxt="n",lwd="2.0",col=2)
lines(x=1990:2020,y=-rev(UB),type = "l",lty=2,lwd="2.0",col=4)
lines(x=1990:2020,y=rep(1.96,31),type = "l",lty=5)
lines(x=1990:2020,y=rep(-1.96,31),type = "l",lty=5)
lines(x=1990:2020,y=rep(0,31),type = "l")
title("冬季M-K检验")
axis(1,1990:2020,1990:2020,las=1)
#####输出突变年份#####
year_mk<-1990:2020
year_point<-year_mk[which((as.numeric(UF)-(-rev(UB)))>0)[1]-1]
print(year_point)

滑动T检验

# coding:utf-8
from tqdm import tqdm
from matplotlib import pyplot as plt
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"D:\python\pygame-master\pygame-master\GluttonousSnake\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 = r'./HTyear.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_ylim(-5, 5)
    # ax.set_xlim((time[0], time[datacount-step]))
    ax.plot(x, t, c='k', linewidth=1.5, label=u'滑动T值', zorder=1,marker ='o')
    # 添加显著水平线
    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'./192021/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=list(dict.keys())[head])

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

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

借鉴其他大佬代码，稍作修改，可以运行