python如何实现MK突变检验方法,代码复制修改可用
作者:David_wangzw 时间:2022-04-10 13:31:18
需求
已知年份和历年最大冻土深度,计算最大冻土深度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()
最后成图以后的样子。
来源:https://blog.csdn.net/weixin_52753312/article/details/127953122
标签:python,MK,突变检验
![](/images/zang.png)
![](/images/jiucuo.png)
猜你喜欢
python中获得当前目录和上级目录的实现方法
2022-01-07 20:24:30
Python可视化tkinter详解
2022-12-31 06:09:19
![](https://img.aspxhome.com/file/2023/9/115509_0s.png)
Python scikit-learn 做线性回归的示例代码
2022-05-03 11:00:54
![](https://img.aspxhome.com/file/2023/1/97661_0s.png)
Oracle 数据 使用游标
2009-07-02 12:14:00
Python3环境安装Scrapy爬虫框架过程及常见错误
2021-10-19 00:01:05
一个简单的ASP生成HTML分页程序
2009-07-05 18:32:00
Django基于Token的验证使用的实现
2023-06-14 18:43:54
![](https://img.aspxhome.com/file/2023/6/123636_0s.jpg)
一分钟带你掌握Python中pip的安装与使用方法
2021-02-10 10:38:12
![](https://img.aspxhome.com/file/2023/1/83921_0s.jpg)
python实现MD5进行文件去重的示例代码
2021-12-13 02:28:23
![](https://img.aspxhome.com/file/2023/0/97550_0s.png)
div + ajax + 分页函数
2009-10-18 11:28:00
Ubuntu16.04 安装多个python版本的问题及解决方法
2021-05-26 05:27:11
![](https://img.aspxhome.com/file/2023/6/90636_0s.jpg)
python爬虫增加访问量的方法
2021-08-23 06:32:23
闲谈CSS3动画
2010-05-07 12:34:00
![](https://img.aspxhome.com/file/UploadPic/20105/7/t1ztvzxxhdxxxxxxxx-226-58-35s.png)
python实现将html表格转换成CSV文件的方法
2023-08-25 00:48:41
详解用python生成随机数的几种方法
2022-01-24 14:17:42
SQL Server 中死锁产生的原因及解决办法
2008-11-25 11:50:00
微信小程序 云开发模糊查询实现解析
2023-08-24 14:47:57
python3使用requests模块爬取页面内容的实战演练
2022-01-08 18:26:57
![](https://img.aspxhome.com/file/2023/5/90625_0s.png)
python中defaultdict用法实例详解
2022-08-09 17:01:10
![](https://img.aspxhome.com/file/2023/9/109339_0s.png)
用Python实现服务器中只重载被修改的进程的方法
2022-06-21 05:11:38