相空間重構求關聯維數——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 τ為時間延遲， τ \tau τ=k Δ t {\rm{\Delta }}t Δt，其中k為整數，為采樣時間間隔；m為嵌入維數；i=1，2，?，N；N為重構后向量的個數， N = n ? ( m ? 1 ) τ N = n - (m - 1)\tau N=n?(m?1)τ，
重構相空間關聯維數為：

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→0lim?lnrlncr??

c r = 1 N 2 {c_r} = \frac{1}{{{N^2}}} cr?=N21? ∑ ∑ H \sum\sum H ∑∑H ( r ? ∣ ∣ x j ? x k ∣ ∣ ) \left( {r - ||{x_j} - {x_k}||} \right) (r?∣∣xj??xk?∣∣)

式中：j≠k；r為m維超球半徑；H為Heaviside函式，

def GP(imf,tau):            #GP演算法求關聯維數
    N=2000
    if (len(imf) != N):
        print('請輸入指定的資料長度！')   # N為指定資料長度
        return
    elif (isinstance(imf, np.ndarray) != True):
        print('資料格式錯誤！')
        return
    else:
        m_max=10                  #最大嵌入維數
        ss=50                     #r的步長
        fig=plt.figure()
        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]
            dist_min, dist_max = np.linalg.norm(kj_m[0] - kj_m[1]), np.linalg.norm(kj_m[0] - kj_m[1])
            Dist_m = np.zeros((i_num, i_num))  # 兩向量之間的距離
            for i in range(i_num):
                for k in range(i_num):
                    D= np.linalg.norm(kj_m[i] - kj_m[k])
                    if(D>dist_max):
                        dist_max=D
                    elif(D>0 and D<dist_min):
                        dist_min=D
                    Dist_m[i][k] = D
            dr=(dist_max-dist_min)/(ss-1)           #r的間距
            r_m=[]
            Cr_m=[]
            for r_index in range(ss):
                r=dist_min+r_index*dr
                r_m.append(r)
                Temp=np.heaviside(r-Dist_m,1)
                for i in range(i_num):
                    Temp[i][i]=0
                Cr_m.append(np.sum(Temp))
            r_m=np.log(np.array((r_m)))
            Cr_m=np.log(np.array((Cr_m))/(i_num*(i_num-1)))
            plt.plot(r_m,Cr_m)
        plt.show()

自相關法確定 τ \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

三連、三連、三連