如何用Python实现MK突变检验方法？

本文共计902个文字，预计阅读时间需要4分钟。

目录需求原理工具和语言代码过程定义函数导入变量形成突变检验图总结需求已知年份和历年最大冻土深度，计算最大冻土深度突变检验。 原理原理+工具+语言+python+jupyter+note代码过程定义函数导入变量形成突变检验图总结

目录

• 需求
• 原理
• 工具和语言
• 代码过程
• 定义函数
• 导入变量 ，形成突变检验图
• 总结

需求

已知年份和历年最大冻土深度，计算最大冻土深度Mk突变检验。

原理

工具和语言

• python
• jupter notebook

代码过程

定义函数

def mktest(inputdata): import numpy as np inputdata = np.array(inputdata) n=inputdata.shape[0] Sk = np.zeros(n) UFk = np.zeros(n) r = 0 for i in range(1,n): for j in range(i): if inputdata[i] > inputdata[j]: r = r+1 Sk[i] = r E = (i+1)*i/4 Var = (i+1)*i*(2*(i+1)+5)/72 UFk[i] = (Sk[i] - E)/np.sqrt(Var) Sk2 = np.zeros(n) UBk = np.zeros(n) inputdataT = inputdata[::-1] r = 0 for i in range(1,n): for j in range(i): if inputdataT[i] > inputdataT[j]: r = r+1 Sk2[i] = r E = (i+1)*(i/4) Var = (i+1)*i*(2*(i+1)+5)/72 UBk[i] = -(Sk2[i] - E)/np.sqrt(Var) UBk2 = UBk[::-1] return UFk, UBk2 定义函数计算变量 ```python def mktest(inputdata): import numpy as np inputdata = np.array(inputdata) n=inputdata.shape[0] s = 0 Sk = np.zeros(n) UFk = np.zeros(n) for i in range(1,n): for j in range(i): if inputdata[i] > inputdata[j]: s = s+1 else: s = s+0 Sk[i] = s E = (i+1)*(i/4) Var = (i+1)*i*(2*(i+1)+5)/72 UFk[i] = (Sk[i] - E)/np.sqrt(Var) Sk2 = np.zeros(n) UBk = np.zeros(n) s = 0 inputdataT = inputdata[::-1] for i in range(1,n): for j in range(i): if inputdataT[i] > inputdataT[j]: s = s+1 else: s = s+0 Sk2[i] = s E = (i+1)*(i/4) Var = (i+1)*i*(2*(i+1)+5)/72 UBk[i] = -(Sk2[i] - E)/np.sqrt(Var) UBk2 = UBk[::-1] return UFk, UBk2

导入变量 ，形成突变检验图

import matplotlib.dates as mdates #處理日期 import matplotlib.pyplot as plt import numpy as np from pylab import mpl from matplotlib.pyplot import MultipleLocator mpl.rcParams['font.sans-serif'] = ['SimHei'] #防止标题出现乱码。 plt.rcParams['axes.unicode_minus'] = False #防止出现图上的负数为方框。 # y值和x值 分别输入六个站点的最大冻土深度值，将值以列表的方式导入 a = [150,150,114,109,96,95,83,76,109,80,115,80,94,86,133,91,110,116,114,128,172,172, 162,121,175,151,110,92,116,156,134,110,89,97,109,157,153,105,76,87,122,78,97,93,141,162, 123,133,161,128,138,104,133,102,140,109,118,86,126,92,121,149,116] #这个部分值可以替换成为要检验的气温、水文等值 x_values=list(range(1961,2022)) uf,ub = mktest(a) plt.figure(figsize=(8,4)) #图片的大小 plt.plot(uf,'r',label='UFk') plt.plot(ub,'b',label='UBk') plt.xticks([0,5,10,15,20,25,30,35,40,45,50,55,60],['1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020',]) #将默认的x轴数值替换为年份的X轴，默认是0-61，一共62个值，代表X轴内容。 # 0.01显著性检验 plt.legend() plt.axhline(1.96) plt.axhline(-1.96) #设置图片的标签（标题） plt.title("富蕴点最大冻土深度突变检验结果")#x轴上的名字 plt.xlabel("年份（1960年-2022年）")#x轴上的名字 plt.ylabel("突变值波动参数")#y轴上的名字 plt.grid() #形成网格线输出 x_major_locator=MultipleLocator(5) plt.show()

最后成图以后的样子。

总结

以上为个人经验，希望能给大家一个参考，也希望大家多多支持自由互联。