1. 程式人生 > 實用技巧 >Python實現 MK檢驗

Python實現 MK檢驗

MK檢驗:時間序列進行檢測,並找出突變點,本文參考網上的matlab程式改寫為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)

#正序計算
#定義累計量序列Sk,長度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 Var = i * (i - 1) * (2 * i + 5) / 72 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') 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()