相空間重構求關聯維數——GP演算法、自相關法求時間延遲tau、最近鄰演算法求嵌入維數m
技術標籤：演算法
相空間重構求關聯維數——GP演算法、自相關法求時間延遲tau、最近鄰演算法求嵌入維數m
GP演算法：
若有一維時間序列為{x1,x2,…,xn},對其進行相空間重構得到高維相空間的一系列向量：
x i ( τ , m ) = ( x i , x i 1 , ⋯ , x i + ( m − 1 ) τ ) {x_i}(\tau ,m) = \left( {{x_i},{x_{i1}}, \cdots ,{x_{i + {{(m - 1)}_\tau }}}} \right) xi(τ,m)=(xi,xi1,⋯,xi+(m−1)τ)
式中：
τ
\tau
重構相空間關聯維數為：
D 2 = lim r → 0 ln c r ln r {D_2} = \mathop {\lim }\limits_{r \to 0} \frac{{\ln {c_r}}}{{\ln r}} D2=r→0limlnrlncr
c
r
=
1
N
2
{c_r} = \frac{1}{{{N^2}}}
式中：j≠k；r為m維超球半徑；H為Heaviside函式。
def GP(imf,tau): #GP演算法求關聯維數
N=2000
if (len(imf) != N):
print('請輸入指定的資料長度！') # N為指定資料長度
return
elif (isinstance(imf, np. 自相關法確定 τ \tau τ：
計算時間序列{x1,x2,…,xn}的自相關函式：
R ( j τ ) = 1 N ∑ R(j\tau )= \frac{1}{{{N}}}\sum R(jτ)=N1∑ x ( i ) x ( i + j τ ) x(i)x(i + j\tau ) x(i)x(i+jτ)
當自相關函式值下降到初始函式值的1- e − 1 {{\rm{e}}^{ - 1}} e−1時。所對應的 τ \tau τ即為時間延遲引數。
# 計算GP演算法的時間延遲引數（自相關法）
def get_tau(imf):
N=2000
if (len(imf) != N):
print('請輸入指定的資料長度！') # N為指定資料長度
return 0
elif (isinstance(imf, np.ndarray) != True):
print('資料格式錯誤！')
return 0
else:
j = 1 # j為固定值
tau_max = 20
Rall = np.zeros(tau_max)
for tau in range(tau_max):
R = 0
for i in range(N - j * tau):
R += imf[i] * imf[i + j * tau]
Rall[tau] = R / (N - j * tau)
for tau in range(tau_max):
if Rall[tau] < (Rall[0] * 0.6321):
break
return tau
假近鄰演算法確定m
對m維相空間每一個向量 X i ( m ) = { x i , x i + τ , ⋯ , x i + ( m − 1 ) τ } {X_{i(m)}} = \left\{ {{x_i},{x_{i + \tau }}, \cdots ,{x_{i + (m - 1)\tau }}} \right\} Xi(m)={xi,xi+τ,⋯,xi+(m−1)τ},i=1,2,…,N,N為向量總數，找出它的最近向量 X j ( m ) X_{j(m)} Xj(m)，計算兩者歐氏距離 R m ( i ) = ∣ ∣ X i ( m ) − X j ( m ) ∣ ∣ {R_{m }}(i) = ||{X_{i(m)}} - {X_{j(m )}}|| Rm(i)=∣∣Xi(m)−Xj(m)∣∣，它們在m+1維空間的距離為：
R m + 1 ( i ) = ∣ ∣ X i ( m + 1 ) − X j ( m + 1 ) ∣ ∣ {R_{m + 1}}(i) = ||{X_{i(m + 1)}} - {X_{j(m + 1)}}|| Rm+1(i)=∣∣Xi(m+1)−Xj(m+1)∣∣
如果 R m + 1 ( i ) {R_{m + 1}}(i) Rm+1(i)>> R m ( i ) {R_{m}}(i) Rm(i),則為虛假近鄰點，定義比值：
R ( i ) = R(i)= R(i)= [ R m + 1 ( i ) ] 2 − [ R m ( i ) ] 2 [ R m ( i ) ] 2 \sqrt {\frac{{{{\left[ {{R_{m + 1}}(i)} \right]}^2} - {{\left[ {{R_m}(i)} \right]}^2}}}{{{{\left[ {{R_m}(i)} \right]}^2}}}} [Rm(i)]2[Rm+1(i)]2−[Rm(i)]2
若 R ( i ) > R 0 R(i)>R_0 R(i)>R0,則稱 X j X_j Xj為 X i X_i Xi的假近鄰點， R 0 R_0 R0為閾值通常取大於10.計算該m下虛假近鄰點佔點比例，直到虛假近鄰點百分比很小或不隨m增大而減少時，此時的m即為所需嵌入維數。
#計算GP演算法的嵌入維數(假近鄰演算法)
def get_m(imf, tau):
N=2000
if (len(imf) != N):
print('請輸入指定的資料長度！') # N為指定資料長度
return 0, 0
elif (isinstance(imf, np.ndarray) != True):
print('資料格式錯誤！')
return 0, 0
else:
m_max = 10
P_m_all = [] # m_max-1個假近鄰點百分率
for m in range(1, m_max + 1):
i_num = N - (m - 1) * tau
kj_m = np.zeros((i_num, m)) # m維重構相空間
for i in range(i_num):
for j in range(m):
kj_m[i][j] = imf[i + j * tau]
if (m > 1):
index = np.argsort(Dist_m)
a_m = 0 # 最近鄰點數
for i in range(i_num):
temp = 0
for h in range(i_num):
temp = index[i][h]
if (Dist_m[i][temp] > 0):
break
D = np.linalg.norm(kj_m[i] - kj_m[temp])
D = np.sqrt((D * D) / (Dist_m[i][temp] * Dist_m[i][temp]) - 1)
if (D > 10):
a_m += 1
P_m_all.append(a_m / i_num)
i_num_m = i_num - tau
Dist_m = np.zeros((i_num_m, i_num_m)) # 兩向量之間的距離
for i in range(i_num_m):
for k in range(i_num_m):
Dist_m[i][k] = np.linalg.norm(kj_m[i] - kj_m[k])
P_m_all = np.array(P_m_all)
m_all = np.arange(1, m_max)
return m_all, P_m_all
三連、三連、三連
