Python实现MK突变检验算法

Anonymous888

MK突变检验
根据《现在气候统计诊断与预测计数》（第2版）P63中算法编写Python代码，并绘图。测试数据集是P65中数据。

# @Author        DFY
# @Time                2020-3-17

import matplotlib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
matplotlib.rcParams['xtick.direction'] = 'in'  #in, out, inout
matplotlib.rcParams['ytick.direction'] = 'in'

def get_uf(data):
    index = np.linspace(1,len(data), len(data), True, dtype=int)
    data = pd.Series(data, index=range(1,92))
    s = pd.Series(0, index=index, dtype=np.float)
    uf = pd.Series(0, index=index, dtype=np.float)
    rsum = 0
    for k in index:
        for j in range(1,k):
            if data[k] > data[j]:
                rsum = rsum + 1
        s[k] = rsum
        if k > 1:
            e = (k*(k-1))/4
            var = (k*(k-1)*(2*k+5))/72
            uf[k]= (s[k] - e)/np.sqrt(var)
    return uf

def get_ub(data):
    data = data[::-1]
    uf_temp = get_uf(data)
    ub = uf_temp * -1
    ub = ub[::-1]
    return ub

#准备数据
datas = np.loadtxt("data.txt")
data = datas[:, 1]
year = datas[:, 0]
year = year.astype(np.int)

uf = get_uf(data)
ub = get_ub(data)
line0 = pd.Series(0, index=year)
line1 = pd.Series(1.96, index=year)

#绘图
fig = plt.figure(figsize=(16,12), dpi=80)
ax = fig.add_subplot()
print(ax)
plt.ylim(-4,8)
plt.xlim(1900,1990)
x_major_locator = MultipleLocator(10)
ax.xaxis.set_major_locator(x_major_locator)
x_minor_locator = MultipleLocator(1)
ax.xaxis.set_minor_locator(x_minor_locator)

ax.plot(year,uf, c="black", ls="-", label="uf")
ax.plot(year,ub, c="black", ls="--", label="ub")
ax.plot(year, line0, c="black")
ax.plot(year, line1,c="black")
ax.plot(year, line1*-1, c="black")
ax.legend(ncol=2, framealpha=0)
plt.savefig("mk_fig.png", bbox_inches = 'tight')
plt.show()

123肆

666666  膜拜方有

Anonymous888

123肆 发表于 2018-7-31 09:50
666666  膜拜方有

完了完了， 被发现啦！！！

寻路人

不知道出图和Matlab哪个爽？

943314446

挺好的，谢谢

地表植被

楼主，为什么同样的程序，我报错了：

东方的一条虫

for i in range(1,l):
8.        r = 0
9.        for j in range(i):
10.            if(data > data[j]):
11.                r = r + 1
12.        s = s[i-1] + r
13.        index = i + 1
14.        e = (index * (index-1)) / 4
15.        v = index * (index - 1) * (2 * index + 5) / 72
16.        uf = (s - e) / np.sqrt(v)
17.        uf = round(uf,2)
代码有点问题
第十行 应该是data>data[j](第一个data加上索引i，可能是论坛编辑器自动删除了，因为我这里也加不上去)
第十二行会提示错误：标量变量索引无效，应该是s[i-1]=s[i-1]+r
第十三行index=i+1的话，e的计算应该是index*(index+1)/4

wqry2580

楼主，本人完全的小白，想了解用fortran得到的结果是一致的么？还有这段代码能加上作图的吗？

Dannyqiu

想请教一下楼主，2和-2两条参考线是什么意义呢？

内马尔

Dannyqiu 发表于 2020-5-27 15:12
想请教一下楼主，2和-2两条参考线是什么意义呢？

好像是95%的显著水平？