python如何實現(xiàn)MK突變檢驗方法,代碼復制修改可用

更新時間：2023年05月08日 10:23:10 作者：David_wangzw

這篇文章主要介紹了python如何實現(xiàn)MK突變檢驗方法，代碼復制修改可用，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教

需求

已知年份和歷年最大凍土深度，計算最大凍土深度Mk突變檢驗。

原理

請?zhí)砑訄D片描述

工具和語言

python
jupter notebook

代碼過程

定義函數(shù)

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
定義函數(shù)計算變量
```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'] #防止標題出現(xiàn)亂碼。
plt.rcParams['axes.unicode_minus'] = False   #防止出現(xiàn)圖上的負數(shù)為方框。
# 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軸數(shù)值替換為年份的X軸，默認是0-61，一共62個值，代表X軸內(nèi)容。
# 0.01顯著性檢驗
plt.legend()
plt.axhline(1.96)
plt.axhline(-1.96)
#設置圖片的標簽（標題）
plt.title("富蘊點最大凍土深度突變檢驗結(jié)果")#x軸上的名字
plt.xlabel("年份（1960年-2022年）")#x軸上的名字
plt.ylabel("突變值波動參數(shù)")#y軸上的名字
plt.grid() #形成網(wǎng)格線輸出
x_major_locator=MultipleLocator(5)
plt.show()