互信息法确定时间序列最佳时延

1 代码实现

最近需要实现对时间序列的相空间重构，参考ChatGPT与相关论文，实现了基于互信息法确定时间序列最佳时延的程序，代码如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as pltN_ft = 1000def delay_time(data, max_delay=10):# 1. 计算自信息和联合信息ts = pd.Series(data)delays = range(1, max_delay+1)# ics = [ts.autocorr(d) for d in delays]   # 自信息ics = [ts.shift(d).autocorr() for d in delays]jcs = []for d in delays:jcs.append(ts.corr(ts.shift(d)))     # 联合信息# 2. 计算互信息print(ics)print(jcs)mis = []for jc, ic in zip(jcs, ics):print(jc, ic)mi = -0.5*np.log(1-jc**2)+0.5*np.log(1-ic**2) # 互信息print(mi)mis.append(mi)# 3. 找到第一个局部极小值并返回其对应的时延diffs = np.diff(mis)print(diffs)i = np.where(diffs > 0)[0][0]delay = delays[i]# 可视化互信息函数plt.plot(delays[0:], mis, 'bo-')plt.xlabel('Delay(τ)')plt.ylabel('Mutual Information(I(τ))')plt.grid(axis='x')plt.grid(axis='y')plt.axvline(x=delay, color='r', linestyle='--')plt.show()return delayt = []
f1 = 25
f2 = 30
for i in range(N_ft):t.append(i * 0.001)
t = np.array(t)
# yu = np.ones(M * N)
AEall = np.sin(t * 2 * np.pi * f1) + np.sin(t * 2 * np.pi * f2)  #在这里直接改信号delay = delay_time(AEall, max_delay=30)
print('Delay time:', delay)

运行结果如图所示：
延时-互信息关系曲线根据论文《混沌时间序列预测研究及应用》，寻找第一个局部极小值点确定为最佳时延，即该序列最佳时延为9.

2 相关说明

程序中相关公式为chatGPT提供，其正确性可能有待进一步确定。使用过程中若有问题，欢迎大家一起讨论！另外，在完成最佳时延的确定后，需要完成最佳嵌入维数的确定，可参考博客，这里，对其中GP算法的实现作出略微修改，代码如下：

import numpy as np
from scipy.fftpack import fft
from scipy import fftpack
import matplotlib.pyplot as plt
import pandas as pdN_ft = 1000# GP算法求关联维数(时频域特征)
def GP(imf, tau):if (len(imf) != N_ft):print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框returnelif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')returnelse:m_max = 10  # 最大嵌入维数ss = 50  # r的步长fig = plt.figure(1)for m in range(1, m_max + 1):i_num = N_ft - (m - 1) * taukj_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 = Delif (D > 0 and D < dist_min):dist_min = DDist_m[i][k] = Ddr = (dist_max - dist_min) / (ss - 1)  # r的间距r_m = []Cr_m = []for r_index in range(ss):r = dist_min + r_index * drr_m.append(r)Temp = np.heaviside(r - Dist_m, 1)for i in range(i_num):Temp[i][i] = 0Cr_m.append(np.sum(Temp))r_m = np.log(np.array((r_m)))print(r_m)Cr_m = np.log(np.array((Cr_m)) / (i_num * (i_num - 1)))print(Cr_m)plt.plot(r_m, Cr_m, label = str(m))plt.xlabel('ln(r)', fontsize=18)plt.ylabel('ln(C)', fontsize=18)plt.xticks(fontsize=16)plt.yticks(fontsize=16)plt.legend(fontsize=15)plt.show()if __name__=='__main__':# 检验关联维数程序t = []f1 = 25f2 = 30for i in range(N_ft):t.append(i * 0.001)t = np.array(t)# yu = np.ones(M * N)AEall = np.sin(t * 2 * np.pi * f1) + np.sin(t * 2 * np.pi * f2)  #在这里直接改信号GP(AEall, 1)