python如何实现MK突变检验方法,代码复制修改可用

需求

已知年份和历年最大冻土深度，计算最大冻土深度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()

最后成图以后的样子。

总结