Python实现 MK检验示例代码

MK检验：时间序列进行检测，并找出突变点，本文参考网上的matlwww.cppcns.comab程序改写为python代码如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#读取时间序列数据
data = pd.read_csv('')
#定义时间和径流数据列
x = list(range(len(data)))
y = data.to_list
#获取样本数据
n = len(y)

#正序计算
#定义累计量序列ShTuLvBvk，长度n，初始值为0
Sk = np.zeros(n)
UFk = np.zeros(n)
#定义Sk序列元素s
s = 0

#i从2开始，根据统计量UFk公式，i=1时，Sk（1）、E（1）、Var（1）均为0，此时UFk无意义，因此公式中，令UFk（1） = 0
for i in range(2,n):
    for j in range(1,i):
        if y[i]>y[j]:
            s += 1
    Sk[i] = s
    E = i * (i - 1)/4
    Var = i * (i - 1) * (2 * i + 5)/72
    UFk[i] = (Sk[i] - E)/np.sqrt(Var)

#逆序计算，构造逆序列y2，长度为n，初值为0
#定义逆累计量序列Sk2，长度n，初始值为0
#定义逆统计量序列Sk2，长度n，初始值为0
y2 = np.zeros(n)
Sk2 = np.zeros(n)
UBk = np.zeros(n)

#s归零
s = 0
#按时间序列逆转样本
y2 = y[::-1]

# i从2开始，根据统计量UBk公式，i=1时，Sk（1）、E（1）、Var（1）均为0，此时UFk无意义，因此公式中，令UBk（1） = 0
for i in range(2, n):
    for j in range(1, i):
        if y2[i] > y2[j]:
            s += 1
    Sk2[i] = s
    E = i * (i - 1) / 4
  hTuLvBv  Var = i * (i - 1) * (2 * i + 5) / 72
www.cppcns.com    UBk[i] = -(Sk2[i] - E) / np.sqrt(Var)

#逆转逆序列
UBk2 = UBk[::-1]

#画图

plt.figure(figsize=(10,5))
plt.plot(range(1 ,n+1),UFk,label = 'UFk',color = 'orange')
plt.plot(range(1 ,n+1),UBk2,label = 'UBk',color = 'cornflowerblue')hTuLvBv
plt.ylabel('UFk-UBk')
x_lim = plt.xlim()
plt.plot(x_lim,[-1.96,-1.96],'m--',color = 'r')
plt.plot(x_lim, [0,0],'m--')
plt.plot(x_lim,[1.96,1.96],'m--',color = 'r')
plt.show()

到此这篇关于Python实现MK检验的文章就介绍到这了,更多相关PythonMK检验内容请搜索我们以前的文章或继续浏览下面的相关文章希望大家以后多多支持我们！