# 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')
评论0