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')
借鉴其他大佬代码,稍作修改,可以运行